2.1 The Koopman operator

We consider nonlinear dynamical systems of the form

(2.1) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\boldsymbol{u}=\boldsymbol{F}(\boldsymbol{u}),\end{eqnarray}$$

with the time-forward map $\boldsymbol{f}^{t}(\boldsymbol{u})=\boldsymbol{u}+\int _{0}^{t}\boldsymbol{F}(\boldsymbol{u})\,\text{d}t^{\prime }$ . The Koopman operator, $\mathscr{K}^{t}$ , is an infinite-dimensional linear operator that propagates functionals $\unicode[STIX]{x1D713}$ of the state vector – or ‘observables’ – forward in time (Koopman Reference Koopman1931; Mezić Reference Mezić2005) along a trajectory of (2.1),

(2.2) $$\begin{eqnarray}\mathscr{K}^{t}\unicode[STIX]{x1D713}(\boldsymbol{u}):=\unicode[STIX]{x1D713}(\boldsymbol{f}^{t}(\boldsymbol{u})).\end{eqnarray}$$

The eigenfunctions of this linear, infinite-dimensional operator are special observables with exponential time dependence

(2.3) $$\begin{eqnarray}\mathscr{K}^{t}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}(\boldsymbol{u})=\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}(\boldsymbol{u})\text{e}^{\unicode[STIX]{x1D706}t}.\end{eqnarray}$$

It then follows that Koopman eigenfunctions can be computed by the relation

(2.4) $$\begin{eqnarray}\text{d}_{t}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}(\boldsymbol{u})=\boldsymbol{F}(\boldsymbol{u})\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\boldsymbol{u}}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}(\boldsymbol{u})=\unicode[STIX]{x1D706}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}(\boldsymbol{u}).\end{eqnarray}$$

The default assumption is then that the Koopman eigenfunctions can be used to expand a vector of observables,

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D74D}(\boldsymbol{u})=\mathop{\sum }_{n}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}_{n}}(\boldsymbol{u})\hat{\unicode[STIX]{x1D74D}}_{n},\end{eqnarray}$$

where the coefficients $\hat{\unicode[STIX]{x1D74D}}_{n}$ are called Koopman modes (Rowley et al. Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009). A common choice is to consider a spatially varying functional so that the vector $\unicode[STIX]{x1D74D}(\boldsymbol{u})$ is just that functional evaluated over a discretisation of space. In this case, the nonlinear evolution of $\unicode[STIX]{x1D74D}(\boldsymbol{u})$ is expressed as a superposition of fixed spatial structures (the Koopman modes) with an exponential time dependence (through the Koopman eigenfunctions), giving the appearance of linearity. Representations like (2.5) as applied to turbulent flows typically feature an additional term associated with a continuous spectrum of eigenvalues (Mezić Reference Mezić2005, Reference Mezić2013). In these applications the Koopman eigenvalues are all purely imaginary, and the discrete terms appearing in (2.5) consist of a long-time average (eigenvalue $\unicode[STIX]{x1D706}_{0}=0$ ) and ‘harmonic averages’ at discrete frequencies $\unicode[STIX]{x1D706}=\text{i}\unicode[STIX]{x1D714}$ . However, individual exact coherent structures have their own associated Koopman eigenfunctions, and our interest is in building local decompositions around these solutions (which may be embedded in a turbulent attractor). It is at present unclear how such an expansion might be constructed for a turbulent trajectory that visits multiple simple invariant solutions. For example, in a system with multiple equilibria, one would expect each fixed point to correspond to a neutral eigenfunction of the Koopman operator leading to a degeneracy of the $\unicode[STIX]{x1D706}=0$ eigenvalue.

2.2 Stuart–Landau equation

For a simple example with multiple equilibria, we first revisit the problem considered in Bagheri (Reference Bagheri2013) – an analytical derivation of the Koopman decomposition for solutions to a Stuart–Landau equation,

(2.6) $$\begin{eqnarray}\frac{\text{d}A}{\text{d}t}=a_{0}A-a_{1}A|A|^{2}.\end{eqnarray}$$

Following Bagheri (Reference Bagheri2013), we write the complex amplitude in polar coordinate, $A(t)=r(t)\text{exp}[\text{i}\unicode[STIX]{x1D703}(t)]$ . In our analysis we neglect the dependence on $\unicode[STIX]{x1D703}(t)$ by assuming $a_{0}:=\unicode[STIX]{x1D707}\in \mathbb{R}$ and without loss of generality $a_{1}:=1$ to focus solely on the evolution of amplitude variable $r(t)$ . The angular dependence is straightforward to incorporate and its inclusion only complicates the presentation. The evolution equation for $r(t)$ is

(2.7) $$\begin{eqnarray}\frac{\text{d}r}{\text{d}t}=\unicode[STIX]{x1D707}r-r^{3},\end{eqnarray}$$

which has a pitchfork bifurcation at $\unicode[STIX]{x1D707}=0$ ; for $\unicode[STIX]{x1D707}>0$ there are attractors at $r=\pm \sqrt{\unicode[STIX]{x1D707}}$ and a repeller at $r=0$ . Similar to Bagheri (Reference Bagheri2013), we consider trajectories for which $r(t=0)>0$ and $r(t\rightarrow \infty )\rightarrow \sqrt{\unicode[STIX]{x1D707}}$ and seek a Koopman representation for an observable $\unicode[STIX]{x1D713}(r)$ . However, rather than inverse engineering the Koopman eigenfunctions, eigenvalues and modes from a Fourier expansion around the limit cycle (appendix A in Bagheri Reference Bagheri2013), we identify them directly from the relationship (2.4). Since this holds universally across the dynamics and not just close to any simple invariant solution, we can construct Koopman representations for the full lifespan of the solution trajectory. Interestingly, two different non-overlapping representations emerge, one centred around the repeller ( $r=0$ ) and the other around the attractor ( $r=\sqrt{\unicode[STIX]{x1D707}}$ ), which meet at a ‘cross-over’ point where both fail simultaneously to converge. Our analysis is complementary to the study by Gaspard et al. (Reference Gaspard, Nicolis, Provata and Tasaki1995), who examined the evolution of probability densities under the dynamics described by (2.7).

2.3 Koopman mode decompositions

Assuming $\unicode[STIX]{x1D707}>0$ , equation (2.7) can be rescaled with $R:=r/\sqrt{\unicode[STIX]{x1D707}}$ and $T:=\unicode[STIX]{x1D707}t$ to

(2.8) $$\begin{eqnarray}\frac{\text{d}R}{\text{d}T}=R-R^{3}=:f(R),\end{eqnarray}$$

which has solution

(2.9) $$\begin{eqnarray}R(T;R_{0})=\frac{1}{\sqrt{1+b(R_{0})\text{e}^{-2T}}},\end{eqnarray}$$

where $b(R_{0}):=(1-R_{0}^{2})/R_{0}^{2}$ so $R(0;R_{0})=R_{0}$ . Our aim is to write the evolution of an observable, $\unicode[STIX]{x1D713}(R)$ , as an expansion in eigenfunctions of the Koopman operator. For this one-dimensional example, equation (2.4) for the Koopman eigenfunctions becomes simply

(2.10) $$\begin{eqnarray}f(R)\frac{\text{d}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}}{\text{d}R}=\unicode[STIX]{x1D706}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}},\end{eqnarray}$$

which has solution

(2.11) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}(R)=\left(\frac{R^{2}}{1-R^{2}}\right)^{\unicode[STIX]{x1D706}/2}=\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}(R_{0})\text{e}^{\unicode[STIX]{x1D706}T},\end{eqnarray}$$

where formally $\unicode[STIX]{x1D706}$ can be any complex number. The unstable eigenfunctions $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}>0}$ are singular at the attracting fixed point $R=1$ and the temporally decaying $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}<0}$ are singular at the repeller $R=0$ . In order to be able to build expansions around the two equilibria $R=0$ and $R=1$ , the eigenfunctions presumably need to be analytic in the neighbourhood of one of these equilibria which means $\unicode[STIX]{x1D706}\in \mathbb{R}$ at least.

There is another set of (weak) eigensolutions in this problem that solve

(2.12) $$\begin{eqnarray}\int f(R)\frac{\text{d}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}}{\text{d}R}w(R)\,\text{d}R=\unicode[STIX]{x1D706}\int \unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}}(R)w(R)\,\text{d}R,\end{eqnarray}$$

where $w(R)$ is a smooth weight function. The eigensolutions are

(2.13) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}_{n}}(R)=\unicode[STIX]{x1D6FF}^{(n)}\left(\frac{1-R^{2}}{R^{2}}\right),\quad \unicode[STIX]{x1D706}_{n}=2(n+1)\quad n\in \mathbb{N}\cup \{0\},\end{eqnarray}$$

and

(2.14) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}_{n}}(R)=\unicode[STIX]{x1D6FF}^{(n)}\left(\frac{R^{2}}{1-R^{2}}\right),\quad \unicode[STIX]{x1D706}_{n}=-2(n+1)\quad n\in \mathbb{N}\cup \{0\},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}^{(n)}$ is the $n\text{th}$ derivative of the delta distribution. These eigenmeasures consist of one family that exhibits exponential growth (2.13) at the stable fixed point(s) and another which is decaying (2.14) at the repeller. They were originally identified by Gaspard et al. (Reference Gaspard, Nicolis, Provata and Tasaki1995), and are important if the goal is to compute the time evolution of densities, since Koopman eigenfunctions are bi-orthogonal to eigenfunctions of the Perron–Frobenius operator (Gaspard Reference Gaspard1998). In that problem, solving the initial value problem for the density is challenging since the inclusion of (2.13) and (2.14) alongside (2.11) introduces a degeneracy. However, these eigensolutions are not required to express the evolution of observables we are interested in (namely the state itself) and we can safely ignore them.

In the absence of the eigenmeasures discussed above, there is a single one-parameter family of Koopman eigenfunctions (2.11) and a continuous spectrum of eigenvalues. A Koopman representation for a general observable would then be

(2.15) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}(R) & = & \displaystyle \int _{-\infty }^{\infty }a_{\unicode[STIX]{x1D713}}(-\unicode[STIX]{x1D706})\unicode[STIX]{x1D711}_{-\unicode[STIX]{x1D706}}(R)\,\text{d}\unicode[STIX]{x1D706},\nonumber\\ \displaystyle & = & \displaystyle \int _{-\infty }^{\infty }a_{\unicode[STIX]{x1D713}}(-\unicode[STIX]{x1D706})\unicode[STIX]{x1D711}_{-\unicode[STIX]{x1D706}}(R_{0})\text{e}^{-\unicode[STIX]{x1D706}T}\,\text{d}\unicode[STIX]{x1D706},\end{eqnarray}$$

where $a_{\unicode[STIX]{x1D713}}(-\unicode[STIX]{x1D706})$ is the Koopman mode density for the observable $\unicode[STIX]{x1D713}$ corresponding to the Koopman eigenvalue $-\unicode[STIX]{x1D706}$ . Writing the integrand in terms of $-\unicode[STIX]{x1D706}$ highlights the fact that the expression (2.15) is a bilateral Laplace transform with $\unicode[STIX]{x1D706}$ playing the role of the time-like variable and $T$ the transform variable. Setting $\unicode[STIX]{x1D713}=R$ , which is often the first observable of interest, we write

(2.16) $$\begin{eqnarray}R(T)=\int _{-\infty }^{\infty }a(-\unicode[STIX]{x1D706})\unicode[STIX]{x1D711}_{-\unicode[STIX]{x1D706}}(R_{0})\text{e}^{-\unicode[STIX]{x1D706}T}\,\text{d}\unicode[STIX]{x1D706},\end{eqnarray}$$

and then the inverse Laplace transform inversion in the complex- $T$ plane

(2.17) $$\begin{eqnarray}a(-\unicode[STIX]{x1D706})\unicode[STIX]{x1D711}_{-\unicode[STIX]{x1D706}}(R_{0})=\frac{1}{2\unicode[STIX]{x03C0}\text{i}}\int _{\unicode[STIX]{x1D6FE}-\text{i}\infty }^{\unicode[STIX]{x1D6FE}+\text{i}\infty }R(T)\text{e}^{\unicode[STIX]{x1D706}T}\,\text{d}T=\frac{1}{2\unicode[STIX]{x03C0}\text{i}}\int _{\unicode[STIX]{x1D6FE}-\text{i}\infty }^{\unicode[STIX]{x1D6FE}+\text{i}\infty }\frac{\text{e}^{\unicode[STIX]{x1D706}T}}{\sqrt{1+b(R_{0})\text{e}^{-2T}}}\,\text{d}T,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}\in \mathbb{R}$ has to be chosen such that

(2.18) $$\begin{eqnarray}\int _{-\infty }^{\infty }\text{e}^{-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D706}}|a(-\unicode[STIX]{x1D706})\unicode[STIX]{x1D711}_{-\unicode[STIX]{x1D706}}(R_{0})|\,\text{d}\unicode[STIX]{x1D706}<\infty .\end{eqnarray}$$

For unilateral Laplace transforms, this just means choosing $\unicode[STIX]{x1D6FE}$ to the right of all singularities in the complex transform variable plane. The convergence condition as $\unicode[STIX]{x1D706}\rightarrow -\infty$ is moot because $a(-\unicode[STIX]{x1D706})=0$ for all negative $\unicode[STIX]{x1D706}$ (the time-like variable). For the bilateral Laplace transform, the condition (2.18) becomes much more stringent. In particular, for $\unicode[STIX]{x1D706}\rightarrow \infty$ , $\unicode[STIX]{x1D6FE}$ must be chosen to the right of all singularities in the $T$ -plane (and the contour closed in the left-hand plane) whereas for $\unicode[STIX]{x1D706}\rightarrow -\infty$ , $\unicode[STIX]{x1D6FE}$ must be to the left of all singularities in the $T$ plane (and the contour closed in the right-hand plane). Clearly these are incompatible unless $a(\unicode[STIX]{x1D706})$ vanishes above or below some $\unicode[STIX]{x1D706}_{crit}$ . We now examine both possibilities.

The singularities of the integrand in (2.17) are the branch points

(2.19) $$\begin{eqnarray}T_{n}={\textstyle \frac{1}{2}}\text{ln}\,b+(n+{\textstyle \frac{1}{2}})\text{i}\unicode[STIX]{x03C0}\quad n\in \mathbb{Z},\end{eqnarray}$$

and require branch cuts. Considering the case of $a(\unicode[STIX]{x1D706})$ vanishing below some $\unicode[STIX]{x1D706}_{crit}$ , these branch cuts are taken out to $-\infty$ parallel to the negative $\text{Re}(T)$ axis so that the Bromwich contour is closed to the left with a large semicircle, which is indented for each of the branch cuts. The contribution on the semicircle vanishes as its radius extends to infinity provided $\unicode[STIX]{x1D706}>\unicode[STIX]{x1D706}_{crit}:=-1$ so that the integral (2.16) reads

(2.20) $$\begin{eqnarray}R_{+}(T)=\int _{-1}^{\infty }a_{+}(\unicode[STIX]{x1D706})\unicode[STIX]{x1D711}_{-\unicode[STIX]{x1D706}}(R_{0})\text{e}^{-\unicode[STIX]{x1D706}T}\,\text{d}\unicode[STIX]{x1D706},\end{eqnarray}$$

which is a representation built upon Koopman eigenfunctions with eigenvalues in $(-\infty ,1)$ . The inverse Laplace transform (2.17) is a sum over keyhole contours, $C_{n}$ , around the branch cuts

(2.21) $$\begin{eqnarray}a_{+}(\unicode[STIX]{x1D706})\unicode[STIX]{x1D711}_{-\unicode[STIX]{x1D706}}(R_{0})=-\frac{1}{2\unicode[STIX]{x03C0}\text{i}}\mathop{\sum }_{n=-\infty }^{\infty }\int _{C_{n}}\frac{\text{e}^{\unicode[STIX]{x1D706}T}}{\sqrt{1+b\text{e}^{-2T}}}\,\text{d}T.\end{eqnarray}$$

Parameterising around each keyhole contour, it can be shown that

(2.22) $$\begin{eqnarray}a_{+}(\unicode[STIX]{x1D706})\unicode[STIX]{x1D711}_{-\unicode[STIX]{x1D706}}(R_{0})=\frac{(\text{e}^{\text{i}\unicode[STIX]{x03C0}}b)^{\unicode[STIX]{x1D706}/2}}{\unicode[STIX]{x03C0}}\int _{0}^{\infty }(u^{2}+1)^{-\unicode[STIX]{x1D706}/2-1}\,\text{d}u\mathop{\sum }_{n=-\infty }^{\infty }\text{e}^{\text{i}\unicode[STIX]{x03C0}\unicode[STIX]{x1D706}n}.\end{eqnarray}$$

The infinite sum of complex exponents is the Fourier representation of a Dirac comb, $\sum _{n}\text{e}^{\text{i}\unicode[STIX]{x03C0}\unicode[STIX]{x1D706}n}=2\sum _{n}\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D706}-2n)$ . After dividing by the eigenfunction $\unicode[STIX]{x1D711}_{-\unicode[STIX]{x1D706}}$ , the Koopman mode density is found to be

(2.23) $$\begin{eqnarray}a_{+}(\unicode[STIX]{x1D706})=\frac{2\text{e}^{\text{i}\unicode[STIX]{x03C0}\unicode[STIX]{x1D706}/2}}{\unicode[STIX]{x03C0}}\underbrace{\int _{0}^{\infty }(u^{2}+1)^{-\unicode[STIX]{x1D706}/2-1}\text{d}u}_{I(\unicode[STIX]{x1D706})}\mathop{\sum }_{n=-\infty }^{\infty }\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D706}-2n).\end{eqnarray}$$

Since $\unicode[STIX]{x1D706}>-1$ and $\unicode[STIX]{x1D706}$ is even, only Koopman eigenfunctions not associated with exponential growth are included in the representation. Moreover, for discrete $\unicode[STIX]{x1D706}_{n}\in \{0,2,4,\ldots \}$ , it can be shown that $I(\unicode[STIX]{x1D706}_{n})=(1-1/2n)I(\unicode[STIX]{x1D706}_{n-1})$ , with $I(\unicode[STIX]{x1D706}_{0})=\unicode[STIX]{x03C0}/2$ , so

(2.24) $$\begin{eqnarray}I(\unicode[STIX]{x1D706}_{n})=\frac{(2n)!}{2^{2n}(n!)^{2}}\frac{\unicode[STIX]{x03C0}}{2}.\end{eqnarray}$$

Use of the Koopman mode density $a_{+}(\unicode[STIX]{x1D706})$ in (2.20) thus picks out a discrete Koopman expansion around the attractor $R=1$ as found earlier by Bagheri (Reference Bagheri2013),

(2.25) $$\begin{eqnarray}R_{+}(T)=\mathop{\sum }_{n=0}^{\infty }\underbrace{\frac{(-1)^{n}(2n)!}{2^{2n}(n!)^{2}}}_{\hat{R}_{-2n}}\unicode[STIX]{x1D711}_{-2n}(R_{0})\text{e}^{-2nT},\end{eqnarray}$$

where $\hat{R}_{-2n}$ is the Koopman mode for the observable $\unicode[STIX]{x1D713}=R$ associated with the Koopman eigenfunction $\unicode[STIX]{x1D711}_{-2n}$ . Note that the Koopman eigenvalues in this expansion, $\unicode[STIX]{x1D706}_{n}=-2n$ , are integer multiples of the stable eigenvalue $-2$ associated with the linearised dynamics around the fixed point $R=1$ . The expansion therefore agrees with the results of Mezic (Reference Mezic2017), who proved that the Koopman expansion should take this form in nonlinear systems with a single stable attracting fixed point. In this case, however, the presence of multiple equilibria complicates matters and (2.25) is not uniformly convergent. In fact, this series can be recognised as just the Taylor expansion of the exact solution (2.9)

(2.26) $$\begin{eqnarray}R(T;R_{0})=\frac{1}{\sqrt{1+{\displaystyle \frac{1-R_{0}^{2}}{R_{0}^{2}}}\text{e}^{-2T}}}\end{eqnarray}$$

in $y:=((1-R_{0}^{2})/R_{0}^{2})$ about $y=0$ valid for dynamics ‘close’ to the attracting fixed point $R=1$ (Bagheri (Reference Bagheri2013) proceeded in the opposite direction starting from the Taylor expansion to deduce the Koopman expansion). Since $((1-R(T)^{2})/R(T)^{2})=((1-R_{0}^{2})/R_{0}^{2})\text{e}^{-2T}$ , we can rewrite this as the identity

(2.27) $$\begin{eqnarray}R=\frac{1}{\sqrt{1+{\displaystyle \frac{1-R^{2}}{R^{2}}}}},\end{eqnarray}$$

which indicates that the representation (2.25) will fail to converge for $(1-R^{2})/R^{2}\geqslant 1$ or $R\leqslant 1/\sqrt{2}$ . In other words, since $R(T)$ increases monotonically with time, the representation (2.25) holds for any solution with initial condition $R_{0}>1/\sqrt{2}$ .

We now turn our attention to the other possible scenario where $a(\unicode[STIX]{x1D706})$ vanishes above some $\unicode[STIX]{x1D706}_{crit}$ . The branch cuts must now be taken out to $+\infty$ parallel to the positive $\text{Re}(T)$ axis so that the Bromwich contour is closed to the right with a large semicircle, which is indented for each of the branch cuts. The contribution on the semicircle vanishes as its radius extends to infinity provided that $\unicode[STIX]{x1D706}<\unicode[STIX]{x1D706}_{crit}:=0$ and we write

(2.28) $$\begin{eqnarray}R_{-}(T)=\int _{-\infty }^{0}a_{-}(\unicode[STIX]{x1D706})\unicode[STIX]{x1D711}_{-\unicode[STIX]{x1D706}}(R_{0})\text{e}^{-\unicode[STIX]{x1D706}T}\,\text{d}T\end{eqnarray}$$

so only Koopman eigenfunctions associated with exponential growth are included. Parameterising around the branch cuts as previously yields the Koopman mode density,

(2.29) $$\begin{eqnarray}a_{-}(\unicode[STIX]{x1D706})=\frac{2\text{e}^{-(1/2)\text{i}\unicode[STIX]{x03C0}(\unicode[STIX]{x1D706}+1)}}{\unicode[STIX]{x03C0}}\int _{0}^{\infty }(u^{2}+1)^{(\unicode[STIX]{x1D706}-1)/2}\text{d}u\mathop{\sum }_{n=-\infty }^{\infty }\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D706}+1-2n).\end{eqnarray}$$

As before, we can evaluate the integral using a recurrence relation and using (2.28) we recover another Koopman expansion

(2.30) $$\begin{eqnarray}R_{-}(T)=\mathop{\sum }_{n=0}^{\infty }\frac{(-1)^{n}(2n)!}{2^{2n}(n!)^{2}}\unicode[STIX]{x1D711}_{2n+1}(R_{0})\text{e}^{(2n+1)T}.\end{eqnarray}$$

This is the Taylor expansion in $z:={\textstyle \frac{R}{\sqrt{1-R^{2}}}}$ around $z=0$ of the exact solution

(2.31) $$\begin{eqnarray}R=\frac{z}{\sqrt{1+z^{2}}}\end{eqnarray}$$

(a simple manipulation of the identity (2.27)) which fails to converge when $z=1$ or $R(T)\geqslant 1/\sqrt{2}$ . So if $R_{0}<1/\sqrt{2}$ , this representation will hold until $R=1/\sqrt{2}$ . Beyond this point in time, the other Koopman expansion can then be used to represent the solution. So the two Koopman decompositions (2.25) and (2.30) together allow (almost) the entire nonlinear evolution to be expressed as a superposition of linear (exponential time dependence) observables. The performance of the two decompositions, truncated at a finite number of Koopman modes, is examined in figure 1. As expected, the two expansions fail as they are pushed beyond the cross-over point $R=1/\sqrt{2}$ .

Figure 1.

(a,b) True evolution (2.9) with $R_{0}=10^{-3}$ (dashed black) with Koopman approximations (2.25) and (2.30) overlaid ((a,c)  $N=2$ modes, (b,d)  $N=5$ ; red/blue correspond to the attracting/repelling expansions respectively). (c,d) The error $\unicode[STIX]{x1D700}_{N}:=|R(t)-R_{\pm }(t;N)|$ , where $N$ is the number of modes included in the expansion. Grey region identifies the cross-over point between repelling and attracting Koopman expansions.

At this point it is interesting to ask what goes wrong in attempting to build a Koopman expansion centred around another point (say, even $R=1/\sqrt{2}$ ) which is not an equilibrium. Here a connection with Carleman linearisation (Carleman Reference Carleman1932) is useful. Carleman linearisation makes a nonlinear system linear by relabelling each nonlinearity as a new dependent variable of the system. Typically, this converts a finite-dimensional nonlinear system into an infinite linear system as additional equations need to be added to describe how the new dependent variables evolve. This generically introduces further nonlinearities and the procedure mushrooms with yet more variables needing to be defined. When this linearisation procedure is carried out around a solution of the system such as an equilibrium, it produces a purely linear system as opposed to the generic affine one – i.e. the time evolution of the system is given by a linear operator. The (adjoint) eigenfunctions and eigenvalues (modulo exponentiation) of this operator would then seem to correspond with the Koopman (eigenfunctions) modes and eigenvalues of the Koopman operator. A Koopman expansion centred at this point clearly makes sense. In contrast, if the Carleman linearisation is performed around a non-solution, the resulting system is then only affine and the temporal evolution cannot be purely expressible as a sum of exponentially time varying Koopman modes: see appendix A for details for the model studied here. Thus, it would only seem to make sense to talk about Koopman expansions about simple invariant solutions or just the equilibria $R=\{0,1\}$ here for $R\geqslant 0$ .

Finally, it is worth emphasising that the breakdown of a given Koopman expansion is associated with a loss of convergence rather than any pathology in the component Koopman eigenfunctions. In fact, the Koopman eigenfunctions exist everywhere away from the fixed points $R=\{0,1\}$ . The point is just that certain subsets cannot be used in a convergent representation at a given point in the dynamics. The jump between Koopman expansions has important consequences for DMD, which we now explore.

2.4 Dynamic mode decomposition

In his examination of the transient collapse of the flow over a cylinder onto the vortex-shedding limit cycle, Bagheri (Reference Bagheri2013) noted a curious dependence of the output of DMD on the length of time window over which data were collected. If the observation window was restricted to a time interval where the velocity, $\boldsymbol{u}$ , is ‘close’ to the periodic orbit, DMD accurately reproduced the Koopman eigenvalues for the attracting expansion. However, when the window was extended to also include the ‘amplification’ region, the DMD eigenvalues appeared as discrete approximations to continuous lines of decaying eigenvalues. This phenomenon is a consequence of the jump between Koopman expansions around the fixed points of equation (2.7).

In fact, the cross-over point between Koopman expansions has a critical impact on the ability of DMD to approximate a Koopman expansion even if the DMD design is perfect in all other respects – i.e. the elements of the user-defined observable vector, $\unicode[STIX]{x1D74D}(\boldsymbol{u})$ , are a suitable basis for the Koopman eigenfunctions and a sufficiently large amount of data is examined (Williams et al. Reference Williams, Kevrekidis and Rowley2015). In particular, if the desire is to obtain a Koopman decomposition of the state variable, the DMD must be restricted to within a neighbourhood of an exact solution inside which the expansion is valid. We demonstrate this behaviour by conducting DMD of the simple one-dimensional (1-D) problem (2.7), successfully obtaining Koopman eigenvalues only when snapshot pairs are restricted to times $\{T:R(T)<1/\sqrt{2}\}$ or $\{T:R(T)>1/\sqrt{2}\}$ , but not for an overlapping interval.

We generate snapshots $\{R_{i}\}$ of the trajectory reported in figure 1 with a spacing $\unicode[STIX]{x0394}t_{s}=0.1$ on the interval $t\in [0,20]$ . For the DMD, we use a snapshot spacing $\unicode[STIX]{x1D6FF}t=1$ , so in total we have available $M_{max}=190$ snapshot pairs ( $t=0$ to $19$ step $0.1$ mapped to $t=1$ to $20$ step $0.1$ ). The observable vector for the DMD is made up of polynomials in $R$ , $\unicode[STIX]{x1D74D}(R)=(R,R^{2},R^{3},R^{4})$ . The DMD reported here is slightly unusual in the sense that $N$ , the dimension of $\unicode[STIX]{x1D74D}$ , is much less than $M$ , the number of snapshots. Typically $N\gg M$ in fluid mechanics, and it will be shown in § 3 that analogous behaviour to that found in this 1-D problem occurs along heteroclinic connections between equilibria of the Navier–Stokes equations.

The DMD methodology is essentially as specified in Tu et al. (Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2014), although the inclusion of polynomials of the state in $\unicode[STIX]{x1D74D}$ makes the current problem an example of extended dynamic mode decomposition (Williams et al. Reference Williams, Kevrekidis and Rowley2015). Given a matrix of snapshots,

(2.32) $$\begin{eqnarray}\unicode[STIX]{x1D733}^{t}=\left[\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D74D}(R(t_{i})) & \unicode[STIX]{x1D74D}(R(t_{j})) & \cdots \end{array}\right],\end{eqnarray}$$

and a corresponding matrix with the observables now evaluated $\unicode[STIX]{x1D6FF}t$ later,

(2.33) $$\begin{eqnarray}\unicode[STIX]{x1D733}^{t+\unicode[STIX]{x1D6FF}t}=\left[\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D74D}(R(t_{i}+\unicode[STIX]{x1D6FF}t)) & \unicode[STIX]{x1D74D}(R(t_{j}+\unicode[STIX]{x1D6FF}t)) & \cdots \end{array}\right],\end{eqnarray}$$

the DMD operator $\hat{\boldsymbol{K}}$ is the linear operator which best maps between corresponding snapshot pairs (in a least-squares sense),

(2.34) $$\begin{eqnarray}\hat{\boldsymbol{K}}:=\unicode[STIX]{x1D733}^{t+\unicode[STIX]{x1D6FF}t}(\unicode[STIX]{x1D733}^{t})^{+},\end{eqnarray}$$

where the $+$ superscript indicates a pseudo (Moore–Penrose) inverse. Note that the snapshot times, $\{t_{i}\}$ , do not need to be sequential, and are drawn randomly from within the time interval of interest. As described in Rowley & Dawson (Reference Rowley and Dawson2017), the right eigenvectors of the DMD operator, $\hat{\unicode[STIX]{x1D74D}}_{j}=(\hat{R}_{j},\hat{R}_{j}^{2},\hat{R}_{j}^{3},\hat{R}_{j}^{4})^{\text{T}}$ , approximate Koopman modes, while the left eigenvectors, $\boldsymbol{w}_{j}$ , can be used to find the Koopman eigenfunctions,

(2.35) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{j}(R)=\boldsymbol{w}_{j}^{H}\unicode[STIX]{x1D74D}(R),\end{eqnarray}$$

under the assumptions that (i) the elements of $\unicode[STIX]{x1D74D}$ constitute a suitable basis for the eigenfunctions and (ii) sufficient data have been collected such that $\boldsymbol{w}\in \text{range}(\unicode[STIX]{x1D733}^{t})$ .

Figure 2.

Eigenvalues obtained from DMD on the evolution shown in figure 1 with $\unicode[STIX]{x1D74D}(R)=(R,R^{2},R^{3},R^{4})^{\text{T}}$ . $M=20$ snapshots and $\unicode[STIX]{x1D6FF}t=1$ obtained on $t\in [\!0,5\!)$ , (b) $t\in [5,15)$ and (c) $t\in [10,20)$ .

Eigenvalues from three DMDs are reported in figure 2. Each calculation was performed on snapshot pairs extracted from a different time window. When the time window is limited to the repelling region, DMD yields eigenvalues $\unicode[STIX]{x1D706}_{n}=n$ (while the expansion for $R$ around the repeller requires only odd integers, the inclusion of powers of $R$ in the observable means a larger set, $n\in \mathbb{N}$ , are uncovered: odd integers can sum to be even). When the time window lies within the region of validity for the attracting expansion, DMD finds the attractor eigenvalues, $\unicode[STIX]{x1D706}_{n}=-2n$ (sums of even integers remain even). On the other hand, DMD on snapshots from a time window which overlaps both expansion regions is unable to find eigenvalues for either expansion. While these results were obtained for relatively short observation windows ( $T_{w}\in \{5,10\}$ ), analogous results are obtained for longer windows. For example, even windows beginning at $T=0$ fail to identify the repelling expansion if the cross-over point is included (not shown).

Figure 3.

Comparison of Koopman eigenfunctions (colours) with eigenfunctions extracted from the repelling and attracting DMDs (black) reported in figure 2, alongside corresponding DMD modes $\hat{\unicode[STIX]{x1D74D}}$ . Note that the $l\text{th}$ component of the $j\text{th}$ DMD mode represents the DMD approximation to $\hat{R}_{j}^{l}$ , the $j\text{th}$ Koopman mode in the Koopman decomposition of $R^{l}$ .

The performance of the DMD can be assessed in more detail by comparing the predicted Koopman eigenfunctions and modes to those derived in § 2. In figure 3 the DMD approximations to the Koopman eigenfunctions are reported for both the ‘repelling’ and ‘attracting’ windows. In both cases, the DMD algorithm is able to build a locally valid approximation to the true Koopman eigenfunction from the polynomials $R^{m}$ in the observable vector $\unicode[STIX]{x1D74D}$ . These locally valid expansions break down as the cross-over point, $R=1/\sqrt{2}$ , is approached. The correspondence between DMD and Koopman also gets progressively worse for the higher-order eigenfunctions – a consequence of the limited number of polynomials in $\unicode[STIX]{x1D74D}$ .

The DMD modes reported in figure 3 should be interpreted in the following way: the $l\text{th}$ component of the DMD mode alongside eigenfunction $\unicode[STIX]{x1D711}_{j}$ is the DMD approximation to the $j\text{th}$ Koopman mode in an expansion of $R^{l}$ i.e. $\hat{R}_{j}^{l}$ . So, for example, component $\hat{\unicode[STIX]{x1D713}}_{l=2}$ alongside eigenfunction $\unicode[STIX]{x1D711}_{4}$ (bottom left corner of figure 3) is the DMD approximation to Koopman mode $\hat{R}_{4}^{2}$ in the expansion $R^{2}=\sum _{m\in \mathbb{N}}\unicode[STIX]{x1D711}_{2m}(R)\hat{R}_{2m}^{2}$ . The DMD approximations to the Koopman modes reported in figure 3 are consistent with the analytical expansions derived in § 2. For example, the repeller decomposition (2.30) indicates that the Koopman eigenvalues required to advance $R(T)$ are the odd integers. The DMD identifies a broader set of Koopman eigenvalues, $\unicode[STIX]{x1D706}\in \mathbb{N}$ , than those needed for $R$ alone, but correctly finds that the Koopman mode $\hat{R}_{2}=0$ while picking up the contributions $\hat{R}_{1}$ and $\hat{R}_{3}$ (the DMD mode for $\hat{R}_{4}$ is non-zero but small – DMD with higher-order polynomials included in $\unicode[STIX]{x1D74D}$ can eliminate this error).

The first non-zero Koopman eigenvalue in both expansions is the growth/decay rate associated with the locally linear dynamics around the repelling and attracting equilibria respectively. The higher-order terms in the Koopman decompositions allow us to propagate observables (in particular the state variable itself) beyond these linear subspaces, and we have demonstrated here that DMD is a robust method for finding these contributions provided that the observation window is contained within a particular ‘expansion region’. In the remainder of this paper we will show how similar behaviour is observed along heteroclinic connections between equilibria of the Navier–Stokes equations, and that DMD can successfully identify modes associated with repelling and attracting expansions along their unstable and stable manifolds, respectively.

3.1 Heteroclinic connection between Nagata solutions

We consider a Nagata (Reference Nagata1990) box of size $(L_{x},L_{y},L_{z})=(5\unicode[STIX]{x03C0}/2,4\unicode[STIX]{x03C0}/3,2)$ initially at $Re=135$ . In this configuration the Navier–Stokes equations support three equilibrium solutions – the constant shear solution $\boldsymbol{u}_{C}$ and the Nagata lower- and upper-branch solutions, $\boldsymbol{u}_{LB}$ and $\boldsymbol{u}_{UB}$ respectively (Nagata Reference Nagata1990). These two solutions are born out of a saddle-node bifurcation at around $Re\sim 125$ (Nagata Reference Nagata1990) for this box. At $Re=135$ both $\boldsymbol{u}_{C}$ and $\boldsymbol{u}_{UB}$ are stable while $\boldsymbol{u}_{LB}$ is the unstable edge state on the dividing manifold between their respective basins of attraction. We compute $\boldsymbol{u}_{LB}$ in the ‘GHC’ box at $Re=400$ by using a snapshot of (transient) turbulence as a guess in the Newton–GMRES–hookstep algorithm described above. This solution is then continued down to the target Reynolds number and target box size. A finite-time approximation to the heteroclinic connection between $\boldsymbol{u}_{LB}$ and $\boldsymbol{u}_{UB}$ is then obtained in the following manner: (i) velocity snapshots are generated along a short trajectory $t\in [0,50]$ with the initial condition $\boldsymbol{u}_{0}=(1+\unicode[STIX]{x1D700})\boldsymbol{u}_{LB}-\unicode[STIX]{x1D700}\boldsymbol{u}_{C}$ , where $\unicode[STIX]{x1D700}=10^{-6}$ ; (ii) the unstable eigenfunction, $\hat{\boldsymbol{u}}_{1}$ , is extracted from this trajectory using DMD; (iii) the new initial condition $\boldsymbol{u}_{0}^{\prime }=\boldsymbol{u}_{LB}+\unicode[STIX]{x1D6FF}\hat{\boldsymbol{u}}_{1}$ , where $\unicode[STIX]{x1D6FF}|\boldsymbol{u}_{LB}|=10^{-8}$ is then used to compute a more accurate approximation to the heteroclinic connection. At $Re=135$ , the first initial condition $\boldsymbol{u}_{0}$ is actually sufficient to obtain a good approximation to the heteroclinic connection since the upper-branch solution is stable. However, at higher Reynolds numbers $\boldsymbol{u}_{UB}$ becomes unstable, and the initial condition described in (iii) can generate trajectories which still spend some time in its vicinity before being flung out along its unstable manifold.

The energy production,

(3.1) $$\begin{eqnarray}I^{\prime }:=\frac{1}{2L_{x}L_{y}}\int _{0}^{L_{x}}\int _{0}^{L_{y}}\frac{1}{Re}\frac{\unicode[STIX]{x2202}u}{\unicode[STIX]{x2202}z}\bigg|_{z=\pm 1}\,\text{d}x\,\text{d}y\end{eqnarray}$$

per unit area and arclength and normalised by its value in laminar flow ( $I:=I^{\prime }/I_{lam}^{\prime }$ ) is computed along the heteroclinic connection and is reported as a function of time in figure 4 to highlight the qualitative similarity with the evolution $R(T)$ in the model problem of § 2. Similar to the behaviour near to the origin there, the Nagata lower-branch solution is a repeller with a single unstable direction. However, the stable subspace around the attractor, $\boldsymbol{u}_{UB}$ , is more complex (four-dimensional), as described below. While analytical construction of Koopman decompositions around these fixed points is not possible here, we employ DMD to identify repelling and attracting Koopman expansions.

Figure 4.

(a) Energy production along finite-time approximation to the $\boldsymbol{u}_{LB}\rightarrow \boldsymbol{u}_{UB}$ heteroclinic connection at $Re=135$ . (b) Arclength along the heteroclinic connection measured as distance from $\boldsymbol{u}_{LB}$ .

Figure 5.

(a) Imaginary and (b) real components of DMD eigenvalues obtained for a time window of length $T_{w}=100$ passed through the trajectory shown in figure 4 (the variable $t_{F}$ is the ‘final’ time of each time window). Each individual calculation is performed on $M=25$ snapshots pairs separated by $\unicode[STIX]{x1D6FF}t=2$ selected randomly within each time window. Red and blue colouring indicate whether the behaviour is classified as locally ‘repelling’ or ‘attracting’ respectively.

For the DMD, 1000 snapshots of the full velocity field, separated by $\unicode[STIX]{x0394}t_{s}=2$ , are stored for the trajectory in figure 4 and the observable vector is

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D74D}(\boldsymbol{u})=\boldsymbol{u}-\boldsymbol{u}_{C}.\end{eqnarray}$$

Initially, we pass a fixed time window of width $T_{w}=100$ along the heteroclinic connection, performing many DMD calculations with the results collated in figure 5. Each individual calculation is performed with $M=25$ snapshot pairs separated by $\unicode[STIX]{x1D6FF}t=2$ extracted randomly from within the interval of interest. Initially, and as anticipated for a trajectory repelled from the edge, the DMD identifies a single unstable eigenvalue $\unicode[STIX]{x1D706}_{1}\approx 0.02$ associated with the unstable linear subspace about $\boldsymbol{u}_{LB}$ . As the time window is passed along the heteroclinic connection, further unstable eigenvalues $\unicode[STIX]{x1D706}_{n}=n\unicode[STIX]{x1D706}_{1}$ are uncovered. This suggests that the DMD algorithm is identifying Koopman eigenfunctions in the same family as $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}_{1}}(\boldsymbol{u})$ , i.e. $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}_{n}}(\boldsymbol{u})=\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}_{1}}^{n}(\boldsymbol{u})$ (higher harmonics of the primary instability) in analogy to the model problem considered in § 2.

The eigenvalues for one particular DMD calculation inside this ‘growing’ region are reported in figure 6, and the corresponding DMD modes, $\{\boldsymbol{v}_{n}\}$ , are shown in figure 7. The neutral DMD mode is Nagata’s lower-branch solution, and the first growing mode is localised at the critical layer where $\boldsymbol{u}_{LB}\boldsymbol{\cdot }\hat{\boldsymbol{x}}=0$ . The higher-order modes are qualitatively similar to the first.

Figure 6.

Eigenvalues obtained from a DMD calculation on $t\in [400,500]$ from the trajectory shown in figure 4 (cf. figure 5). The number of snapshots is $M=25$ and $\unicode[STIX]{x1D6FF}t=2$ . The dashed red lines identify integer multiples of the first unstable mode, $\unicode[STIX]{x1D706}_{1}\approx 0.020$ .

As the DMD window is pushed further along the heteroclinic connection, pairs of (unstable) complex-conjugate eigenvalues emerge (beyond $t_{F}\sim 700$ in figure 5) with growth rates/frequencies that are inconsistent from calculation to calculation. However, beyond $t_{F}\sim 800$ a new picture emerges, and DMD identifies a variety of decaying modes that are consistent over many time windows. This behaviour is analogous to the cross-over to the ‘attracting’ expansion observed in the Stuart–Landau equation in § 2. Furthermore, the fact that the time interval where the DMD output is inconsistent (highlighted in purple in figure 5) is roughly equal to the length of the DMD time window itself, $T_{w}=100$ , hints that there may also be a single cross-over point between the two decompositions identified in the DMD rather than a finite patch of state space where neither expansion holds. At late times the DMD identifies a single complex-conjugate pair of decaying modes in addition to a neutral eigenvalue, which indicates that trajectories spiral into the upper branch.

An example eigenvalue spectrum from the ‘decaying’ region of the heteroclinic connection is reported in figure 8. There is a neutral mode which is the upper-branch (stable) equilibrium. The modes highlighted in blue also include the complex-conjugate pair of modes commented on above, $\unicode[STIX]{x1D706}_{1}^{\pm }\approx -0.017\pm 0.031\text{i}$ . The other blue eigenvalues then sit in a ‘lattice’ structure in the complex plain, as predicted for an isolated stable fixed point by Mezic (Reference Mezic2017): these eigenvalues are built from linear combinations of slowest-decaying (linear) pair, i.e. $\unicode[STIX]{x1D706}_{1}^{+}+\unicode[STIX]{x1D706}_{1}^{-}$ , $2\unicode[STIX]{x1D706}_{1}^{+}$ and $2\unicode[STIX]{x1D706}_{1}^{-}$ . If the Koopman eigenfunctions associated with the least-decaying pair are $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}_{1}}^{\pm }(\boldsymbol{u})$ , then the eigenfunctions corresponding to the higher-order modes are $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}_{1}}^{+}(\boldsymbol{u})\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}_{1}}^{-}(\boldsymbol{u})$ , $(\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}_{1}}^{+}(\boldsymbol{u}))^{2}$ and $(\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}_{1}}^{-}(\boldsymbol{u}))^{2}$ respectively. The corresponding DMD modes are reported in figure 9.

Figure 7.

DMD modes corresponding to the four eigenvalues reported in figure 6, in order of increasing growth rate from top to bottom. The visualisation of each mode shows contours of the streamwise velocity at the midplane $z=0$ (a) and contours of streamwise velocity with the streamfunction overlaid on a cross-stream plane at $x=0$  (b).

Figure 8.

Eigenvalues obtained from a DMD calculation on $t\in [1000,1100]$ from the trajectory shown in figure 4 (cf. figure 5). The number of snapshots is $M=25$ and $\unicode[STIX]{x1D6FF}t=2$ . The dashed blue lines identify integer multiples of the slowest-decaying mode, $\unicode[STIX]{x1D706}_{1}\approx -0.017+0.031\text{i}$ .

Figure 9.

DMD modes (real part shown) corresponding to the four blue eigenvalues with $\unicode[STIX]{x1D706}_{i}\geqslant 0$ reported in figure 8, in order of increasing $|\unicode[STIX]{x1D706}|$ from top to bottom. The visualisation of each mode shows contours of the streamwise velocity at the midplane $z=0$ (a) and contours of streamwise velocity with the streamfunction overlaid on a cross-stream plane at $x=0$ (b).

In addition to the family of Koopman eigenfunctions linked to the least-damped linear behaviour, there is an additional complex-conjugate pair of eigenvalues $\unicode[STIX]{x1D701}_{1}^{\pm }\approx 0.03\pm 0.10\text{i}$ highlighted in orange in figure 8. This pair of modes is consistent across the DMD calculations reported in figure 5, although it is a little difficult to distinguish in that figure due to the closeness of the decay rate to $2\unicode[STIX]{x1D706}_{1}$ . These eigenvalues also describe a decaying spiral and indicate that the stable subspace around the Nagata upper branch solution is actually four-dimensional. Their importance as the Reynolds number is increased is discussed in more detail below.

The existence of a cross-over point between the two Koopman decompositions can be explored further by using the output of the DMD calculations to construct approximations to the true trajectory. To that end, we use the two DMD calculations reported in figures 6 and 8 to construct approximations to the true heteroclinic connection. We seek a low-dimensional representation of the flow from DMD, $\boldsymbol{u}_{D}$ , by summing over a subset ${\mathcal{V}}_{\pm }$ of the DMD modes from either the repelling or attracting regions,

(3.3) $$\begin{eqnarray}\boldsymbol{u}_{D}(\boldsymbol{x},t)=u_{C}(\boldsymbol{x})+\mathop{\sum }_{\unicode[STIX]{x1D706}_{j}\in {\mathcal{V}}}a_{j}\boldsymbol{v}_{j}(\boldsymbol{x})\text{e}^{\unicode[STIX]{x1D706}_{j}t}.\end{eqnarray}$$

For each expansion, the modes in ${\mathcal{V}}_{\pm }$ are exactly those reported in figures 6 and 8. For the repelling expansion, this includes the neutral mode and the three unstable eigenvalues. For the attracting expansion, the neutral mode, the five stable (blue) modes associated with the slowest-decaying spiral and the complex-conjugate (orange) pair of modes associated with the second, more rapidly decaying spiral are included.

Figure 10.

Error in the DMD approximation(s) (3.3) versus the true evolution, $\unicode[STIX]{x1D700}:=\Vert \boldsymbol{u}_{D}-\boldsymbol{u}\Vert /\Vert \boldsymbol{u}\Vert$ . Red and blue lines identify approximations to the attracting and repelling expansions respectively, the grey regions identify windows where the DMD calculations and fitting were performed.

The amplitudes, $\{a_{j}\}$ , assigned to the DMD modes are determined by a least-squares fit to the true trajectory within the DMD time window. Taking $M$ equally spaced snapshots along the fitting window separated by a time $\unicode[STIX]{x1D6FF}t$ , the function to be minimised is

(3.4) $$\begin{eqnarray}J(\boldsymbol{a}):=\frac{1}{M}\mathop{\sum }_{m=0}^{M-1}\left|\unicode[STIX]{x1D74D}(\boldsymbol{u}(\boldsymbol{x},m\unicode[STIX]{x1D6FF}t))-\mathop{\sum }_{j}a_{j}\boldsymbol{v}_{j}\text{e}^{m\unicode[STIX]{x1D706}_{j}\unicode[STIX]{x1D6FF}t}\right|^{2}.\end{eqnarray}$$

The solution to the least-squares problem for $\boldsymbol{a}$ is then

(3.5) $$\begin{eqnarray}\hat{\boldsymbol{a}}=\bigg(\mathop{\sum }_{m}(\unicode[STIX]{x1D726}^{H})^{m}\unicode[STIX]{x1D651}^{H}\unicode[STIX]{x1D651}\unicode[STIX]{x1D726}^{m}\bigg)^{-1}\mathop{\sum }_{m}(\unicode[STIX]{x1D726}^{H})^{m}\unicode[STIX]{x1D651}^{H}\unicode[STIX]{x1D74D}(\boldsymbol{u}(\boldsymbol{x},m\unicode[STIX]{x1D6FF}t)),\end{eqnarray}$$

where $\unicode[STIX]{x1D726}$ is a diagonal matrix where the $i\text{th}$ entry is $\text{e}^{\unicode[STIX]{x1D706}_{i}\unicode[STIX]{x1D6FF}t}$ and $\unicode[STIX]{x1D651}$ is a matrix whose $i\text{th}$ column is the $i\text{th}$ DMD mode.

The error between the repelling and attracting approximations and the true solution is reported in figure 10. Unsurprisingly, the error in each case is smallest in the fitting windows themselves. The error also remains vanishingly small as each expansion is pushed towards its respective equilibrium solution, but both expansions blow up as they are pushed beyond an apparent cross-over point at $t\approx 750$ .

In addition to computing the error against the true trajectory, an alternative way of assessing the output of the DMD calculations in connection to the Koopman operator is to examine the numerical approximation to the Koopman eigenfunctions. As described in § 2, these objects are obtained from the left eigenvectors of the DMD operator, $\{\boldsymbol{w}_{j}\}$ as follows: (Rowley & Dawson Reference Rowley and Dawson2017),

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D711}_{j}(\boldsymbol{u})=\boldsymbol{w}_{j}^{H}\unicode[STIX]{x1D74D}(\boldsymbol{u}).\end{eqnarray}$$

The ‘performance’ of DMD can be examined by evaluating this inner product for points on the trajectory beyond the DMD time window, which we do in figure 11. It is clear that the DMD calculations on the relatively short time window $T_{w}=100$ have been able to accurately build locally valid representations of the Koopman eigenfunctions which remain reasonably accurate for two to three hundred advective time units beyond the observation window. These local approximations become increasingly poor around the ‘cross-over point’ inferred from earlier figures (e.g. figures 5 and 10), a behaviour which again is analogous to DMD of the Stuart–Landau equation (e.g. figure 3).

Figure 11.

DMD approximations to Koopman eigenfunctions, including their time evolution, obtained in both the repelling (a – also see figures 6 and 7) and attracting (b – also see figures 8 and 9) regions. The DMD time window is highlighted in grey, and the dashed lines identify temporal behaviour ${\sim}\text{e}^{\unicode[STIX]{x1D706}_{j}t}$ . Note all repelling eigenfunctions are normalised to unity at $t=400$ ; the attracting eigenfunctions are normalised to unity at $t=1100$ .

3.2 Higher Reynolds numbers

In both the one-dimensional Stuart–Landau equation (§ 2) and the example discussed above at $Re=135$ , there are only two fixed points: repelling and attracting equilibria. This results in a pair of Koopman expansions that extend beyond the respective repelling/attracting linear subspaces to a cross-over point in state space. Here, we increase the Reynolds number in the Nagata box to examine the consequences for Koopman decompositions and DMD when the structure of state space is complicated by the presence of additional invariant sets. The motivation here is to explore the possibility of applying DMD to turbulent trajectories as a method of locating nearby coherent structures and their associated Koopman mode expansions.

The eigenvalue spectra obtained in figures 5 and 8 indicate that the stable subspace around the upper branch is four-dimensional (two orthogonal spirals). The DMD also revealed the higher-order Koopman eigenvalues required to propagate $\boldsymbol{u}$ beyond the linear subspace. The second, more rapidly decaying spiral (orange squares in figure 8) becomes increasingly dominant in the dynamics as the Reynolds number is increased. This behaviour is apparent in figure 12, where we report DMD eigenvalues from trajectories very close to the upper branch $\boldsymbol{u}_{UB}$ at $Re\in \{140,150\}$ .

At $Re=140$ the eigenvalues associated with the second spiral in the linear subspace around $\boldsymbol{u}_{UB}$ ( $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D701}_{1}^{\pm }\approx =-0.015\pm 0.091\text{i}$ ; orange in figure 12), which had only a weak effect on the dynamics at $Re=135$ , have become the most slowly decaying to dominate the linearised dynamics. In addition, higher-order Koopman eigenvalues in the same family as $\unicode[STIX]{x1D701}_{1}^{\pm }$ are also obtained from the DMD (e.g. $2\unicode[STIX]{x1D701}_{1}$ associated with $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D701}_{1}}^{2}(\boldsymbol{u})$ ) and are highlighted with dashed lines. The decay rate of the first spiral (blue triangles) has approximately doubled. Note that, in addition to the two families of Koopman eigenvalues associated with the dynamics in the linear subspace, there are also eigenvalues which are connected to products of Koopman eigenfunctions from these families. For example, in figure 12 the green diamonds identify eigenvalues $\unicode[STIX]{x1D706}_{1}^{\pm }+\unicode[STIX]{x1D701}_{1}^{\pm }$ associated with eigenfunctions $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D706}_{1}}^{\pm }(\boldsymbol{u})\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D701}_{1}}^{\pm }(\boldsymbol{u})$ .

Figure 12.

DMD eigenvalue spectra obtained in the vicinity of $\boldsymbol{u}_{UB}$ at $Re=140$ (a) and $Re=150$ (b). At $Re=140$ , the upper-branch solution is stable and the spectrum was obtained in a similar manner to that shown in figure 8 at $Re=135$ ( $M=25$ snapshots in a time window of length $T_{w}=100$ , with $\unicode[STIX]{x1D6FF}t=2$ ). At $Re=150$ , the DMD was performed on a trajectory $\boldsymbol{f}^{t}(\tilde{\boldsymbol{u}}_{UB})$ over $t\in [0,700)$ , where $\tilde{\boldsymbol{u}}_{UB}$ is the numerical approximation to the upper-branch equilibrium converged using Newton–GMRES over a time interval $T=4$ . $M=25$ snapshots were used with spacing $\unicode[STIX]{x1D6FF}t=1$ .

At around $Re\approx 145$ the complex-conjugate pair of eigenvalues associated with the dominant spiral, $\unicode[STIX]{x1D701}_{1}^{\pm }$ , cross the imaginary axis (not shown) to become unstable, and a stable period orbit (SPO) is born in a Hopf bifurcation off $\boldsymbol{u}_{UB}$ . This behaviour is apparent in the results of DMD at $Re=150$ in figure 12. In addition to the unstable pair of eigenvalues associated with the dynamics in the linear subspace, $\unicode[STIX]{x1D701}_{1}^{\pm }\approx 0.007\pm 0.095\text{i}$ , higher-order Koopman eigenvalues are again observed, and correspond to products of the Koopman eigenfunctions $\unicode[STIX]{x1D711}_{\unicode[STIX]{x1D701}_{1}}^{\pm }(\boldsymbol{u})$ .

The destabilisation of the upper-branch solution and the emergence of a SPO has further consequences for DMD. There are now three cross-over points associated with Koopman decompositions around each of the three invariant sets ( $\boldsymbol{u}_{LB}$ , $\boldsymbol{u}_{UB}$ and the SPO), and DMD will only ‘work’ if it is restricted to a particular expansion zone. The results highlight the care that must be taken in more complex flows with many exact coherent states buried in the turbulent attractor.

Figure 13.

(a,b) Imaginary and (c,d) real components of DMD eigenvalues obtained on shorter time windows passed through a trajectory running from the linear subspace of $\boldsymbol{u}_{LB}$ to the SPO at $Re=150$ (the variable $t_{F}$ is the ‘final’ time of each time window). (a,c) Time windows of length $T_{w}=50$ are used with $M=50$ snapshots. (b,d)  $T_{w}=200$ with $M=150$ snapshots. Snapshot spacing is $\unicode[STIX]{x1D6FF}t=1$ for both sets of calculations. The colouring serves as a guide for the eye, with red, orange and blue identifying expansions around $\boldsymbol{u}_{LB}$ , $\boldsymbol{u}_{UB}$ and the SPO respectively.

Figure 14.

DMD eigenvalues obtained from the $\boldsymbol{u}_{LB}\rightarrow$ SPO heteroclinic connection at $Re=150$ . (a) Time window $t\in [450,500]$ , with $M=50$ snapshots, $\unicode[STIX]{x1D6FF}t=1$ (cf. upper-branch spectrum in figure 12). (b) Time window $t\in [650,1000]$ , with $M=200$ snapshots, $\unicode[STIX]{x1D6FF}t=1$ . The modes highlighted in green are the Floquet multipliers $\unicode[STIX]{x1D707}\approx -0.017\pm 0.026\text{i}$ (from the form $\text{e}^{\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FF}t}$ ). Dashed lines identify harmonics of the fundamental frequency of the periodic orbit.

To demonstrate the restrictions placed on DMD, we consider again a trajectory beginning in the linear subspace around $\boldsymbol{u}_{LB}$ , which now collapses into the SPO as $t\rightarrow \infty$ . Similar to our approach for the $\boldsymbol{u}_{LB}\rightarrow \boldsymbol{u}_{UB}$ connection at $Re=135$ , we perform many DMDs in a fixed time window which is passed along the finite time approximation to the heteroclinic connection. The results of these calculations for two DMD time windows, $T_{w}\in \{50,200\}$ , are reported in figure 13.

For the shorter DMD time window, $T_{w}=50$ , three distinct trends in the eigenvalues are observed. As the trajectory moves away from the lower branch solution, the DMD locates the positive, real eigenvalue associated with the growth rate along the single unstable direction, before identifying the integer multiples of this growth rate corresponding to the higher-order Koopman modes. This behaviour is analogous to that found at $Re=135$ (see figure 5); like that earlier $\boldsymbol{u}_{LB}\rightarrow \boldsymbol{u}_{UB}$ connection, there is also a breakdown in the DMD/Koopman eigenvalues, here at $t_{F}\sim 400$ . The region where there is inconsistency between successive DMD calculations is roughly twice the length of the DMD time window, which suggests the presence of a cross-over point. Beyond $t_{F}\sim 450$ , there is a clear repeated frequency in the DMD eigenvalues, $\unicode[STIX]{x1D706}_{i}\approx 0.064$ , consistent over approximately 100 advective time units. This frequency is close to that associated with the stable spiral into the upper branch ( $\unicode[STIX]{x1D706}_{i}=0.069$ , see figure 12). However, the DMD is unable to resolve the associated decay rate correctly, or obtain the complex conjugate pair of unstable modes associated with $\boldsymbol{u}_{UB}$ at this Reynolds number. An individual eigenvalue spectrum from this region is reported in figure 14 and should be contrasted with those obtained on trajectories starting in the linear subspace of $\boldsymbol{u}_{UB}$ (figure 12). A pair of unstable eigenvalues are found, though their growth rate and frequency do not correspond to the unstable directions identified in figure 12. It is likely that the trajectory simply does not go close enough to the upper-branch equilibrium to accurately distinguish the correct form of the neutral mode ( $\boldsymbol{u}_{UB}$ itself) from the slowly growing eigenvalues, and this error contaminates the rest of the spectrum. Finally, at around $t_{F}\sim 600$ the output of the DMD calculations jumps again. There are clear repeated harmonics of a fundamental frequency, $\unicode[STIX]{x1D714}_{f}=0.085$ ( $\unicode[STIX]{x1D706}_{r}$ is very close to zero), which corresponds to a periodic orbit with period $T=73.9$ . Occasionally, the DMD erroneously identifies growing modes, while the lattice of decaying eigenvalues one would expect to find around a stable limit cycle is absent (e.g. as observed by Bagheri (Reference Bagheri2013) and proven for systems with a single attracting limit cycle by Mezic (Reference Mezic2017)).

To accurately determine the Koopman eigenvalues around the SPO, a much longer time window is required. For example, the longer time window considered in figure 13, $T_{w}=200$ , no longer shows exponentially unstable modes in the collapse onto the SPO. Instead, there is an array of decaying eigenvalues, all with decay rate $\unicode[STIX]{x1D706}_{r}\approx -0.017$ . Each harmonic of the SPO is flanked by a pair of decaying modes, $\unicode[STIX]{x1D706}=n\unicode[STIX]{x1D714}_{f}\text{i}+(\unicode[STIX]{x1D707}_{r}\pm \unicode[STIX]{x1D707}_{i}\text{i})$ , indicating the presence of a pair of stable Floquet multipliers $\text{e}^{\unicode[STIX]{x1D707}T}$ , with $\unicode[STIX]{x1D707}=-0.017\pm 0.026\text{i}$ . To obtain higher-order Koopman modes associated with the SPO (see Bagheri Reference Bagheri2013), even longer time windows are required. For example, the decay rate $2\unicode[STIX]{x1D707}_{r}$ is observed in figure 14 with a slightly longer time window $T_{w}=350$ , although only approximately and there are eigenvalues missing.

A consequence of the longer time window required for a more accurate resolution of the Koopman eigenvalues around the SPO is the loss of any indication of the presence of the upper branch in figure 13. The time window is longer than the residence time in the upper-branch expansion region, so all DMD calculations that see $\boldsymbol{u}_{UB}$ contain at least one cross-over point. In addition to a large time window, $T_{w}>T$ , DMD calculations which are able to accurately resolve the Koopman eigenvalues around the periodic orbit require many snapshot pairs. In a turbulent flow with unstable periodic orbits, each with many different Floquet multipliers, these requirements are unlikely to be achievable in practice. However, the fact that DMD time windows which are shorted than the period, $T_{w}<T$ , are still able to identify the fundamental frequencies and associated mode shapes indicates that DMD may be a useful alternative to a recurrent flow analysis in generating guesses for unstable periodic orbits.