2.1 Space-only POD

The most commonly employed form of POD generates spatial modes $\unicode[STIX]{x1D753}(\boldsymbol{x})$ . This is accomplished by defining the stochastic ensemble to consist of snapshots of the flow field at different time instances. In other words, the flow at each instant is treated as a realization of a stochastic process. The appropriate inner product is then

(2.4) $$\begin{eqnarray}\langle \boldsymbol{u},\boldsymbol{v}\rangle _{x}=\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{v}^{\ast }(\boldsymbol{x},t)\unicode[STIX]{x1D652}(\boldsymbol{x})\boldsymbol{u}(\boldsymbol{x},t)\,\text{d}\boldsymbol{x},\end{eqnarray}$$

where $\boldsymbol{u}$ and $\boldsymbol{v}$ are any two elements in ${\mathcal{H}}$ , $\unicode[STIX]{x1D6FA}$ denotes the spatial domain over which the flow is defined, and the weight $\unicode[STIX]{x1D652}$ is a positive-definite Hermitian tensor of appropriate dimension. We will restrict our attention to bounded spatial domains, but note that unbounded homogeneous dimensions can be accommodated by transforming those directions to Fourier space (Lumley Reference Lumley, Yaglom and Tatarski1967, Reference Lumley1970; George Reference George2017). The expectation operator for this definition of the stochastic ensemble is simply a time average, so we are restricted to statistically stationary flows.

The quantity to maximize is

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D706}=\frac{E\{|\langle \boldsymbol{q}(\boldsymbol{x},t),\unicode[STIX]{x1D753}(\boldsymbol{x})\rangle _{x}|^{2}\}}{\langle \unicode[STIX]{x1D753}(\boldsymbol{x}),\unicode[STIX]{x1D753}(\boldsymbol{x})\rangle _{x}}\end{eqnarray}$$

and the resulting Fredholm eigenvalue problem is

(2.6) $$\begin{eqnarray}\int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x}^{\prime })\unicode[STIX]{x1D652}(\boldsymbol{x}^{\prime })\unicode[STIX]{x1D753}(\boldsymbol{x}^{\prime })\,\text{d}\boldsymbol{x}^{\prime }=\unicode[STIX]{x1D706}\unicode[STIX]{x1D753}(\boldsymbol{x}),\end{eqnarray}$$

where

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x}^{\prime })=E\{\boldsymbol{q}(\boldsymbol{x},t)\boldsymbol{q}^{\ast }(\boldsymbol{x}^{\prime },t)\}\end{eqnarray}$$

is the two-point spatial correlation tensor. This tensor is a nuclear kernel, i.e. it is compact and $\int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x})\,\text{d}\boldsymbol{x}<\infty$ . As a result, Hilbert–Schmidt theory guarantees that the eigenmodes satisfying (2.6) have a number of special properties. First, there exists a countably infinite set of eigenmodes, $\{\unicode[STIX]{x1D753}_{j},\unicode[STIX]{x1D706}_{j}\}$ , that can be ranked according to their eigenvalue, $\unicode[STIX]{x1D706}_{1}\geqslant \unicode[STIX]{x1D706}_{2}\geqslant \cdots \geqslant 0$ . Each eigenvalue gives the average energy captured by that mode, measured in the spatial norm induced by the inner product (2.4), and the total energy of the flow is given by the sum of the eigenvalues.

The eigenvectors are orthogonal, $\langle \unicode[STIX]{x1D753}_{j},\unicode[STIX]{x1D753}_{k}\rangle _{x}=\unicode[STIX]{x1D6FF}_{jk}$ , and provide a complete basis for $\boldsymbol{q}$ . Accordingly, the flow field can be expanded as

(2.8) $$\begin{eqnarray}\boldsymbol{q}(\boldsymbol{x},t)=\mathop{\sum }_{j=1}^{\infty }a_{j}(t)\unicode[STIX]{x1D753}_{j}(\boldsymbol{x}),\end{eqnarray}$$

with $a_{j}(t)=\langle \boldsymbol{q}(\boldsymbol{x},t),\unicode[STIX]{x1D753}_{j}(\boldsymbol{x})\rangle _{x}$ . This expansion is optimal in its ability to capture the flow energy; if the expansion is truncated at order $n$ , any other orthogonal expansion of the same order will capture less energy. The expansion coefficients are uncorrelated at zero time lag,

(2.9) $$\begin{eqnarray}E\{a_{j}(t)a_{k}^{\ast }(t)\}=\unicode[STIX]{x1D706}_{j}\unicode[STIX]{x1D6FF}_{jk}.\end{eqnarray}$$

Finally, the eigenmodes provide a diagonal representation of the two-point spatial correlation tensor

(2.10) $$\begin{eqnarray}\unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x}^{\prime })=\mathop{\sum }_{j=1}^{\infty }\unicode[STIX]{x1D706}_{j}\unicode[STIX]{x1D753}_{j}(\boldsymbol{x})\unicode[STIX]{x1D753}_{j}^{\ast }(\boldsymbol{x}^{\prime })\end{eqnarray}$$

and are therefore its principal components. Accordingly, space-only POD modes optimally represent spatial correlations within the flow.

2.2 Spectral proper orthogonal decomposition

Alternatively, we can seek modes that depend on both space and time. This is accomplished by defining the stochastic ensemble to consist of a collection of realizations of the time-dependent flow. For example, different runs of the same experiment are considered to be realizations of a stochastic process. The appropriate inner product is then

(2.11) $$\begin{eqnarray}\langle \boldsymbol{u},\boldsymbol{v}\rangle _{x,t}=\int _{-\infty }^{\infty }\int _{\unicode[STIX]{x1D6FA}}\boldsymbol{v}^{\ast }(\boldsymbol{x},t)\unicode[STIX]{x1D652}(\boldsymbol{x})\boldsymbol{u}(\boldsymbol{x},t)\,\text{d}\boldsymbol{x}\,\text{d}t\end{eqnarray}$$

and the expectation operator is an ensemble average over different stochastic realizations of the flow. The quantity to maximize is then

(2.12) $$\begin{eqnarray}\unicode[STIX]{x1D706}=\frac{E\{|\langle \boldsymbol{q}(\boldsymbol{x},t),\unicode[STIX]{x1D753}(\boldsymbol{x},t)\rangle _{x,t}|^{2}\}}{\langle \unicode[STIX]{x1D753}(\boldsymbol{x},t),\unicode[STIX]{x1D753}(\boldsymbol{x},t)\rangle _{x,t}},\end{eqnarray}$$

which leads to the eigenvalue problem

(2.13) $$\begin{eqnarray}\int _{-\infty }^{\infty }\int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x}^{\prime },t,t^{\prime })\unicode[STIX]{x1D652}(\boldsymbol{x}^{\prime })\unicode[STIX]{x1D753}(\boldsymbol{x}^{\prime },t^{\prime })\,\text{d}\boldsymbol{x}^{\prime }\,\text{d}t^{\prime }=\unicode[STIX]{x1D706}\unicode[STIX]{x1D753}(\boldsymbol{x},t),\end{eqnarray}$$

where

(2.14) $$\begin{eqnarray}\unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x}^{\prime },t,t^{\prime })=E\{\boldsymbol{q}(\boldsymbol{x},t)\boldsymbol{q}^{\ast }(\boldsymbol{x}^{\prime },t^{\prime })\}\end{eqnarray}$$

is the two-point space–time correlation tensor.

Because they persist indefinitely, statistically stationary flows have infinite energy in a space–time norm. As a result, the space–time correlation tensor (2.14) is not nuclear and the eigenmodes of (2.13) do not posses the properties generally associated with POD. To remedy this, a new eigenvalue problem can be obtained in spectral space from which modes with useful properties can be obtained. In what follows, we focus exclusively on stationary flows.

To derive the spectral eigenvalue problem, we recall that the correlation tensor of a wide-sense stationary flow depends only on the difference between two times,

(2.15) $$\begin{eqnarray}\unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x}^{\prime },t,t^{\prime })\rightarrow \unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x}^{\prime },t-t^{\prime }).\end{eqnarray}$$

Then, the cross-spectral density tensor $\unicode[STIX]{x1D64E}$ can be defined as the Fourier transform pair of the correlation tensor,

(2.16) $$\begin{eqnarray}\unicode[STIX]{x1D64E}(\boldsymbol{x},\boldsymbol{x}^{\prime },f)=\int _{-\infty }^{\infty }\unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x}^{\prime },\unicode[STIX]{x1D70F})\text{e}^{-\text{i}2\unicode[STIX]{x03C0}f\unicode[STIX]{x1D70F}}\,\text{d}\unicode[STIX]{x1D70F}\end{eqnarray}$$

and

(2.17) $$\begin{eqnarray}\unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x}^{\prime },\unicode[STIX]{x1D70F})=\int _{-\infty }^{\infty }\unicode[STIX]{x1D64E}(\boldsymbol{x},\boldsymbol{x}^{\prime },f)\text{e}^{\text{i}2\unicode[STIX]{x03C0}f\unicode[STIX]{x1D70F}}\,\text{d}f.\end{eqnarray}$$

Using these definitions, the following result can be derived: for any frequency $f^{\prime }$ , the function $\unicode[STIX]{x1D753}(\boldsymbol{x},t)=\unicode[STIX]{x1D74D}(\boldsymbol{x},f^{\prime })\text{e}^{\text{i}2\unicode[STIX]{x03C0}f^{\prime }t}$ is a solution of the eigenvalue problem (2.13) with eigenvalue $\unicode[STIX]{x1D706}(f^{\prime })$ , where $\unicode[STIX]{x1D74D}(\boldsymbol{x},f^{\prime })$ and $\unicode[STIX]{x1D706}(f^{\prime })$ satisfy the spectral eigenvalue problem

(2.18) $$\begin{eqnarray}\int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D64E}(\boldsymbol{x},\boldsymbol{x}^{\prime },f^{\prime })\unicode[STIX]{x1D652}(\boldsymbol{x}^{\prime })\unicode[STIX]{x1D74D}(\boldsymbol{x}^{\prime },f^{\prime })\,\text{d}\boldsymbol{x}^{\prime }=\unicode[STIX]{x1D706}(f^{\prime })\unicode[STIX]{x1D74D}(\boldsymbol{x},f^{\prime }).\end{eqnarray}$$

This result was first given by Lumley (Reference Lumley, Yaglom and Tatarski1967, Reference Lumley1970); we offer an alternative derivation in appendix A.

The cross-spectral density tensor is nuclear, so the eigenmodes of (2.18) at each frequency inherit properties analogous to those of space-only POD modes. There is a countably infinite set of eigenfunctions $\unicode[STIX]{x1D74D}_{j}(\boldsymbol{x},f)$ at each frequency that are orthogonal to all other modes at the same frequency in the spatial inner product (2.4), i.e. $\langle \unicode[STIX]{x1D74D}_{j},\unicode[STIX]{x1D74D}_{k}\rangle _{x}=\unicode[STIX]{x1D6FF}_{jk}$ . The Fourier modes of each flow realization are optimally expanded as

(2.19) $$\begin{eqnarray}\hat{\boldsymbol{q}}(\boldsymbol{x},f)=\mathop{\sum }_{j=1}^{\infty }a_{j}(f)\unicode[STIX]{x1D74D}_{j}(\boldsymbol{x},f),\end{eqnarray}$$

with $a_{j}(f)=\langle \hat{\boldsymbol{q}}(\boldsymbol{x},f),\unicode[STIX]{x1D74D}_{j}(\boldsymbol{x},f)\rangle _{x}$ , and the expansion coefficients are uncorrelated,

(2.20) $$\begin{eqnarray}E\{a_{j}(f)a_{k}^{\ast }(f)\}=\unicode[STIX]{x1D706}_{j}(f)\unicode[STIX]{x1D6FF}_{jk}.\end{eqnarray}$$

The cross-spectral density tensor has the diagonal representation

(2.21) $$\begin{eqnarray}\unicode[STIX]{x1D64E}(\boldsymbol{x},\boldsymbol{x}^{\prime },f)=\mathop{\sum }_{j=1}^{\infty }\unicode[STIX]{x1D706}_{j}(f)\unicode[STIX]{x1D74D}_{j}(\boldsymbol{x},f)\unicode[STIX]{x1D74D}_{j}^{\ast }(\boldsymbol{x}^{\prime },f),\end{eqnarray}$$

so the SPOD modes are its principal components. Furthermore, the modes $\unicode[STIX]{x1D74D}_{j}(\boldsymbol{x},f)\text{e}^{\text{i}2\unicode[STIX]{x03C0}ft}$ are orthogonal in the space–time inner product (2.11), so each mode at each frequency can be viewed as a distinct space–time mode. The space–time correlation tensor can be written as

(2.22) $$\begin{eqnarray}\unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x}^{\prime },\unicode[STIX]{x1D70F})=\int _{-\infty }^{\infty }\mathop{\sum }_{j=1}^{\infty }\unicode[STIX]{x1D706}_{j}(f)\unicode[STIX]{x1D74D}_{j}(\boldsymbol{x},f)\unicode[STIX]{x1D74D}_{j}^{\ast }(\boldsymbol{x}^{\prime },f)\text{e}^{\text{i}2\unicode[STIX]{x03C0}f\unicode[STIX]{x1D70F}}\,\text{d}f.\end{eqnarray}$$

In summary, for stationary flows, the space–time POD formulation leads to spectral POD modes that each oscillate at a single frequency and optimally represent the second-order space–time flow statistics.

2.3 Spectral versus space-only POD

The essential difference between spectral and space-only POD is that the former yields modes that are coherent in space and time, whereas the later gives modes that are only spatially coherent. First, consider space-only POD. The time dependence of the flow field $\boldsymbol{q}(\boldsymbol{x},t)$ is treated as a stochastic parameter, with different time instances taken to represent different members in an ensemble of spatially dependent fields. Once interpreted in this way, the flow snapshots that make up the ensemble lose any concept of sequential ordering, so the time-dependent evolution of the flow has no impact on the definition of the POD modes. Therefore, the POD modes are impervious to temporal correlation within the data, which is an essential feature of physical coherent structures. Because of this, POD modes do not necessarily represent structures that evolve coherently in space and time. This can be shown explicitly by writing the space–time correlation tensor in terms of space-only POD modes,

(2.23a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x}^{\prime },\unicode[STIX]{x1D70F}) & = & \displaystyle E\left\{\left(\mathop{\sum }_{j=1}^{\infty }a_{j}(t)\unicode[STIX]{x1D753}_{j}(\boldsymbol{x})\right)\left(\mathop{\sum }_{k=1}^{\infty }a_{k}(t+\unicode[STIX]{x1D70F})\unicode[STIX]{x1D753}_{k}(\boldsymbol{x}^{\prime })\right)^{\ast }\right\}\end{eqnarray}$$

(2.23b ) $$\begin{eqnarray}\displaystyle & = & \displaystyle \mathop{\sum }_{j=1}^{\infty }\mathop{\sum }_{k=1}^{\infty }\unicode[STIX]{x1D60A}_{a_{j}a_{k}}^{POD}(\unicode[STIX]{x1D70F})\unicode[STIX]{x1D753}_{j}(\boldsymbol{x})\unicode[STIX]{x1D753}_{k}^{\ast }(\boldsymbol{x}^{\prime }),\end{eqnarray}$$

(2.24) $$\begin{eqnarray}\unicode[STIX]{x1D60A}_{a_{j}a_{k}}^{POD}(\unicode[STIX]{x1D70F})=E\{a_{j}(t)a_{k}^{\ast }(t+\unicode[STIX]{x1D70F})\}.\end{eqnarray}$$

When $\unicode[STIX]{x1D70F}=0$ , (2.9) ensures that $\unicode[STIX]{x1D60A}_{a_{j}a_{k}}^{POD}(\unicode[STIX]{x1D70F})=\unicode[STIX]{x1D706}_{j}\unicode[STIX]{x1D6FF}_{jk}$ and (2.23) reduces to the spatial correlation tensor. On the other hand, (2.9) is not applicable when $\unicode[STIX]{x1D70F}\neq 0$ , and, as a result, POD theory does not guarantee any special properties for $\unicode[STIX]{x1D60A}_{a_{j}a_{k}}^{POD}(\unicode[STIX]{x1D70F})$ ; thus, the temporal correlation of two terms in the POD expansion is not known a priori. This means that the part of the flow described by a given POD mode is not necessarily correlated with the part of the flow described by the same POD mode at a later time, nor is it necessarily uncorrelated with the part of the flow described by a different mode at a later time. Therefore, contrary to some previous statements, space-only POD modes do not necessarily represent flow structures that evolve coherently. A recent analysis of space-only POD by George (Reference George2017) erroneously assumed the expansion coefficients to be uncorrelated at different times.

In contrast, SPOD modes do represent structures that evolve coherently. The space–time correlation tensor is written in terms of SPOD modes in (2.22). This form of the space–time correlation tensor is the result of special properties of the correlations

(2.25) $$\begin{eqnarray}\unicode[STIX]{x1D60A}_{a_{j}a_{k}}^{SPOD}(f,f^{\prime })\triangleq E\{a_{j}(f)a_{k}^{\ast }(f^{\prime })\}=\unicode[STIX]{x1D706}_{j}\unicode[STIX]{x1D6FF}_{jk}\unicode[STIX]{x1D6FF}(f-f^{\prime }),\end{eqnarray}$$

which govern the correlation between the part of the flow described by individual SPOD modes. Specifically, (2.22) can be obtained by inserting (2.25) along with the inverse Fourier transform of (2.19) into the definition of the space–time correlation tensor given by (2.14) and (2.15). The first two terms in the final form of (2.25) follow from (2.20) and ensure that two terms in the SPOD expansion at the same frequency are uncorrelated. The final term is a consequence of the fact that the frequency components from the Fourier transform of a stationary random process are uncorrelated (Lumley Reference Lumley1970; George Reference George1988).

In sum, (2.22) and (2.25) show that the part of the flow described by a particular SPOD mode is perfectly correlated with the part of the flow described by that same mode at all times and entirely uncorrelated with the part of the flow described by all other modes at all times. In other words, each SPOD mode describes a structure that evolves coherently in space and time.

The preceding analysis does not imply that space-only POD modes can never exhibit space–time coherence. For example, Rowley et al. (Reference Rowley, Colonius and Murray2004) observed that some of the leading space-only POD modes of a compressible cavity flow capture single-frequency Rossiter modes. Rather, our analysis shows that space-only POD modes do not have this property by construction, in contrast to SPOD modes which evolve coherently in space and time by construction.

It is also possible to derive equations relating space-only and spectral POD modes and eigenvalues. Using the fact that the spatial correlation tensor is equivalent to the zero-time-lag space–time correlation tensor, (2.17) implies that

(2.26) $$\begin{eqnarray}\unicode[STIX]{x1D63E}(\boldsymbol{x},\boldsymbol{x}^{\prime })=\int _{-\infty }^{\infty }\unicode[STIX]{x1D64E}(\boldsymbol{x},\boldsymbol{x}^{\prime },f)\,\text{d}f.\end{eqnarray}$$

Expanding $\unicode[STIX]{x1D63E}$ and $\unicode[STIX]{x1D64E}$ in terms of POD and SPOD modes respectively gives

(2.27) $$\begin{eqnarray}\mathop{\sum }_{j=1}^{\infty }\unicode[STIX]{x1D706}_{j}\unicode[STIX]{x1D753}_{j}(\boldsymbol{x})\unicode[STIX]{x1D753}_{j}(\boldsymbol{x}^{\prime })^{\ast }=\int _{-\infty }^{\infty }\mathop{\sum }_{k=1}^{\infty }\unicode[STIX]{x1D706}_{k}(f)\unicode[STIX]{x1D74D}_{k}(\boldsymbol{x},f)\unicode[STIX]{x1D74D}_{k}(\boldsymbol{x}^{\prime },f)^{\ast }\,\text{d}f.\end{eqnarray}$$

Applying the operation $\langle \unicode[STIX]{x1D753}_{j},\cdot \rangle _{x}$ to both sides of (2.27) and dividing by $\unicode[STIX]{x1D706}_{j}$ leads to the expression

(2.28) $$\begin{eqnarray}\unicode[STIX]{x1D753}_{j}(\boldsymbol{x})=\int _{-\infty }^{\infty }\mathop{\sum }_{k=1}^{\infty }\frac{\unicode[STIX]{x1D706}_{k}(f)}{\unicode[STIX]{x1D706}_{j}}c_{jk}(f)\unicode[STIX]{x1D74D}_{k}(\boldsymbol{x},f)\,\text{d}f,\end{eqnarray}$$

where $c_{jk}(f)=\langle \unicode[STIX]{x1D753}_{j}(\boldsymbol{x}),\unicode[STIX]{x1D74D}_{k}(\boldsymbol{x},f)\rangle _{x}$ . Taking the same inner product again and moving $\unicode[STIX]{x1D706}_{j}$ back to the left-hand side yields

(2.29) $$\begin{eqnarray}\unicode[STIX]{x1D706}_{j}=\int _{-\infty }^{\infty }\mathop{\sum }_{k=1}^{\infty }\unicode[STIX]{x1D706}_{k}(f)|c_{jk}(f)|^{2}\,\text{d}f.\end{eqnarray}$$

Equations (2.28) and (2.29) show that each space-only POD mode is potentially made up of many SPOD modes. Physically, this means that the spatially coherent structures represented by space-only POD are composed of contributions from spatiotemporal coherent structures at many frequencies. Practically, this is manifested as broadband frequency content within the coefficients $a_{j}(t)$ . This highlights the fact that each space-only POD mode typically represents flow phenomena at many different time scales, which muddies their interpretation. In contrast, SPOD modes decouple flow phenomena at different time scales, which can be helpful for understanding the flow dynamics and deriving non-empirical models.

4.1 DMD and the DFT

For zero-mean data that are uniformly sampled in time, Chen, Tu & Rowley (Reference Chen, Tu and Rowley2012) showed that DMD is formally equivalent to the DFT (the flow snapshots must also be linearly independent for this to hold, which is true in most applications). As a result, each DMD mode has zero growth/decay rate. This formal DMD–DFT equivalence does not hold for data with a non-zero mean, but in practice the physically relevant DMD modes tend to have nearly zero growth/decay rates for stationary flows (e.g. Rowley et al. Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009; Chen et al. Reference Chen, Tu and Rowley2012; Schmid, Violato & Scarano Reference Schmid, Violato and Scarano2012; Semeraro, Bellani & Lundell Reference Semeraro, Bellani and Lundell2012), which is logical since stationary flows are persistent by definition.

This tendency is explained more rigorously by the connection between DMD and Koopman operator theory (Mezić Reference Mezić2005; Rowley et al. Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009). The Koopman operator is an infinite-dimensional linear operator that describes the evolution of scalar observables of a nonlinear dynamical system on a finite manifold. The eigenvectors of the Koopman operator can be used to decompose vector-valued observables (equivalent to our $\boldsymbol{q}$ ) into modes that evolve with a single frequency and growth/decay rate. Mezić (Reference Mezić2005) showed that for any dynamical system with a Borel probability measure, the growth/decay rate is zero and Koopman modes are equivalent to Fourier modes. Stationary flows possess an ergodic measure by definition, so their Koopman modes are simply Fourier modes.

Rowley et al. (Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009) showed that DMD modes approximate Koopman modes when the flow snapshots used to compute the modes are linearly independent, and Tu et al. (Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2014) showed that this holds under a slightly weaker condition on the data termed linear consistency. In light of this, it is not surprising that DMD modes tend to be similar to DFT modes for stationary flows. Moreover, deviations of the DMD modes from Fourier modes can be viewed as artefacts of the DMD approximation of the underlying Koopman operator. This suggests that, contrary to prevailing wisdom, it is advantageous to subtract the mean when applying DMD to stationary flows to ensure that the DMD modes mimic the zero-growth-rate property of the underlying Koopman modes.

4.2 An ensemble DMD problem for stationary flows

To conceptually relate SPOD and DMD, imagine that we have an ensemble of $N_{e}$ data sets, each representing a stochastic realization of the same stationary flow. There are at least two approaches one could adopt for applying DMD to this problem. The first would be to simply perform separate DMD computations for each realization of the flow. If we agree to subtract the mean to ensure zero growth rates as described previously, each DMD mode will exactly correspond to a DFT mode, but in general the mode at a given frequency will be different for each realization.

A more sophisticated approach for applying DMD to this problem would be to use the approach of Tu et al. (Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2014) to combine multiple flow realization into a single DMD calculation. Their variant of DMD, which they call ‘exact DMD’, is defined in terms of the operator

(4.1) $$\begin{eqnarray}\mathbf{A}\triangleq \mathbf{Y}\mathbf{X}^{+},\end{eqnarray}$$

where $\mathbf{X}^{+}$ is the pseudo-inverse of $\mathbf{X}$ and the columns of the matrices $\mathbf{X}$ and $\mathbf{Y}$ are input–output data pairs that are related by a linear operator that is to be approximated by the matrix $\mathbf{A}$ . The DMD modes and eigenvalues are then given by the eigenvectors and eigenvalues of $\mathbf{A}$ respectively.

For a standard application in which the flow data consist of sequential snapshots of the flow, the input and output matrices are

(4.2a ) $$\begin{eqnarray}\displaystyle \mathbf{X} & \triangleq & \displaystyle [\mathbf{q}_{1},\mathbf{q}_{2},\ldots ,\mathbf{q}_{M-1}]=\mathbf{Q}\mathbf{T}_{\text{X}},\end{eqnarray}$$

(4.2b ) $$\begin{eqnarray}\displaystyle \mathbf{Y} & \triangleq & \displaystyle [\mathbf{q}_{2},\mathbf{q}_{3},\ldots ,\mathbf{q}_{M}]=\mathbf{Q}\mathbf{T}_{\text{Y}},\end{eqnarray}$$

$\mathbf{q}_{j}\in \mathbb{R}^{N}$

$\mathbf{Q}\in \mathbb{R}^{N\times M}$

$\mathbf{T}_{\text{X}},\mathbf{T}_{\text{Y}}\in \mathbb{R}^{M\times M-1}$

$\mathbf{Q}$

(4.3a,b ) $$\begin{eqnarray}\mathbf{T}_{\text{X}}\triangleq \left[\begin{array}{@{}cccc@{}}1 & 0 & \cdots \, & 0\\ 0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0\\ 0 & \cdots \, & 0 & 1\\ 0 & 0 & \cdots \, & 0\end{array}\right],\quad \mathbf{T}_{\text{Y}}\triangleq \left[\begin{array}{@{}cccc@{}}0 & 0 & \cdots \, & 0\\ 1 & 0 & \cdots \, & 0\\ 0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0\\ 0 & \cdots \, & 0 & 1\end{array}\right].\end{eqnarray}$$

As described by Tu et al. (Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2014), multiple realizations of a flow can be accommodated within the exact DMD framework by arranging the realizations together into ensemble input and output matrices,

(4.4a ) $$\begin{eqnarray}\displaystyle \mathbf{X} & \triangleq & \displaystyle [\mathbf{Q}^{(1)}\mathbf{T}_{\text{X}},\ldots ,\mathbf{Q}^{(N_{e})}\mathbf{T}_{\text{X}}],\end{eqnarray}$$

(4.4b ) $$\begin{eqnarray}\displaystyle \mathbf{Y} & \triangleq & \displaystyle [\mathbf{Q}^{(1)}\mathbf{T}_{\text{Y}},\ldots ,\mathbf{Q}^{(N_{e})}\mathbf{T}_{\text{Y}}],\end{eqnarray}$$

$\mathbf{Q}^{(n)}\in \mathbb{R}^{N\times M}$

$n$

The DMD/DFT equivalence for zero-mean data proven by Chen et al. (Reference Chen, Tu and Rowley2012) was derived within the context of the original Arnoldi-based formulation of DMD given by Schmid (Reference Schmid2010) and was restricted to the case of sequential linearly independent snapshots. Accordingly, it does not apply to the ensemble DMD problem defined by (4.1) and (4.4).

In appendix B, we prove that the DMD modes obtained from this ensemble formulation are precisely the DFT modes of each realization of the flow if each realization has zero mean (i.e. we have subtracted the mean in line with our earlier arguments for stationary flows) and the ensemble input/output data matrices are linearly consistent. Under these conditions, (4.1) and (4.4) reduce to

(4.5) $$\begin{eqnarray}\mathbf{A}\left[\tilde{\mathbf{Q}}^{(1)},\ldots ,\tilde{\mathbf{Q}}^{(N_{e})}\right]=\left[\tilde{\mathbf{Q}}^{(1)},\ldots ,\tilde{\mathbf{Q}}^{(N_{e})}\right]\left[\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D6B2}_{\text{DMD}} & & \\ & \ddots & \\ & & \unicode[STIX]{x1D6B2}_{\text{DMD}}\end{array}\right],\end{eqnarray}$$

where the columns of $\tilde{\mathbf{Q}}^{(n)}$ contain the DFT modes of $\mathbf{Q}^{(n)}$ excluding the mean component,

(4.6) $$\begin{eqnarray}\unicode[STIX]{x1D6B2}_{\text{DMD}}=\left[\begin{array}{@{}cccc@{}}z & & & \\ & z^{2} & & \\ & & \ddots & \\ & & & z^{M-1}\end{array}\right],\end{eqnarray}$$

and $z=\text{e}^{-\text{i}2\unicode[STIX]{x03C0}/M}$ is the primitive $M$ th root of unity. Therefore, each DFT mode $\hat{\mathbf{q}}_{k}^{(n)}$ , $k=2,\ldots ,M$ , is an eigenvector of $\mathbf{A}$ with eigenvalue $z^{k-1}$ . These eigenvalues have unit magnitude, so the growth/decay rate is zero, as expected. The frequency corresponding to each DFT mode is recovered from the eigenvalue as $f_{k}=\text{Re}[\text{ln}(z^{k-1})/(-\text{i}2\unicode[STIX]{x03C0}\unicode[STIX]{x0394}t)]$ , and it is easy to verify that these frequencies match those given by (3.5). Therefore, each DFT mode from each realization of the flow is a DMD mode that oscillates at the corresponding DFT frequency. Since each eigenvalue is repeated $N_{e}$ times, any linear combination of the DFT modes at a given frequency is also a DMD mode.

4.3 DMD and SPOD

Overall, for either the simple approach in which a separate $\mathbf{A}$ operator is defined for each flow realization or the more sophisticated approach in which multiple realizations are incorporated into a single ensemble DMD computation, the end result is a set of DMD modes at each frequency consisting of the DFT mode of each realization (or a linear combination thereof). Because of the random nature of turbulent flows and the uncertainty inherent to the DFT, the mode for a given frequency obtained from each flow realization will in general be different. This statistical variability may seem undesirable, but in fact it exposes an important characteristic of turbulent flows – the behaviour at a given frequency cannot be described by a single deterministic mode, i.e. the cross-spectral density tensor is not rank one. In light of this, how can useful information be extracted from this set of differing DMD modes? A sensible approach is to search for functions that best represent the ensemble of DMD modes at each frequency; these are given precisely by the SPOD modes at that frequency. In other words, SPOD modes provide the optimal basis, as defined by (2.12), for describing the variability within the ensemble of DMD modes. Said the other way around, each DMD mode in the ensemble at a particular frequency represents one of the possible ways in which the coherent structures represented by the SPOD modes can coexist in one particular realization of the turbulent flow.

The preceding analysis of the ensemble DMD problem enables us to make a still stronger connection between DMD and SPOD for stationary flows. As already mentioned, any linear combination of the set of DFT modes at a given frequency is also a DMD mode. Recall from (3.11) that each SPOD mode can be written as a linear combination of the DFT modes used to estimate the cross-spectral density tensor. Therefore, each SPOD mode is itself a DMD mode. Specifically, SPOD modes are given by the particular linear combination, or weighted average, of DFT/DMD modes that yields uncorrelated modes that optimally capture the flow energy, as defined by (2.12). Therefore, SPOD modes are ensemble DMD modes that have been optimally averaged to provide the best possible representation of the second-order space–time flow statistics. This means that SPOD modes are dynamically significant in the same sense as standard DMD modes, but at the same time optimally describe the random nature of turbulent flows. This connection between SPOD and DMD suggests that SPOD could be directly related to the concept of the stochastic Koopman operator (Mezić Reference Mezić2005), but this remains to be explored in detail.

5.1 Resolvent analysis

We begin with nonlinear flow equations of the form

(5.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D74C}}{\unicode[STIX]{x2202}t}=\boldsymbol{{\mathcal{F}}}(\unicode[STIX]{x1D74C}),\end{eqnarray}$$

where $\unicode[STIX]{x1D74C}$ is a state vector of flow variables. The compressible Navier–Stokes equations are naturally written in the form of (5.1), whereas the incompressible Navier–Stokes equations can be written in this form by eliminating the pressure by projecting the velocity field onto a divergence-free basis. Substituting the standard Reynolds decomposition

(5.2) $$\begin{eqnarray}\unicode[STIX]{x1D74C}(\boldsymbol{x},t)=\bar{\unicode[STIX]{x1D74C}}(\boldsymbol{x})+\unicode[STIX]{x1D74C}^{\prime }(\boldsymbol{x},t)\end{eqnarray}$$

into (5.1) and isolating the terms that are linear in $\unicode[STIX]{x1D74C}^{\prime }$ yields an equation of the form

(5.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D74C}^{\prime }}{\unicode[STIX]{x2202}t}-{\mathcal{A}}(\bar{\unicode[STIX]{x1D74C}})\unicode[STIX]{x1D74C}^{\prime }=\boldsymbol{h}(\bar{\unicode[STIX]{x1D74C}},\unicode[STIX]{x1D74C}^{\prime }),\end{eqnarray}$$

where

(5.4) $$\begin{eqnarray}{\mathcal{A}}(\bar{\unicode[STIX]{x1D74C}})=\frac{\unicode[STIX]{x2202}\boldsymbol{{\mathcal{F}}}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D74C}}(\bar{\unicode[STIX]{x1D74C}})\end{eqnarray}$$

is the linearized flow operator and $\boldsymbol{h}$ contains the remaining nonlinear terms as well as any exogenous inputs such as incoming fluctuations at boundaries or environmental disturbances.

To make the ensuing analysis as flexible as possible, it is useful to introduce two new variables: an input variable $\boldsymbol{\unicode[STIX]{x1D702}}$ and an output variable $\boldsymbol{y}$ , defined by the relations

(5.5) $$\begin{eqnarray}\boldsymbol{h}(\boldsymbol{x},t)={\mathcal{B}}\boldsymbol{\unicode[STIX]{x1D702}}(\boldsymbol{x},t)\end{eqnarray}$$

and

(5.6) $$\begin{eqnarray}\boldsymbol{y}(\boldsymbol{x},t)={\mathcal{C}}\unicode[STIX]{x1D74C}^{\prime }(\boldsymbol{x},t)\end{eqnarray}$$

respectively. The linear operator ${\mathcal{C}}={\mathcal{C}}(\boldsymbol{x})$ can be used to select certain flow variables and/or spatial regions that are of particular interest, or in general any linear function of the state perturbation. Similarly, the linear operator ${\mathcal{B}}={\mathcal{B}}(\boldsymbol{x})$ limits the allowable form of $\boldsymbol{h}$ , which can be useful for enforcing known properties of the nonlinearity in certain flows, e.g. that the incompressible continuity equation is linear so that the corresponding component of $\boldsymbol{h}$ must be zero. Since the flow is stationary, all quantities can be Fourier decomposed, and (5.3), (5.5) and (5.6) can be written in the frequency domain as

(5.7) $$\begin{eqnarray}(\text{i}2\unicode[STIX]{x03C0}f{\mathcal{I}}-{\mathcal{A}})\hat{\unicode[STIX]{x1D74C}}=\hat{\boldsymbol{h}}={\mathcal{B}}\hat{\boldsymbol{\unicode[STIX]{x1D702}}}\end{eqnarray}$$

and

(5.8) $$\begin{eqnarray}\hat{\boldsymbol{y}}={\mathcal{C}}\hat{\unicode[STIX]{x1D74C}}.\end{eqnarray}$$

Following McKeon & Sharma (Reference McKeon and Sharma2010), the basic objective of resolvent analysis is to find pairs of inputs and outputs at each frequency that are optimal in terms of their linear gain

(5.9) $$\begin{eqnarray}\unicode[STIX]{x1D70E}^{2}=\frac{\langle \,\hat{\boldsymbol{y}},\hat{\boldsymbol{y}}\rangle _{y}}{\langle \hat{\boldsymbol{\unicode[STIX]{x1D702}}},\hat{\boldsymbol{\unicode[STIX]{x1D702}}}\rangle _{\unicode[STIX]{x1D702}}},\end{eqnarray}$$

where

(5.10) $$\begin{eqnarray}\langle \,\hat{\boldsymbol{y}}_{1},\hat{\boldsymbol{y}}_{2}\rangle _{y}=\int _{\unicode[STIX]{x1D6FA}}\hat{\boldsymbol{y}}_{2}^{\ast }(\boldsymbol{x},f)\unicode[STIX]{x1D652}_{y}(\boldsymbol{x})\hat{\boldsymbol{y}}_{1}(\boldsymbol{x},f)\,\text{d}\boldsymbol{x}\end{eqnarray}$$

and

(5.11) $$\begin{eqnarray}\langle \hat{\boldsymbol{\unicode[STIX]{x1D702}}}_{1},\hat{\boldsymbol{\unicode[STIX]{x1D702}}}_{2}\rangle _{\unicode[STIX]{x1D702}}=\int _{\unicode[STIX]{x1D6FA}}\hat{\boldsymbol{\unicode[STIX]{x1D702}}}_{2}^{\ast }(\boldsymbol{x},f)\unicode[STIX]{x1D652}_{\unicode[STIX]{x1D702}}(\boldsymbol{x})\hat{\boldsymbol{\unicode[STIX]{x1D702}}}_{1}(\boldsymbol{x},f)\,\text{d}\boldsymbol{x}\end{eqnarray}$$

are inner products on the output and input spaces respectively.

Using (5.7) and (5.8), the input and output are related as

(5.12) $$\begin{eqnarray}\hat{\boldsymbol{y}}(\boldsymbol{x},f)={\mathcal{R}}(\boldsymbol{x},f)\hat{\boldsymbol{\unicode[STIX]{x1D702}}}(\boldsymbol{x},f),\end{eqnarray}$$

where

(5.13) $$\begin{eqnarray}{\mathcal{R}}(\boldsymbol{x},f)={\mathcal{C}}(\boldsymbol{x})(\text{i}2\unicode[STIX]{x03C0}f{\mathcal{I}}-{\mathcal{A}}(\boldsymbol{x}))^{-1}{\mathcal{B}}(\boldsymbol{x})\end{eqnarray}$$

is termed the resolvent operator. Note that $(\text{i}2\unicode[STIX]{x03C0}f{\mathcal{I}}-{\mathcal{A}})$ is invertible for all $f$ if ${\mathcal{A}}$ is stable; i.e. the real part of all of its eigenvalues is negative.

The optimal inputs and outputs are defined by the Schmidt decomposition of the resolvent operator

(5.14) $$\begin{eqnarray}{\mathcal{R}}(\boldsymbol{x},f)=\mathop{\sum }_{j=1}^{\infty }\unicode[STIX]{x1D70E}_{j}(f)\hat{\boldsymbol{u}}_{j}(\boldsymbol{x},f)\otimes \hat{\boldsymbol{v}}_{j}(\boldsymbol{x},f),\end{eqnarray}$$

defined in terms of the inner products (5.10) and (5.11). The input modes $\hat{\boldsymbol{v}}_{j}$ and output modes $\hat{\boldsymbol{u}}_{j}$ provide complete bases for their respective spaces and are orthogonal in their respective inner products (i.e. $\langle \hat{\boldsymbol{u}}_{j},\hat{\boldsymbol{u}}_{k}\rangle _{y}=\langle \hat{\boldsymbol{v}}_{j},\hat{\boldsymbol{v}}_{k}\rangle _{\unicode[STIX]{x1D702}}=\unicode[STIX]{x1D6FF}_{jk}$ ). The modes are ordered by their singular value $\unicode[STIX]{x1D70E}_{j}$ , and the gain between $\hat{\boldsymbol{v}}_{j}$ and $\hat{\boldsymbol{u}}_{j}$ is $\unicode[STIX]{x1D70E}_{j}^{2}$ . The resolvent operator is said to be low rank if the magnitude of the singular values falls off rapidly with increasing  $j$ .

Since the bases are complete, the output can be expanded as

(5.15) $$\begin{eqnarray}\hat{\boldsymbol{y}}(\boldsymbol{x},f)=\mathop{\sum }_{j=1}^{\infty }\hat{\boldsymbol{u}}_{j}(\boldsymbol{x},f)\unicode[STIX]{x1D70E}_{j}(f)\unicode[STIX]{x1D6FD}_{j}(f),\end{eqnarray}$$

where

(5.16) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{j}(f)=\langle \hat{\boldsymbol{\unicode[STIX]{x1D702}}}(\boldsymbol{x},f),\hat{\boldsymbol{v}}_{j}(\boldsymbol{x},f)\rangle _{\unicode[STIX]{x1D702}}\end{eqnarray}$$

is the projection of the forcing onto the $j$ th input mode (McKeon & Sharma Reference McKeon and Sharma2010). In some flows, reduced-order models have been obtained by retaining a limited number of terms in the expansion when the resolvent operator is low rank (e.g. Beneddine et al. Reference Beneddine, Sipp, Arnault, Dandois and Lesshafft2016; Gómez et al. Reference Gómez, Blackburn, Rudman, Sharma and McKeon2016b ).

5.2 Spectral densities in terms of resolvent modes

A connection between resolvent analysis and SPOD can be obtained by writing the cross-spectral density tensor in terms of resolvent modes. To do so, it is convenient to write the cross-spectral density tensor as

(5.17) $$\begin{eqnarray}\unicode[STIX]{x1D64E}_{yy}(\boldsymbol{x},\boldsymbol{x}^{\prime },f)=E\{\,\hat{\boldsymbol{y}}(\boldsymbol{x},f)\hat{\boldsymbol{y}}^{\ast }(\boldsymbol{x}^{\prime },f)\},\end{eqnarray}$$

which can be proven to be equivalent to the previous definition given by (2.16) (Bendat & Piersol Reference Bendat and Piersol2000). Inserting the resolvent-mode expansion (5.15) leads to the expression

(5.18) $$\begin{eqnarray}\unicode[STIX]{x1D64E}_{yy}(\boldsymbol{x},\boldsymbol{x}^{\prime },f)=\mathop{\sum }_{j=1}^{\infty }\mathop{\sum }_{k=1}^{\infty }\hat{\boldsymbol{u}}_{j}(\boldsymbol{x},f)\hat{\boldsymbol{u}}_{k}^{\ast }(\boldsymbol{x}^{\prime },f)\unicode[STIX]{x1D70E}_{j}(f)\unicode[STIX]{x1D70E}_{k}(f)\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}(f),\end{eqnarray}$$

where

(5.19) $$\begin{eqnarray}\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}(f)=E\{\unicode[STIX]{x1D6FD}_{j}(f)\unicode[STIX]{x1D6FD}_{k}^{\ast }(f)\}\end{eqnarray}$$

is the scalar-valued cross-spectral density between the $j$ th and $k$ th expansion coefficients. In obtaining (5.19), the output resolvent modes and singular values were moved outside of the expectation operator, which is permitted because they are deterministic quantities, depending only on the resolvent operator and the inner products. On the other hand, the expansion coefficients depend on the forcing term $\hat{\boldsymbol{\unicode[STIX]{x1D702}}}$ , which is stochastic due to the random nature of turbulent flows, so the expansion coefficients must be described by their cross-spectral density. This fact is central to the remainder of this section.

A connection between SPOD and resolvent analysis is obtained by choosing the SPOD quantity of interest $\boldsymbol{q}$ and the resolvent output $\boldsymbol{y}$ to be the same, with the same inner products, $\langle \cdot ,\cdot \rangle _{x}=\langle \cdot ,\cdot \rangle _{y}$ . Then, the cross-spectral density tensors $\unicode[STIX]{x1D64E}$ and $\unicode[STIX]{x1D64E}_{yy}$ are identical and the SPOD expansion (2.21) and resolvent-mode expansion (5.18) can be equated,

(5.20a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64E}_{yy}(\boldsymbol{x},\boldsymbol{x}^{\prime },f) & = & \displaystyle \mathop{\sum }_{j=1}^{\infty }\unicode[STIX]{x1D706}_{j}(f)\unicode[STIX]{x1D74D}_{j}(\boldsymbol{x},f)\unicode[STIX]{x1D74D}_{j}^{\ast }(\boldsymbol{x}^{\prime },f)\end{eqnarray}$$

(5.20b ) $$\begin{eqnarray}\displaystyle & = & \displaystyle \mathop{\sum }_{j=1}^{\infty }\mathop{\sum }_{k=1}^{\infty }\hat{\boldsymbol{u}}_{j}(\boldsymbol{x},f)\hat{\boldsymbol{u}}_{k}^{\ast }(\boldsymbol{x}^{\prime },f)\unicode[STIX]{x1D70E}_{j}(f)\unicode[STIX]{x1D70E}_{k}(f)\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}(f).\end{eqnarray}$$

5.3 Special case: uncorrelated expansion coefficients

The resolvent-mode expansion of the cross-spectral density tensor can be simplified in the special case in which the $\unicode[STIX]{x1D6FD}_{j}$ expansion coefficients are uncorrelated from one another, which is expressed as $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}(f)=\unicode[STIX]{x1D707}_{j}(f)\unicode[STIX]{x1D6FF}_{jk}$ . Substituting this into (5.20b ) leads to the simplified expression

(5.21a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64E}_{yy}(\boldsymbol{x},\boldsymbol{x}^{\prime },f) & = & \displaystyle \mathop{\sum }_{j=1}^{\infty }\unicode[STIX]{x1D74D}_{j}(\boldsymbol{x},f)\unicode[STIX]{x1D74D}_{j}^{\ast }(\boldsymbol{x}^{\prime },f)\unicode[STIX]{x1D706}_{j}(f)\end{eqnarray}$$

(5.21b ) $$\begin{eqnarray}\displaystyle & = & \displaystyle \mathop{\sum }_{j=1}^{\infty }\hat{\boldsymbol{u}}_{j}(\boldsymbol{x},f)\hat{\boldsymbol{u}}_{j}^{\ast }(\boldsymbol{x}^{\prime },f)\unicode[STIX]{x1D70E}_{j}^{2}(f)\unicode[STIX]{x1D707}_{j}(f).\end{eqnarray}$$

$\unicode[STIX]{x1D70E}_{j}^{2}\unicode[STIX]{x1D707}_{j}$

$\unicode[STIX]{x1D707}_{j}=1$

$j$

$\unicode[STIX]{x1D70E}_{j}^{2}(f)=\unicode[STIX]{x1D706}_{j}(f)$

$\hat{\boldsymbol{u}}_{j}(\boldsymbol{x},f)=\unicode[STIX]{x1D74D}_{j}(\boldsymbol{x},f)$

The conditions for which the expansion coefficients are uncorrelated can be elucidated by manipulating the definition of $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}$ using properties of inner products,

(5.22) $$\begin{eqnarray}\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}=E\{\langle \hat{\boldsymbol{\unicode[STIX]{x1D702}}},\hat{\boldsymbol{v}}_{j}\rangle _{\unicode[STIX]{x1D702}}\langle \hat{\boldsymbol{\unicode[STIX]{x1D702}}},\hat{\boldsymbol{v}}_{k}\rangle _{\unicode[STIX]{x1D702}}^{\ast }\}=E\{\langle \langle \hat{\boldsymbol{\unicode[STIX]{x1D702}}},\hat{\boldsymbol{v}}_{j}\rangle _{\unicode[STIX]{x1D702}}\hat{\boldsymbol{\unicode[STIX]{x1D702}}}^{\ast },\hat{\boldsymbol{v}}_{k}\rangle _{\unicode[STIX]{x1D702}}^{\ast }\}=\langle \langle \unicode[STIX]{x1D64E}_{\unicode[STIX]{x1D702}\unicode[STIX]{x1D702}},\hat{\boldsymbol{v}}_{j}\rangle _{\unicode[STIX]{x1D702}}^{\ast },\hat{\boldsymbol{v}}_{k}\rangle _{\unicode[STIX]{x1D702}}^{\ast },\end{eqnarray}$$

where $\unicode[STIX]{x1D64E}_{\unicode[STIX]{x1D702}\unicode[STIX]{x1D702}}(\boldsymbol{x},\boldsymbol{x}^{\prime },f)=E\{\hat{\boldsymbol{\unicode[STIX]{x1D702}}}(\boldsymbol{x},f)\hat{\boldsymbol{\unicode[STIX]{x1D702}}}^{\ast }(\boldsymbol{x}^{\prime },f)\}$ is the cross-spectral density tensor of the input $\boldsymbol{\unicode[STIX]{x1D702}}$ . Recalling the orthogonality of the resolvent output modes, the last form of (5.22) shows that $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}=\unicode[STIX]{x1D707}_{j}\unicode[STIX]{x1D6FF}_{jk}$ only if $\langle \unicode[STIX]{x1D64E}_{\unicode[STIX]{x1D702}\unicode[STIX]{x1D702}}(\boldsymbol{x},\boldsymbol{x}^{\prime },f),\hat{\boldsymbol{v}}_{j}(\boldsymbol{x}^{\prime },f)\rangle _{\unicode[STIX]{x1D702}}^{\ast }=\unicode[STIX]{x1D707}_{j}\hat{\boldsymbol{v}}_{j}(\boldsymbol{x},f)$ . Using the definition of the inner product, this can be written as

(5.23) $$\begin{eqnarray}\int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D64E}_{\unicode[STIX]{x1D702}\unicode[STIX]{x1D702}}(\boldsymbol{x},\boldsymbol{x}^{\prime },f^{\prime })\unicode[STIX]{x1D652}_{\unicode[STIX]{x1D702}}(\boldsymbol{x}^{\prime })\boldsymbol{v}_{j}(\boldsymbol{x}^{\prime },f^{\prime })\,\text{d}\boldsymbol{x}^{\prime }=\unicode[STIX]{x1D707}(f^{\prime })\boldsymbol{v}_{j}(\boldsymbol{x},f^{\prime }).\end{eqnarray}$$

This is precisely the eigenvalue problem that defines the SPOD modes of the input. Thus, uncorrelated expansion coefficients imply that the resolvent input modes are the same as the SPOD modes of the input (if there are repeated eigenvalues/singular values, then the corresponding modes must span the same subspace). Conversely, if the SPOD modes of the input are the same as the input resolvent modes, inserting the SPOD expansion of $\unicode[STIX]{x1D64E}_{\unicode[STIX]{x1D702}\unicode[STIX]{x1D702}}$ into (5.22) shows that the expansion coefficients are uncorrelated. Therefore, the expansion coefficients are uncorrelated if and only if the resolvent input modes correspond exactly with the SPOD modes of the input.

On inspection of (5.23), the stronger condition that $\unicode[STIX]{x1D707}_{j}=1$ for every $j$ requires that

(5.24) $$\begin{eqnarray}\unicode[STIX]{x1D64E}_{\unicode[STIX]{x1D702}\unicode[STIX]{x1D702}}(\boldsymbol{x},\boldsymbol{x}^{\prime },f)\unicode[STIX]{x1D652}_{\unicode[STIX]{x1D702}}(\boldsymbol{x}^{\prime })=\unicode[STIX]{x1D644}\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}^{\prime }).\end{eqnarray}$$

When $\unicode[STIX]{x1D652}_{\unicode[STIX]{x1D702}}=\unicode[STIX]{x1D644}$ , this condition takes on physical significance: $\unicode[STIX]{x1D64E}_{\unicode[STIX]{x1D702}\unicode[STIX]{x1D702}}(\boldsymbol{x},\boldsymbol{x}^{\prime },f)=\unicode[STIX]{x1D644}\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}^{\prime })$ corresponds to an input that is completely uncorrelated in space and time and has unit amplitude everywhere, i.e. unit-amplitude white noise. While the nonlinear forcing terms in any real flow are unlikely to be white, this approximation has been shown to be reasonable in some flows and is frequently used for the construction of low-order models (Farrell & Ioannou Reference Farrell and Ioannou1993, Reference Farrell and Ioannou1996, Reference Farrell and Ioannou2001), including resolvent-based models (Jovanović & Bamieh Reference Jovanović and Bamieh2001; Bagheri et al. Reference Bagheri, Henningson, Hoepffner and Schmid2009; Sipp et al. Reference Sipp, Marquet, Meliga and Barbagallo2010; Moarref & Jovanović Reference Moarref and Jovanović2012; Dergham et al. Reference Dergham, Sipp and Robinet2013). In these models, then, resolvent modes can be understood as approximations of SPOD modes. The equivalence of SPOD and resolvent modes in the case of white-noise forcing has been recently pointed out by multiple authors working in the area of jet-noise modelling (Towne et al. Reference Towne, Colonius, Jordan, Cavalieri and Brès2015; Semeraro et al. Reference Semeraro, Jaunet, Jordan, Cavalieri and Lesshafft2016a ; Towne, Brès & Lele Reference Towne, Brès and Lele2016).

5.4 General case: flow reconstruction with correlated expansion coefficients

The nonlinear forcing terms in real flows are not white (Towne, Brès & Lele Reference Towne, Brès and Lele2017; Zare, Jovanović & Georgiou Reference Zare, Jovanović and Georgiou2017), which leads to correlated resolvent-mode expansion coefficients. In this case, a more general relationship between SPOD and resolvent modes can be derived that provides insight into how the resolvent-mode expansion coefficients can be chosen to optimally reconstruct the flow.

A number of strategies for choosing the resolvent-mode expansion coefficients have been proposed and employed in the last few years. Perhaps the most straightforward approach, used for example by Jeun, Nichols & Jovanović (Reference Jeun, Nichols and Jovanović2016), is to use simulation data to compute $\boldsymbol{\unicode[STIX]{x1D702}}$ , estimate $\hat{\boldsymbol{\unicode[STIX]{x1D702}}}$ using a single DFT of the time series, and calculate the expansion coefficients using their definition in (5.16). However, this method suffers from the uncertainty inherent in the DFT, as discussed earlier. Methods using converged statistics of the flow must be used to eliminate this uncertainty.

Moarref et al. (Reference Moarref, Jovanović, Tropp, Sharma and McKeon2014) formulated a convex optimization problem to find a vector of expansion coefficients that optimally reproduce experimentally measured streamwise energy spectra, averaged in the other directions, for a turbulent channel flow. The coefficients that were optimal in the sense they defined led to qualitative agreement with the statistics, but notably the error did not tend towards zero as the number of resolvent modes was increased.

Gómez et al. (Reference Gómez, Blackburn, Rudman, Sharma and McKeon2016a ) derived an equation for the power spectral densities of the expansion coefficients, which define their magnitude. These are obtained by projecting the time-varying forcing term $\boldsymbol{\unicode[STIX]{x1D702}}$ onto the resolvent input basis for a range of frequencies, and computing the power spectral densities of these terms. In a separate paper (Gómez et al. Reference Gómez, Blackburn, Rudman, Sharma and McKeon2016b ), the same authors suggested a somewhat different method using a time series of data at one location in the flow. In this approach, the resolvent expansion is inverse Fourier transformed back into the time domain and the coefficients for all frequencies are simultaneously determined by matching the time series in a least-squares sense.

Beneddine et al. (Reference Beneddine, Sipp, Arnault, Dandois and Lesshafft2016) used two different approaches for finding the coefficient for (only) the leading resolvent mode at each frequency for the flow over a backward-facing step. In the first approach, the amplitudes of the expansion coefficients were chosen so that the rank-one resolvent-mode reconstruction matched the power spectral density of the spanwise-averaged streamwise velocity at a single location in the flow. In the second approach, the expansion coefficient for the leading mode at each frequency was chosen so that the rank-one resolvent expansion matched the leading SPOD mode, again at a single location. They found this second method to provide improved estimates of the power spectral density at other spatial locations in the flow. Their results were further improved by attempting to match the SPOD mode at several locations using a least-squares approach.

All of these methods assume that each expansion coefficient can be described by a single deterministic amplitude and phase. In contrast, (5.18) shows that the cross-spectral densities of the expansion coefficients are required to properly reconstruct the second-order statistics of the flow. Notably, this includes the power spectral density, which is obtained by setting $\boldsymbol{x}^{\prime }=\boldsymbol{x}$ . Treating the expansion coefficients as deterministic quantities imposes a fundamental limitation on the quality of the flow approximation that can be achieved by the resolvent model. For clarity, we will use the symbol $b_{j}$ in place of $\unicode[STIX]{x1D6FD}_{j}$ to denote these deterministic expansion coefficients. Regardless of how the $b_{j}$ values are chosen, the approximation of $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}$ implied by these coefficients is $\unicode[STIX]{x1D61A}_{b_{j}b_{k}}=b_{j}b_{k}^{\ast }$ . Making this substitution in (5.18), the resulting cross-spectral density tensor can be factored into the form

(5.25) $$\begin{eqnarray}\unicode[STIX]{x1D64E}_{yy}(\boldsymbol{x},\boldsymbol{x}^{\prime },f)=\left(\mathop{\sum }_{j=1}^{\infty }\hat{\boldsymbol{u}}_{j}(\boldsymbol{x},f)\unicode[STIX]{x1D70E}_{j}(f)b_{j}(f)\right)\left(\mathop{\sum }_{k=1}^{\infty }\hat{\boldsymbol{u}}_{k}(\boldsymbol{x}^{\prime },f)\unicode[STIX]{x1D70E}_{k}(f)b_{k}(f)\right)^{\ast }.\end{eqnarray}$$

Since the right-hand side is the product of two vectors, this approximation of the cross-spectral density tensor has rank equal to one, whereas the correct cross-spectral density tensor has rank equal to the number of non-zero SPOD eigenvalues at each frequency. This limits the quality of the flow approximation that can be achieved with the resolvent model using deterministic expansion coefficients, no matter how their amplitudes and phases are chosen. Specifically, we will see that convergent approximations of the power spectral and cross-spectral densities are not possible when the statistical nature of the expansion coefficients is not respected.

A few specific choices of $b_{j}$ are worth a closer look. First, we consider the case where the $b_{j}$ values are chosen to have the correct magnitude, so that their power spectral densities $\unicode[STIX]{x1D61A}_{b_{j}b_{j}}=|b_{j}|^{2}$ match the correct values $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{j}}$ . This could be accomplished, for example, using the approach of Gómez et al. (Reference Gómez, Blackburn, Rudman, Sharma and McKeon2016a ). Equation (5.18) shows that the cross-spectral and power spectral densities depend also on off-diagonal terms $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}$ . These will be replaced by the values $\unicode[STIX]{x1D61A}_{b_{j}b_{k}}=b_{j}b_{k}^{\ast }$ , and there is no reason to expect this to provide a good approximation. For example, in the case of uncorrelated expansion coefficients, the off-diagonal terms of $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}$ are zero, but $b_{j}b_{k}^{\ast }$ is non-zero unless either mode $j$ or mode $k$ has zero magnitude. Thus, even though the power spectral densities of the expansion coefficients are correct, the power spectral densities of the output will be incorrect. To make the error explicit, (5.25) can be manipulated into the form

(5.26) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64E}_{yy}(\boldsymbol{x},\boldsymbol{x}^{\prime },f) & = & \displaystyle \mathop{\sum }_{j=1}^{\infty }\hat{\boldsymbol{u}}_{j}(\boldsymbol{x},f)\hat{\boldsymbol{u}}_{j}^{\ast }(\boldsymbol{x}^{\prime },f)\unicode[STIX]{x1D70E}_{j}^{2}(f)|b_{j}(f)|^{2}\nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{j=1}^{\infty }\mathop{\sum }_{\substack{ k=1 \\ k\neq j}}^{\infty }\hat{\boldsymbol{u}}_{j}(\boldsymbol{x},f)\hat{\boldsymbol{u}}_{k}^{\ast }(\boldsymbol{x}^{\prime },f)\unicode[STIX]{x1D70E}_{j}(f)\unicode[STIX]{x1D70E}_{k}(f)b_{j}(f)b_{k}^{\ast }(f).\end{eqnarray}$$

Using the power-spectral-density-based expansion coefficients, the first term in (5.26) is correct, while the second term is incorrect no matter how many terms are included in the expansion.

An alternative approach would be to choose the $b_{j}$ coefficients such that the resolvent-mode expansion optimally matches the data, as defined by the inner product $\langle \cdot ,\cdot \rangle _{y}$ , in the limit of high $N_{r}$ . Spectral POD theory tells us that the optimal representation of the data under the rank-one limitation is given by the leading SPOD mode. Therefore, the optimal $b_{j}$ values are those that reconstruct the leading SPOD mode at each frequency. These values can be obtained by directly projecting the leading SPOD mode onto the resolvent modes, giving

(5.27) $$\begin{eqnarray}b_{j}(f)=\frac{\sqrt{\unicode[STIX]{x1D706}_{1}(f)}}{\unicode[STIX]{x1D70E}_{j}(f)}\langle \unicode[STIX]{x1D74D}_{1},\hat{\boldsymbol{u}}_{j}\rangle _{y}.\end{eqnarray}$$

By construction, these expansion coefficients guarantee that the first SPOD mode at each frequency is recovered when enough terms are retained in the resolvent-mode expansion. Even though this is the optimal expansion possible using a deterministic vector of expansion coefficients, it will provide a good approximation of the flow only if $\unicode[STIX]{x1D706}_{1}(f)\gg \unicode[STIX]{x1D706}_{2}(f)$ over the range of relevant frequencies. Therefore, the quality of the resolvent-mode expansion using optimal (or any other) deterministic coefficients is dependent first and foremost on the low-rank nature of the cross-spectral density tensor rather than of the resolvent operator.

In contrast, convergent approximations of the flow statistics can be achieved if the statistical nature of the resolvent-mode expansion coefficients is incorporated at the outset. Accordingly, we propose a practical method for approximating the required cross-spectral densities $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}$ using the leading SPOD modes. The resolvent-mode expansion coefficients are related to SPOD modes by comparing the SPOD and resolvent expansions of the cross-spectral density tensor given in (5.20). By taking two successive inner products of both expansions with respect to $\hat{\boldsymbol{u}}_{j}$ and $\hat{\boldsymbol{u}}_{k}$ and dividing by $\unicode[STIX]{x1D70E}_{j}(f)$ and $\unicode[STIX]{x1D70E}_{k}(f)$ , we obtain the expression

(5.28) $$\begin{eqnarray}\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}=\mathop{\sum }_{n=1}^{\infty }\frac{\unicode[STIX]{x1D706}_{n}(f)}{\unicode[STIX]{x1D70E}_{j}(f)\unicode[STIX]{x1D70E}_{k}(f)}\unicode[STIX]{x1D6FE}_{nj}(f)\unicode[STIX]{x1D6FE}_{nk}^{\ast }(f),\end{eqnarray}$$

where

(5.29) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{jk}(f)=\langle \unicode[STIX]{x1D74D}_{j},\hat{\boldsymbol{u}}_{k}\rangle _{y}\end{eqnarray}$$

is the projection between the $j$ th SPOD mode and the $k$ th resolvent output mode.

Using the first $N_{s}$ terms in (5.28) to approximate $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}$ and a resolvent-mode expansion with $N_{r}$ terms leads to the approximation of the second-order flow statistics,

(5.30) $$\begin{eqnarray}\unicode[STIX]{x1D64E}_{yy}(\boldsymbol{x},\boldsymbol{x}^{\prime },f)\approx \mathop{\sum }_{n=1}^{N_{s}}\unicode[STIX]{x1D706}_{n}(f)\left(\mathop{\sum }_{j=1}^{N_{r}}\unicode[STIX]{x1D6FE}_{nj}(f)\hat{\boldsymbol{u}}_{j}(\boldsymbol{x},f)\right)\left(\mathop{\sum }_{k=1}^{N_{r}}\unicode[STIX]{x1D6FE}_{nk}(f)\hat{\boldsymbol{u}}_{k}(\boldsymbol{x}^{\prime },f)\right)^{\ast }.\end{eqnarray}$$

From the definition of $\unicode[STIX]{x1D6FE}_{jk}$ and the fact that the input resolvent modes are complete, we have that

(5.31) $$\begin{eqnarray}\mathop{\sum }_{j=1}^{N_{r}}\unicode[STIX]{x1D6FE}_{nj}(f)\hat{\boldsymbol{u}}_{j}(\boldsymbol{x},f)\rightarrow \unicode[STIX]{x1D74D}_{n}(\boldsymbol{x},f)\quad \text{as}\;N_{r}\rightarrow \infty .\end{eqnarray}$$

Therefore, using the first $N_{s}$ SPOD modes to approximate $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}$ leads to a resolvent-mode expansion capable of capturing the first $N_{s}$ SPOD modes, if enough resolvent modes are retained. The approximation is therefore convergent as $N_{s}$ and $N_{r}$ are increased, and if the SPOD spectrum falls off rapidly, accurate approximations of the flow are possible using only a few SPOD modes to compute the expansion coefficients. If a sufficiently high value of $N_{s}$ is used, then the quality of the approximation that can be obtained with a limited number of resolvent modes is determined by the convergence of the limit in (5.31). It is notable that the resolvent gains do not appear explicitly in (5.30) or (5.31). Instead, the efficacy of the resolvent expansion is determined by the efficiency of the leading resolvent modes in providing a basis for the leading SPOD modes.

6.1 Complex Ginzburg–Landau equation

Our first example consists of the linearized complex Ginzburg–Landau equation, which is frequently used as a model for instabilities in spatially evolving flows. The equation can be written in the form of (5.3) with

(6.1) $$\begin{eqnarray}{\mathcal{A}}=-\unicode[STIX]{x1D708}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}x}+\unicode[STIX]{x1D6FE}\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}+\unicode[STIX]{x1D707}(x).\end{eqnarray}$$

The spatial dependence of the solution is a consequence of the parameter $\unicode[STIX]{x1D707}(x)$ , for which we adopt the quadratic form

(6.2) $$\begin{eqnarray}\unicode[STIX]{x1D707}(x)=(\unicode[STIX]{x1D707}_{0}-c_{\unicode[STIX]{x1D707}}^{2})+\frac{\unicode[STIX]{x1D707}_{2}}{2}x^{2}\end{eqnarray}$$

used previously by several authors (Hunt & Crighton Reference Hunt and Crighton1991; Bagheri et al. Reference Bagheri, Henningson, Hoepffner and Schmid2009; Chen & Rowley Reference Chen and Rowley2011). We take $\unicode[STIX]{x1D707}_{0}=0.23$ , and all other parameters in (6.1) and (6.2) are set to the values used by Bagheri et al. (Reference Bagheri, Henningson, Hoepffner and Schmid2009). The resulting model is globally stable (all eigenvalues of ${\mathcal{A}}$ have negative real part) but is susceptible to non-modal growth due to the non-normality of ${\mathcal{A}}$ , which within the context of resolvent analysis leads to gains much larger than one. The input and output operators and inner-product matrices ( ${\mathcal{B}}$ , ${\mathcal{C}}$ , $\unicode[STIX]{x1D652}_{y}$ and $\unicode[STIX]{x1D652}_{\unicode[STIX]{x1D702}}$ respectively) are all set to unity. The value of $\unicode[STIX]{x1D707}_{0}$ was chosen so that the gain of the leading resolvent mode at its peak frequency is 100 times larger than the gain of the second mode, which is a typical value for real flows.

The equations are discretized with a pseudo-spectral approach using Hermite polynomials, as in Bagheri et al. (Reference Bagheri, Henningson, Hoepffner and Schmid2009) and Chen & Rowley (Reference Chen and Rowley2011). The collocation points correspond to the roots of the first $N$ Hermite polynomials; following Bagheri et al. (Reference Bagheri, Henningson, Hoepffner and Schmid2009), we use $N=220$ . This leads to a computational domain $x\in [-85,85]$ , which is large enough to mimic an infinite domain. Since the collocation points are unevenly distributed, the discretized form of the inner product would contain matrices different from the identity. To avoid this, we define the input $\unicode[STIX]{x1D702}$ and output $y$ on a different uniformly spaced grid. This is accomplished by setting the discretized forms of the input and output matrices $B$ and $C$ so that they represent an interpolation from the uniform grid to the Hermite grid and from the Hermite grid to the uniform grid respectively. On the uniform input/output grid, the discrete inner-product matrices $\mathbf{W}_{\mathbf{y}}$ and $\mathbf{W}_{\boldsymbol{\unicode[STIX]{x1D702}}}$ are the identity matrix.

To generate data for our analysis, the discretized equations are stochastically excited in the time domain using forcing terms with prescribed statistics. In the following two subsections, we consider cases of white-noise and correlated forcing respectively.

6.1.1 White-noise forcing

We begin by forcing the linearized Ginzburg–Landau equations with band-limited white noise that is spatially uncorrelated. This forcing is realized by setting the value at each grid point and at each discrete time instance to a random complex number with uniformly distributed phase (between 0 and $2\unicode[STIX]{x03C0}$ ) and normally distributed amplitude with unit variance. This signal is then low-pass filtered in time using a 10th-order finite-impulse-response filter with a cutoff frequency equal to 60 % of the Nyquist frequency. Additionally, the forcing is spatially limited to the interior portion of the domain using an exponential envelope of the form $\exp [(x/L)^{p}]$ with $L=60$ and $p=10$ .

The equations are integrated using a fourth-order embedded Runge–Kutta method (Shampine & Reichelt Reference Shampine and Reichelt1997), and a total of 40 000 snapshots of the solution are collected with spacing $\unicode[STIX]{x0394}t=0.5$ , leading to a Nyquist frequency of $f_{Nyquist}=1$ . A large number of snapshots are required for this problem because of the slow statistical convergence of the white-noise forcing; problems with correlated forcing, as is typical in real flows, tend to require fewer snapshots. Spectra and SPOD modes are computed as outlined in § 3 with a Hann window function, 75 % overlap and $N_{f}=384$ .

Figure 2. Comparison of the first four resolvent gains $\unicode[STIX]{x1D70E}_{j}^{2}$ (lines) and SPOD eigenvalues $\unicode[STIX]{x1D706}_{j}$ (symbols) as a function of frequency for the Ginzburg–Landau model forced with white noise: (○, blue) mode 1; (▫, green) mode 2; (▵, red) mode 3; (▿, cyan) mode 4. For each mode, every other SPOD eigenvalue has been omitted from the plot to improve readability.

Since the forcing is white and the input inner product has unit weight, the conditions are met under which resolvent and SPOD modes are theoretically equivalent. The first four SPOD eigenvalues and resolvent gains are shown in figure 2 as a function of frequency. For consistency with previous publications, frequencies are reported in terms of the angular frequency $\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}f$ . Overall, the SPOD eigenvalues agree with the resolvent gains, but two minor differences can be observed which are related to the finite level of convergence of the spectral estimates used to approximate the SPOD modes. First, the SPOD eigenvalues are not completely smooth as a function of frequency. This is indicative of the remaining uncertainty in the spectral estimates, which can be further reduced by averaging over more data ensembles. Second, there is a small overshoot in the value of the second SPOD eigenvalue at the peak frequency of the first mode. This overshoot is related to spectral leakage, which can be reduced by increasing the frequency resolution of the spectral estimates by increasing the length of the data ensembles. When estimating the cross-spectral density using a data set of finite length, increasing the number of ensembles reduces the length of each ensemble and vice versa. Therefore, decreasing one of the two types of errors tends to increase the other. The proper compromise will in general be problem-dependent. As a general rule, SPOD modes are more difficult to converge for flows that exhibit sharp spectral peaks.

Figure 3. Weighted mode shapes in $\unicode[STIX]{x1D714}$ – $x$ space for the Ginzburg–Landau equation forced with white noise. (a,d,g) Resolvent modes $\unicode[STIX]{x1D70E}_{j}(\unicode[STIX]{x1D714})|\hat{\boldsymbol{u}}_{j}(x,\unicode[STIX]{x1D714})|$ . (b,e,h) Spectral POD modes $\sqrt{\unicode[STIX]{x1D706}_{j}(\unicode[STIX]{x1D714})}|\unicode[STIX]{x1D74D}_{j}(x,\unicode[STIX]{x1D714})|$ . (c,f,i) Space-only POD modes $|\hat{a}_{j}(\unicode[STIX]{x1D714})|\,|\unicode[STIX]{x1D753}_{j}(x)|$ .

The structure of the first three resolvent and SPOD modes is depicted in figure 3 as a function of $x$ and $\unicode[STIX]{x1D714}$ (space-only POD results are also shown and will be discussed later). Specifically, we plot the weighted quantities $\unicode[STIX]{x1D70E}_{j}(\unicode[STIX]{x1D714})|\hat{\boldsymbol{u}}_{j}(x,\unicode[STIX]{x1D714})|$ and $\sqrt{\unicode[STIX]{x1D706}_{j}(\unicode[STIX]{x1D714})}|\unicode[STIX]{x1D74D}_{j}(x,\unicode[STIX]{x1D714})|$ to clearly show where each mode is active in $\unicode[STIX]{x1D714}$ – $x$ space. The contour levels in each plot are distributed between zero and the maximum value of the weighted mode.

Figure 4. Spectral POD/resolvent-mode projection coefficient (a) and resolvent-mode expansion-coefficient cross-spectral density (b) at the frequency of highest gain for the Ginzburg–Landau model forced with white noise. The first 12 SPOD and resolvent modes are nearly identical (diagonal $\unicode[STIX]{x1D6FE}_{jk}$ ) because the resolvent-mode expansion coefficients are almost completely uncorrelated (diagonal $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}$ ).

The first resolvent and SPOD modes are nearly indistinguishable. The suboptimal SPOD modes likewise closely mimic their corresponding resolvent modes, although the limited smoothness as a function of frequency is again evident. The match between the resolvent and SPOD modes can be quantified using the projection $\unicode[STIX]{x1D6FE}_{jk}$ , defined in (5.29). This is shown for the peak frequency of the leading resolvent mode, $\unicode[STIX]{x1D714}=-0.6$ , in figure 4(a). The diagonal ( $j=k$ ) entries dominate up to $j=k=12$ , indicating that the first 12 resolvent and SPOD modes are nearly the same. Finally, figure 4(b) shows the expansion-coefficient cross-spectral density tensor $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}$ for the same frequency, normalized by its maximum value. As expected, it is nearly diagonal, which is the root cause of the approximate equivalence of the resolvent and SPOD modes for this problem.

We also briefly compare the SPOD results with those obtained using space-only POD. Figure 3(c,f,i) shows the $\unicode[STIX]{x1D714}$ – $x$ distribution of the first three POD modes. The frequency distribution was obtained by computing the power spectral densities of the time-dependent expansion coefficients $a_{j}(t)$ using Welch’s method, which is similar to the method of Cammilleri et al. (Reference Cammilleri, Guéniat, Carlier, Pastur, Mémin, Lusseyran and Artana2013) except that we use proper spectral averaging rather than a single DFT to obtain frequency information. The quantity plotted is then $|\hat{a}_{j}(\unicode[STIX]{x1D714})|\,|\unicode[STIX]{x1D753}_{j}(x)|$ . As expected, each space-only POD mode represents a range of frequencies. The first space-only POD mode is relatively similar to the leading SPOD mode, but subsequent POD modes represent structures that are quite different. Two things stand out. First, since the spatial shape of each space-only POD mode does not depend on frequency, they are unable to capture the frequency-dependent shape of the resolvent modes that is captured by the SPOD modes. Second, the frequency-dependent magnitude of each suboptimal space-only POD mode exhibits a clear dip at the frequencies at which previous modes peak. Overall, the SPOD modes provide a better representation than the space-only POD modes of the underlying dynamics and non-normal growth mechanisms, which are completely described by the resolvent modes for this problem since the forcing is white.

Finally, we demonstrate the result from § 2.3 that space-only POD modes do not represent structures that are uncorrelated from one another in time. Figure 5 shows the correlation $\unicode[STIX]{x1D60A}_{a_{1}a_{j}}^{POD}(\unicode[STIX]{x1D70F})=E\{a_{1}(t)a_{j}^{\ast }(t+\unicode[STIX]{x1D70F})\}$ of the first POD expansion coefficient with the first three expansion coefficients ( $j=1,2,3$ ) as a function of the temporal lag $\unicode[STIX]{x1D70F}$ , scaled by $\sqrt{\unicode[STIX]{x1D706}_{1}\unicode[STIX]{x1D706}_{j}}$ . The solid and broken lines show the real part and magnitude of the correlations respectively. At zero time lag, the scaled autocorrelation $\unicode[STIX]{x1D60A}_{a_{1}a_{1}}^{POD}(0)/\unicode[STIX]{x1D706}_{1}$ is one and the cross-correlations $\unicode[STIX]{x1D60A}_{a_{1}a_{k}}^{POD}(0)$ are zero. Both of these values follow theoretically from (2.9). However, for $\unicode[STIX]{x1D70F}>0$ , the autocorrelation quickly drops and the cross-correlations become non-zero. Accordingly, the space-only POD modes do not represent structures that evolve coherently in time.

Figure 5. Temporal correlation of the first three space-only POD expansion coefficients with the first coefficient (see (2.24)) for the Ginzburg–Landau model forced with white noise. Real part of the correlation for (——, blue) mode 1; (– – –, green) mode 2; (– ⋅ – ⋅ –, red) mode 3. The light dotted lines show the magnitude of each correlation. Since the autocorrelation of the first mode decays and the cross-correlation with the suboptimal modes becomes non-zero for $\unicode[STIX]{x1D70F}>0$ , the space-only POD modes do not represent structures that evolve coherently in space and time.

6.1.2 Correlated forcing

Next, we force the linearized Ginzburg–Landau equations with forcing terms that are spatially correlated. This forcing is generated by convolving the same band-limited white-noise signal used previously with a kernel of the form

(6.3) $$\begin{eqnarray}g(x,x^{\prime })=\frac{1}{\sqrt{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D702}}}\exp \left[-\frac{1}{2}\left(\frac{x-x^{\prime }}{\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D702}}}\right)^{2}\right]\exp \left[\text{i}2\unicode[STIX]{x03C0}\frac{x-x^{\prime }}{\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D702}}}\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D702}}$ is the standard deviation of the envelope and $\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D702}}$ is the wavelength of the filter. This leads to a forcing that is white-in-time up to the cutoff frequency but that has non-zero spatial correlation in the form of (6.3) but with $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D702}}$ replaced with $\sqrt{2}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D702}}$ . This form of the forcing statistics is qualitatively similar to those of the nonlinear forcing terms in real flows, such as a turbulent jet (Towne et al. Reference Towne, Brès and Lele2017). We use $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D702}}=4$ and $\unicode[STIX]{x1D706}_{\unicode[STIX]{x1D702}}=20$ . The statistics of the response converge more rapidly in this case due to the correlated forcing, so we use a smaller number of snapshots, $N=10\,000$ , than in the uncorrelated case.

Figure 6. Spectral POD/resolvent projection coefficients $\unicode[STIX]{x1D6FE}_{jk}$ at two frequencies for the Ginzburg–Landau model forced by spatially correlated input. Since the tensors are not diagonal, the resolvent and SPOD modes are not the same in this case.

The correlation of the forcing leads to differences between the output resolvent and SPOD modes. This is clearly demonstrated by the $\unicode[STIX]{x1D6FE}_{jk}$ projection coefficients, which are shown for $\unicode[STIX]{x1D714}=-0.6$ and 0.4 in figure 6. The first of these frequencies is near the peak frequency of the leading resolvent mode where $\unicode[STIX]{x1D70E}_{1}^{2}\approx 560$ and $\unicode[STIX]{x1D70E}_{1}^{2}/\unicode[STIX]{x1D70E}_{2}^{2}\approx 100$ , while the second is away from the peak and $\unicode[STIX]{x1D70E}_{1}^{2}\approx 20$ and $\unicode[STIX]{x1D70E}_{1}^{2}/\unicode[STIX]{x1D70E}_{2}^{2}\approx 10$ . In both cases, the diagonal character of $\unicode[STIX]{x1D6FE}_{jk}$ observed for the white-noise forcing (figure 4 a) is no longer observed, indicating that the SPOD and resolvent modes are no longer the same.

Figure 7. Cross-spectral density of the resolvent-mode expansion coefficients for (a,b,c) $\unicode[STIX]{x1D714}=-0.6$ and (d,e,f) $\unicode[STIX]{x1D714}=0.4$ . (a,d) Correct values. (b,e) Values implied by deterministic coefficients based on a single long-time DFT of the solution. (c,f) Values implied by deterministic coefficients with correct power spectral density.

As discussed in § 5, this can be attributed to the fact that the correlated forcing produces correlated resolvent-mode expansion coefficients. The magnitude of the cross-spectral density of the expansion coefficients $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}$ is shown in figure 7(a,d) for the same two frequencies. The strictly diagonal correlation tensor observed for the white-noise forcing (figure 4 b) has been replaced by a more complicated banded structure, indicating substantial correlations between different expansion coefficients.

The expansion-coefficient correlations are not captured by standard approaches to computing the expansion coefficients that treat them as deterministic quantities. For example, figures 7(b,e) and 7(c,f) show the correlations implied using the DFT- and power spectral density (PSD)-based methods discussed in § 5.4, both of which lead to large errors away from the main correlation bands.

As indicated by (5.18), these errors in the expansion-coefficient cross-spectral density lead to errors in the resolvent-mode reconstruction of the output statistics. Figure 8 compares the true power spectral density of the output as a function of $\unicode[STIX]{x1D714}$ and $x$ with the power spectral densities obtained from resolvent-mode reconstructions of the flow with $N_{r}=1$ , 5, 10 and 30 modes and with the expansion coefficients specified in four different ways. Figure 8(a) shows the true power spectral density, which does not depend on $N_{r}$ . Figure 8(b) shows the power spectral density obtained using the correct values of $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}$ , computed using (5.28). Figures 8(c) and 8(d) show the power spectral densities obtained using the DFT- and PSD-based deterministic expansion coefficients respectively, and figure 8(e) shows results for the optimal $b_{j}$ values given by (5.27). The contour levels in all of the panels are the same and span six orders of magnitude, with the upper bound set to the maximum value of the true power spectral density.

Figure 8. Resolvent-mode reconstructions of the PSD of the solution of the Ginzburg–Landau equation for spatially correlated forcing using different expansion coefficients. Rows: the number of resolvent modes in the expansion increases from top to bottom. Columns: (a) True PSD (same in every row). (b) Reconstruction with correct statistical expansion coefficients. (c) Reconstruction with deterministic expansion coefficients based on a single long-time DFT. (d) Reconstruction with deterministic PSD-based expansion coefficients. (e) Reconstruction with optimal deterministic expansion coefficients.

All four resolvent-mode reconstructions produce similar power spectral densities for $N_{r}=1$ , but the DFT-based expansion is noticeably noisy due to the uncertainty implicit in this approach. As the number of resolvent modes included in the expansion increases, the power spectral density obtained using the proper $\unicode[STIX]{x1D61A}_{\unicode[STIX]{x1D6FD}_{j}\unicode[STIX]{x1D6FD}_{k}}$ values converges to the true power spectral density. In contrast, none of the three methods using deterministic expansion coefficients converge to the true PSD due to their incorrect implied correlations. For the DFT- and PSD-based methods, the reconstructed power spectral densities improve with increasing $N_{r}$ in low-energy regions of $\unicode[STIX]{x1D714}$ – $x$ space but do not improve in high-energy regions. In fact, the reconstruction of the high-energy regions actually worsens with increasing $N_{r}$ for the PSD-based expansion coefficients due to accumulation of the error terms from (5.26) as more modes are included in the expansion. On the other hand, the inclusion of more terms in the resolvent-mode expansion using the optimal $b_{j}$ values has less effect on the low-energy regions but does not degrade the approximation of the high-energy regions. This is because the resolvent-mode reconstruction using these expansion coefficients is converging towards the leading SPOD mode at each frequency, and these modes effectively represent high-energy regions by construction. The expansion with optimal $b_{j}$ values is well converged using 30 resolvent modes, so this is the best approximation of the second-order output statistics that can be obtained using deterministic expansion coefficients.

6.2 Turbulent jet

Our second example is a subsonic turbulent jet issued from a round convergent-straight nozzle. The Mach number, temperature ratio and Reynolds number of the jet are $M=U_{j}/c_{j}=0.4$ , $T_{j}/T_{\infty }=1$ and $Re=\unicode[STIX]{x1D70C}_{j}U_{j}D/\unicode[STIX]{x1D707}_{j}\approx 10^{6}$ respectively, where $U$ is the velocity, $c$ is the speed of sound, $T$ is the temperature, $\unicode[STIX]{x1D70C}$ is the density, $D$ is the nozzle diameter, $\unicode[STIX]{x1D707}$ is the dynamic viscosity, and the subscripts $j$ and $\infty$ indicate mean conditions at the nozzle exit and in the ambient far field respectively.

We use data from a large-eddy simulation (LES) as input for the various empirical modal decomposition methods discussed in this paper. The simulation was performed using the compressible flow solver ‘Charles’ developed at Cascade Technologies (Brès et al. Reference Brès, Ham, Nichols and Lele2017a ), which solves the spatially filtered compressible Navier–Stokes equations on unstructured grids using a finite-volume method and third-order Runge–Kutta time integration. The simulation used approximately sixteen million control volumes and was run for a duration of 2000 acoustic time units ( $tc_{\infty }/D$ ). The numerical methods, geometry and grid are identical to those used in a previous simulation at a higher Mach number, which was validated using double-blind comparisons with an extensive set of experimental measurements (Brès et al. Reference Brès, Jaunet, Le Rallic, Jordan, Colonius and Lele2015, Reference Brès, Jaunet, Le Rallic, Jordan, Towne, Schmidt, Colonius, Cavalieri and Lele2016, Reference Brès, Jordan, Jaunet, Le Rallic, Cavalieri, Towne, Lele, Colonius and Schmidt2017b ).

The available LES database consists of 10 000 snapshots of the jet (velocities, density and pressure) sampled every 0.2 acoustic time units on a structured cylindrical output grid that approximately mirrors the underlying LES resolution and extends a distance of 30 jet diameters in the streamwise direction and six jet diameters in the radial direction. The azimuthal vorticity magnitude and pressure fluctuations from one of the snapshots are shown in figure 9. These snapshots are used to compute SPOD, space-only POD, DMD and DFT modes. Since the azimuthal coordinate is periodic, the snapshots can be decomposed into azimuthal Fourier modes, and each azimuthal mode can treated independently within the context of the modal decompositions. For the sake of brevity, we exclusively present results for the axisymmetric component, and all modes are visualized using the pressure field.

Figure 9. Instantaneous snapshot of the Mach 0.4 turbulent jet. Greyscale: pressure fluctuations. Colour: vorticity magnitude.

6.2.1 Spectral POD modes

First, we study the SPOD modes of the jet. Spectral POD has been applied to jets before using data from both experiments (Glauser et al. Reference Glauser, Leib and George1987; Arndt et al. Reference Arndt, Long and Glauser1997; Citriniti & George Reference Citriniti and George2000; Gudmundsson & Colonius Reference Gudmundsson and Colonius2011; Sinha et al. Reference Sinha, Rodriguez, Bres and Colonius2014; Semeraro et al. Reference Semeraro, Jaunet, Jordan, Cavalieri and Lesshafft2016a ) and simulations (Towne et al. Reference Towne, Colonius, Jordan, Cavalieri and Brès2015; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2017b ). These studies have used a variety of different definitions of the quantity of interest $\boldsymbol{q}$ , spatial domain $\unicode[STIX]{x1D6FA}$ and inner-product weight $\unicode[STIX]{x1D652}$ . Here, we use the full state vector $\boldsymbol{q}=[\unicode[STIX]{x1D70C},u_{x},u_{r},u_{\unicode[STIX]{x1D703}},p]$ in the domain described previously, and we define the inner-product weight such that the induced norm is equivalent to the compressible energy norm proposed by Chu (Reference Chu1965). The SPOD modes are estimated using the procedure outlined in § 3 using blocks containing $N_{f}=256$ snapshots, leading to $N_{b}=78$ realizations of the jet. A standard Hann window is used to reduce spectral leakage. Frequencies are reported in terms of the Strouhal number $St=fD/U_{j}$ .

Figure 10. Spectral POD modes of the Mach 0.4 turbulent jet: (a) SPOD eigenvalues at $St=0.6$ , normalized by the total energy at that frequency; (b) SPOD eigenvalues as a function of frequency, normalized by the total flow energy; (c–l) real part of the pressure field of the first (c,e,g,i,k) and second (d,f,h,j,l) SPOD modes at the five indicated frequencies.

The SPOD modes are summarized in figure 10. The eigenvalues are depicted in two different ways in figure 10(a,b). Figure 10(a) shows the eigenvalues as a function of the mode number for $St=0.6$ . The eigenvalues are normalized by the total energy at this frequency, so each scaled eigenvalue represents the fraction of the energy at $St=0.6$ described by that mode. The leading mode is substantially more energetic than the suboptimal modes and captures ${\sim}22\,\%$ of the flow energy at this frequency. This is indicative of low-rank dynamics at this frequency (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2017b ). The full spectrum of SPOD eigenvalues is shown in figure 10(b) as a function of the frequency. All eigenvalues have been normalized by the total flow energy integrated over all frequencies, so each one can be interpreted as the fraction of the total flow energy described by that mode. The shading of the curves varies linearly from black to white as the mode number increases from 1 to $N_{b}=78$ . The low-rank behaviour observed at $St=0.6$ , indicated by a large gap between the first and second SPOD eigenvalues, persists in the range $0.3\lesssim St\lesssim 1$ .

The pressure fields of the first two SPOD modes are plotted for several frequencies in figure 10(c–l). All of the modes take the form of coherent wavepacket structures. The wavelength and spatial support of the wavepackets are strongly dependent on the frequency; both quantities increase with decreasing frequency. The wavepacket described by the first mode at each frequency has a single streamwise maximum. These wavepackets have been linked to the Kelvin–Helmholtz instability of the annular jet shear layer for $St\gtrsim 0.3$ (Suzuki & Colonius Reference Suzuki and Colonius2006; Gudmundsson & Colonius Reference Gudmundsson and Colonius2011) and to a superposition of non-normal modes for $St\lesssim 0.3$ (Jordan et al. Reference Jordan, Zhang, Lehnasch and Cavalieri2017; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2017b ). This link between the real turbulent jet and these simple concepts from linear stability theory provides a starting point for the construction of non-empirical reduced-order models. The wavepackets described by the second mode at each frequency have two streamwise maxima, and subsequent modes have increasingly complex structure. These modes optimally describe the variability in the shape of the wavepackets observed in the jet at different times (see § 6.2.3). These variations are central to the acoustic characteristics of the jet (Cavalieri et al. Reference Cavalieri, Jordan, Agarwal and Gervais2011; Cavalieri & Agarwal Reference Cavalieri and Agarwal2014).

Overall, the SPOD modes provide valuable insights that enhance physical understanding of the jet dynamics and motivate reduced-order models. In particular, the wavepackets represented by the SPOD modes have been shown to play a key role in generating acoustic radiation and provide a rigorous starting point for modelling and mitigating jet noise.

6.2.2 Space-only POD modes

Second, we examine the space-only POD modes of the jet. We will see that, in contrast to SPOD, the space-only POD modes provide little insight into the jet dynamics. The space-only POD eigenvalues are shown in figure 11(a). The eigenvalues decay slowly and give no indication of the low-rank behaviour of the jet revealed by SPOD. The power spectral density of each expansion coefficient is shown in figure 11(b). The contour levels are logarithmically distributed over three orders of magnitude, with the upper limit equal to the maximum value. The first approximately 10 modes are dominated by low frequencies that are below the relevant range for jet-noise research; even though these are the highest-energy structures in the jet, they are of little interest. All of the modes, and especially for $j\gtrsim 10$ , contain contributions from a wide range of frequencies – each mode represents behaviour at many different time scales.

Figure 11. Space-only POD modes of the Mach 0.4 turbulent jet: (a) space-only POD eigenvalues, normalized by the total flow energy; (b) PSD of the space-only POD expansion coefficients; (c–j) real part of the pressure field of space-only POD modes 1, 10, 20, 40, 60, 80, 100 and 200.

Because of this mixing of different time scales, the space-only POD modes cannot reproduce the orderly wavepacket structures identified by SPOD. The pressure fields of several modes are shown in figure 11(c–h). Modes 1 and 10 are dominated by large structures associated with the low frequencies that these modes represent. The higher modes have a disorganized appearance caused by the superposition of many different spatial scales associated with the wide range of time scales contained in each mode. Fragments of different wavepackets can be observed in these modes, but a clean separation of distinct structures is not achieved; as indicated by (2.28), the space-only POD modes are made up of a combination of many SPOD modes at different frequencies. This scale mixing makes the space-only POD modes far less helpful for understanding and modelling the jet dynamics. This is perhaps unsurprising in light of the analysis in § 2.3 showing that space-only POD modes do not represent physically meaningful coherent structures.

6.2.3 DMD modes

Third, we compare four different types of DMD modes. In all cases, we use the formulation given by (4.1) with different choices of the matrices $\mathbf{X}$ and $\mathbf{Y}$ . The first set of DMD modes is obtained using all 10 000 snapshots of the jet to define $\mathbf{X}$ and $\mathbf{Y}$ as in (4.2) without subtracting the mean – this is the standard way in which DMD would be applied to the jet database. The second set of DMD modes is obtained again using all 10 000 snapshots in the same way, but the mean is subtracted – this is equivalent to taking a single long-time DFT of the database. The third set of DMD modes is obtained from the ensemble DMD problem defined by (4.4) and the mean is subtracted from each flow realization – the DMD modes are thus the short-time DFT modes of each realization, as shown in § 4. The details of the ensemble definition are the same as described in § 6.2.1 for the SPOD computation. The fourth and final set of DMD modes is the SPOD modes constructed from the ensemble of short-time DFT modes.

Figure 12. Dynamic mode decomposition/DFT modes of the Mach 0.4 turbulent jet: (a) eigenvalues for DMD without mean subtraction; (b) eigenvalues in the vicinity of $St=0.8$ for (●) DMD without mean subtraction; (○, green) DMD with mean subtraction, i.e. a long-time DFT; (▪, blue) ensemble DMD, i.e. short-time DFT; (c–l) real part of the pressure field of the different types of DMD modes; (c,d) DMD without mean subtraction; (e,f) long-time DFT; (g–l) three realizations of the short-time DFT. Panels (c,e,g,i,k) and (d,f,h,j,l) show the modes nearest to $St=0.4$ and 0.8 respectively.

The DMD eigenvalues associated with the different sets of DMD modes are depicted in figure 12(a,b). Figure 12(a) shows the full set of eigenvalues for the standard DMD calculation in which the mean is not removed. With the exception of a small number of outliers, the eigenvalues are tightly clustered along the unit circle and therefore have nearly zero growth/decay rates. This is expected since the jet is a stationary flow. Figure 12(b) focuses on the portion of the unit circle corresponding to $St\approx 0.8$ and confirms that the growth/decay rates of the standard DMD modes are indeed small. Also shown in this plot are the DMD eigenvalues associated with the long-time and short-time DFT (only one short-time eigenvalue is visible). Most of the eigenvalues from the standard DMD calculation lie near one of the eigenvalues from the long-time DFT, which shows that the standard DMD eigenvalues are nearly evenly spaced and suggests that each standard DMD mode is probably similar to one of the long-time DFT modes, even without mean subtraction.

The pressure fields of some of the DMD modes are shown in figure 12(c–l). Specifically, we show the modes from each different type of DMD that are nearest to the frequencies $St=0.4$ and 0.8. The standard DMD and long-time DFT modes are shown in figures 12(c,d) and 12(e,f) respectively. We see that these modes are quite similar, which is expected both from examination of the associated DMD eigenvalues and from the theoretical connection between DMD modes, Koopman modes and Fourier modes discussed in § 4. Three of the short-time DFT modes associated with the ensemble DMD problem are shown for the two frequencies in figure 12(g–l). Significant variability is observed between the different realizations, which is a consequence of the physical flow variability in different short-time intervals and the uncertainty inherent to the DFT. The SPOD modes shown in figure 10 provide an optimal basis for describing the ensemble of short-time DFT modes at each frequency. It is thus not surprising that traces of the wavepackets described by the first two SPOD modes can be observed in the short-time DFT modes and indeed in all of the different types of DMD modes. However, none of these other DMD/DFT modes clearly isolate the coherent wavepacket structures. Instead, they each show one of the many ways in which the coherent structures represented by the SPOD modes can coexist in one particular realization of the turbulent flow. Clearly, the elementary coherent structures provided by SPOD are more useful for understanding and modelling the jet than these random realizations provided by the other forms of DMD.

6.2.4 Resolvent modes

Fourth, we compute resolvent modes of the jet and compare them with the SPOD modes. Several authors have computed resolvent modes for jets in recent years (Garnaud et al. Reference Garnaud, Lesshafft, Schmid and Huerre2013; Jeun et al. Reference Jeun, Nichols and Jovanović2016; Semeraro et al. Reference Semeraro, Jaunet, Jordan, Cavalieri and Lesshafft2016a ,Reference Semeraro, Lesshafft, Jaunet and Jordan b ; Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2017b ), and their properties have been shown to predict certain properties of the jet. Our focus here is on making comparisons with SPOD, as suggested by our analysis in § 5.

The jet is governed by the compressible Navier–Stokes equations, so the flow operator ${\mathcal{A}}$ in (5.3) corresponds to these equations linearized about the mean flow computed from the LES data. The output operator ${\mathcal{C}}$ selects the domain $\unicode[STIX]{x1D6FA}=\{x,r\in [0,\,30]\times [0,\,6]\}$ in order to match the support of the LES snapshots used for the empirical methods. The input operator ${\mathcal{B}}$ multiplies each component of $\boldsymbol{\unicode[STIX]{x1D702}}$ by the turbulent kinetic energy of the jet (again computed from the LES data). This choice is motivated by the analysis of Towne et al. (Reference Towne, Brès and Lele2017) and emphasizes the forcing terms in regions of the jet where turbulent fluctuations are of high enough amplitude to interact nonlinearly. The same compressible energy inner product (Chu Reference Chu1965) as used for the POD modes is adopted for both the input and output spaces.

The equations are discretized using fourth-order finite differences on a grid with $950\times 250$ points in $x$ and $r$ respectively. The computational domain includes the physical domain described above plus a surrounding sponge region that prevents waves from being reflected back into the computational domain at its boundaries. The overall numerical scheme is the same as in Schmidt et al. (Reference Schmidt, Towne, Colonius, Cavalieri, Jordan and Brès2017a ); additional details are available in that reference. As indicated by the above definition of $\unicode[STIX]{x1D6FA}$ , the sponge regions carry zero weight in the discretized inner products. The numerical approximations of the leading resolvent output modes $\tilde{\mathbf{U}}$ and singular values $\unicode[STIX]{x1D6BA}$ are calculated from the eigenvalue decomposition $\tilde{\mathbf{R}}\tilde{\mathbf{R}}^{\ast }=\tilde{\mathbf{U}}\unicode[STIX]{x1D6BA}^{2}\tilde{\mathbf{U}}^{\ast }$ using a standard Arnoldi method (see appendix C for a description of the preceding notation for the discretized equations).

Figure 13. Resolvent modes of the Mach 0.4 turbulent jet: gain as a function of mode number for (a) $St=0.6$ and (b) $St=0.2$ ; (c–l) comparison between the leading SPOD mode (c,e,g,i,k) and the resolvent output mode (d,f,h,j,l) at the five indicated frequencies. The real part of the pressure field of each mode is shown.

The theory developed in § 5 provides a rigorous basis for comparing SPOD and resolvent modes. Figure 13(a,b) shows the gains of the first 100 resolvent modes for $St=0.6$ and 0.2. A large gain separation between the first and second modes is observed for $St=0.6$ but not for $St=0.2$ . This mirrors the behaviour of the SPOD eigenvalues observed in figure 10(b). The leading resolvent modes provide an accurate approximation of the leading SPOD modes at frequencies that exhibit this low-rank behaviour. For example, the pressure fields of the leading SPOD and resolvent output modes are compared for several frequencies in figure 13(c–l). The low-rank behaviour is observed in the resolvent gain and SPOD eigenvalue spectra for all but the lowest frequency shown in the figure, and for these frequencies the resolvent modes approximate the corresponding SPOD modes with striking accuracy. The only discernible difference is that the resolvent-mode wavepackets tend to be slightly longer. On the other hand, the leading resolvent mode at $St=0.2$ is quite different from the corresponding SPOD mode. Similarly, the suboptimal resolvent modes at all frequencies (not shown) do not have a one-to-one correspondence with the suboptimal SPOD modes.

The theoretical connection between SPOD and resolvent modes enables us to glean useful information from these observations. The close match between the leading SPOD and resolvent modes for $St\gtrsim 0.3$ implies that the wavepackets described by these modes are robust and insensitive to correlations within the nonlinear forcing field. Conversely, the mismatch of lower-frequency and suboptimal modes indicates that it is necessary to account for correlated forcing in order to model the low-frequency jet dynamics as well as the wavepacket variability described by the suboptimal SPOD modes at all frequencies. Leveraging these ideas to better understand the jet dynamics and construct accurate reduced-order models is the subject of ongoing research (Schmidt et al. Reference Schmidt, Towne, Rigas, Colonius and Brès2017b ; Towne et al. Reference Towne, Brès and Lele2017).