J. Fluid Mech.

, page 1 of 13 
c Cambridge University Press 2009 1
doi:10.1017/S0022112009992059

Spectral analysis of nonlinear ows

C L A R E N C E W. R O W L E Y1 , I G O R M E Z I C2 ,
S H E R V I N B A G H E R I3 , P H I L I P P S C H L A T T E R 3
A N D D A N S. H E N N I N G S O N3
1
Department of Mechanical and Aerospace Engineering, Princeton University, NJ 08544, USA
2
Department of Mechanical Engineering, University of California, Santa Barbara, CA 93106-5070, USA
3
Linne Flow Centre, Department of Mechanics, Royal Institute of Technology (KTH), SE-10044
Stockholm, Sweden

(Received 15 May 2009; revised 8 September 2009; accepted 9 September 2009)

We present a technique for describing the global behaviour of complex nonlinear

ows by decomposing the ow into modes determined from spectral analysis of
the Koopman operator, an innite-dimensional linear operator associated with the
full nonlinear system. These modes, referred to as Koopman modes, are associated
with a particular observable, and may be determined directly from data (either
numerical or experimental) using a variant of a standard Arnoldi method. They
have an associated temporal frequency and growth rate and may be viewed as a
nonlinear generalization of global eigenmodes of a linearized system. They provide
an alternative to proper orthogonal decomposition, and in the case of periodic data
the Koopman modes reduce to a discrete temporal Fourier transform. The Arnoldi
method used for computations is identical to the dynamic mode decomposition
recently proposed by Schmid & Sesterhenn (Sixty-First Annual Meeting of the APS
Division of Fluid Dynamics, 2008), so dynamic mode decomposition can be thought of
as an algorithm for nding Koopman modes. We illustrate the method on an example
of a jet in crossow, and show that the method captures the dominant frequencies
and elucidates the associated spatial structures.

1. Introduction
Many uid ows exhibit complex phenomena that occur on a wide range of scales
in both space and time. Even with large amounts of information available from
simulations, and comprehensive experimental measurements such as time-resolved
particle image velocimetry (PIV), analysis of complex ow phenomena directly from
raw time histories of the dynamics is usually not fruitful. In practice, one often analyses
ow structures by decomposing them into modes. Common techniques include global
eigenmodes for linearized dynamics (see e.g. Bagheri et al. 2009b), discrete Fourier
transforms, proper orthogonal decomposition (POD) for nonlinear ows (Holmes,
Lumley & Berkooz 1996), balancing modes for linear systems (Rowley 2005) and
many variants of these techniques, such as using shift modes (Noack et al. 2003) in
conjunction with POD modes.
Here, we present a modal decomposition for nonlinear ows based on spectral
analysis of a linear operator, called the Koopman operator, that is dened for any

Email address for correspondence: [email protected]

2 C. W. Rowley, I. Mezic, S. Bagheri, P. Schlatter and D. S. Henningson
nonlinear system. Even if the governing dynamics are nite dimensional, the Koopman
operator is innite dimensional, and does not rely on linearization of the dynamics:
indeed, it captures the full information of the nonlinear system. This operator has been
used to analyse nonlinear dynamical systems, for instance in Mezic & Banaszuk (2004)
and Mezic (2005), and in these works it was shown that for nonlinear systems evolving
on an attractor, modes corresponding to eigenvalues of the Koopman operator may
be computed using harmonic averages or discrete Fourier transforms.
The paper is organized as follows. In 2, we dene the Koopman operator and
its modes associated with a particular observable. In 3, we show that one may
compute approximations to these eigenvalues and their associated modes using a
version of the familiar Arnoldi algorithm that does not require knowledge of an
underlying linear operator. This algorithm is the same as that referred to as dynamic
mode decomposition by Schmid & Sesterhenn (2008). Finally, in 4, we illustrate the
method on an example of a jet in crossow.

2. Koopman operator and Koopman modes

Consider a dynamical system evolving on a manifold M such that, for x k M,
x k+1 = f (x k ), (2.1)
where f is a map from M to itself and k is an integer index. Note that one could
equivalently study continuous-time systems of the form x(t) = f (x(t)), but here we
adopt the discrete-time setting, as we are ultimately interested in analysing discrete-
time data. The Koopman operator is a linear operator U that acts on scalar-valued
functions on M in the following manner: for any scalar-valued function g : M R,
U maps g into a new function Ug given by
Ug(x) = g( f (x)). (2.2)
Note that U is a linear operator, since U (g1 + g2 )(x) = Ug1 (x) + Ug2 (x) for any
functions g1 , g2 and scalars , . Although the dynamical system is nonlinear and
evolves on a nite-dimensional manifold M, the Koopman operator U is linear, but
innite dimensional.
The idea is to analyse the ow dynamics governed by (2.1) only from available
data collected either numerically or experimentally using the eigenfunctions and
eigenvalues of U . To this end, let j : M R denote eigenfunctions and j C
denote eigenvalues of the Koopman operator,
U j (x) = j j (x), j = 1, 2, . . . , (2.3)
and consider a vector-valued observable g : M Rp . For instance, if x M contains
the full information about a ow eld at a particular time, g(x) is a vector of any
quantities of interest, such as a velocity measurements at various points in the ow. If
each of the p components of g lies within the span of the eigenfunctions j , then, as
in Mezic (2005), we may expand the vector-valued g in terms of these eigenfunctions
as

g(x) = j (x)v j . (2.4)
j =1

Here, we have assumed that components of g lie within the span of the eigenfunctions
of U . If this is not the case, then to proceed rigorously, one may split U into
 Spectral analysis of nonlinear ows 3
regular and singular components, and project components of g onto the span of the
eigenfunctions, as done in Mezic (2005).
We typically think of the expression (2.4) as expanding g(x) as a linear combination
of the vectors v j , but we may alternatively think of this expression as expanding the
function g(x) as a linear combination of the eigenfunctions j of U , where now
v j are the (vector) coecients in the expansion. In this paper, we will refer to
the eigenfunctions j as Koopman eigenfunctions, and the corresponding vectors v j
in (2.4) the Koopman modes of the map f , corresponding to the observable g. Note
that iterates of x 0 are then given by



g(x k ) = U k j (x 0 )v j = kj j (x 0 )v j . (2.5)
j =1 j =1

The Koopman eigenvalues, j C, therefore characterize the temporal behaviour of

the corresponding Koopman mode v j : the phase of j determines its frequency, and
the magnitude determines the growth rate. Note that, as described in Mezic (2005),
for a system evolving on an attractor, the Koopman eigenvalues always lie on the
unit circle.
The following examples illustrate that eigenvalues and eigenfunctions of the
Koopman operator are related to objects we routinely use in uid mechanics, such
as global eigenmodes (for linear systems) and the discrete Fourier transform (for
periodic solutions of (2.1)).

2.1. Example: Koopman modes for linear systems

Suppose M is an n-dimensional linear space, and suppose the map f is linear, given
by
f (x) = Ax. (2.6)
It turns out that eigenvalues of A are also eigenvalues of U , and the eigenvectors
of A are related to eigenfunctions of U as well.
Let v j and j denote eigenvectors and eigenvalues of A:
Av j = j v j , j = 1, . . . , n, (2.7)

and let wj be corresponding eigenfunctions of the adjoint A (i.e. A w j = j wj ),

normalized so that v j , wk  = j k , where ,  denotes an inner product on M. Next,
dene scalar-valued functions
j (x) = x, wj  , j = 1, . . . , n. (2.8)
Then j are eigenfunctions of U , with eigenvalues j , since
U j (x) = j ( Ax) =  Ax, wj  = x, A w j  = j x, wj  = j j (x). (2.9)
(Note that, unlike A, the operator U has a countably innite number of eigenvalues,
since kj is also an eigenvalue, with eigenfunction j (x)k , for any integer k.) Now, for
any x M, as long as A has a full set of eigenvectors, we may write

n 
n
x= x, wj  v j = j (x)v j . (2.10)
j =1 j =1

For linear systems, then, the Koopman modes coincide with the eigenvectors of A.
4 C. W. Rowley, I. Mezic, S. Bagheri, P. Schlatter and D. S. Henningson
2.2. Example: Koopman modes for periodic solutions
Returning to the nonlinear setting, suppose we have a set of distinct vectors
S = {x 0 , . . . , x m1 } that form a periodic solution of (2.1), such that x k+m = x k for
all k. A common way to analyse such a solution is to take its discrete Fourier
transform, which denes a new set of vectors { x 0 , . . . , x m1 } that satisfy

m1
xk = e2ij k/m x j , k = 0, . . . , m 1. (2.11)
j =0

Now, dene a set of functions j : S C by

j (x k ) = e2ij k/m , j, k = 0, . . . , m 1. (2.12)
Then j are eigenfunctions of the Koopman operator U , with eigenvalues e2ij/m ,
since
U j (x k ) = j ( f (x k )) = j (x k+1 ) = e2ij (k+1)/m = e2ij/m j (x k ). (2.13)
Therefore, we may write the expansion (2.11) as

m1
xk = j (x k ) x j . (2.14)
j =0

Note the similarity with (2.10). Thus, if we restrict our phase space to the periodic
orbit S, the Koopman modes dened as in the previous subsection are the vectors
x j given by the discrete Fourier transform, and the phases of the corresponding
eigenvalues are the frequencies 2j/m. This result in fact applies more generally to
non-periodic systems, as discussed in Mezic & Banaszuk (2004) and Mezic (2005):
when the dynamics are restricted to any attractor, the Koopman modes may be
calculated by harmonic averages, which for nite-time datasets reduce to discrete
Fourier transforms.

3. Computation of Koopman modes from snapshots

Here, we present an algorithm for computing the Koopman modes dened in the
previous section, given only values of a particular observable (snapshots), sampled
at regular times. As before, we will assume that the dynamics are governed by (2.1),
and that for any state x, we can measure a vector-valued observable g(x) Rp . For
instance, if we have access to a grid of velocity vectors in a limited region of space (e.g.
obtained via PIV), then p is the number of grid points times velocity components.
In particular, below we show that the commonly used Arnoldi algorithm, when
applied to a nonlinear system, actually produces approximations to eigenvalues of the
Koopman operator, and their corresponding modes as dened in the previous section.
We rst consider linear systems, and present a version of the Arnoldi algorithm that
does not require explicit knowledge of the underlying operator A. This variant of the
algorithm is described by Saad (Saad 1980, p. 287), and is the same as that referred
to as dynamic mode decomposition by Schmid & Sesterhenn (2008). We then provide
an alternative interpretation of the algorithm that applies to nonlinear systems, and
connects with the Koopman modes.
3.1. Arnoldi algorithm for linear systems
Assume one has a linear dynamical system
x k+1 = Ax k , (3.1)
 Spectral analysis of nonlinear ows 5
where x k R , and n is so large that we cannot compute eigenvalues of A directly.
n

A standard method for computing estimates of these eigenvalues is a Krylov method,

in which one starts with an initial vector x 0 (often chosen to be a random vector),
and computes iterates of x 0 . After m 1 iterations, one has a collection of m vectors
that span a Krylov subspace, given by span{x 0 , Ax 0 , . . . , Am1 x 0 }. One then nds
approximate eigenvalues and eigenvectors by projecting A onto this m-dimensional
subspace, and computing eigenvectors and eigenvalues of the resulting low-rank
operator. If we stack the data vectors into an n m matrix

K = [x 0 Ax 0 A2 x 0 Am1 x 0 ] (3.2)
= [x 0 x 1 x 2 x m1 ], (3.3)

then we wish to nd approximate eigenvectors of A as linear combinations of

the columns of K . The Arnoldi algorithm is a type of Krylov method in which
one orthonormalizes the iterates at each step, and it therefore involves computing the
action of A on arbitrary vectors. A variant of this algorithm that does not require
explicit knowledge of A is given below.
First, consider the special case where the mth iterate x m is a linear combination of
the previous iterates. Following 2 of Ruhe (1984), we may write

x m = Ax m1 = c0 x 0 + + cm1 x m1 = K c, (3.4)

where c = (c0 , . . . , cm1 ). Thus, we have

AK = K C, (3.5)

where C is a companion matrix given by

0 0 0 c0
1 0 0
c1

C=
0 1 0 c2 .
(3.6)
.. .. ..
. . .
0 0 1 cm1

The eigenvalues of C are then a subset of the eigenvalues of A: if

Ca = a,

then using (3.5), one veries that v = K a is an eigenvector of A, with eigenvalue .

More generally, if the mth iterate is not a linear combination of the previous
iterates, then instead of the equality (3.4), we have a residual

r = x m K c,

which is minimized when c is chosen such that r is orthogonal to span{x 0 , . . . , x m1 }.

In this case, the relation (3.5) becomes AK = K C + r eT , where e = (0, . . . , 1) Rm .
The eigenvalues of C are then approximations to the eigenvalues of A, called Ritz
values, and the corresponding approximate eigenvectors are given by v = K a, called
Ritz vectors. Note that the full Arnoldi method is more numerically stable than this
method, and reduces A to an upper Hessenberg matrix, rather than a companion
matrix.
6 C. W. Rowley, I. Mezic, S. Bagheri, P. Schlatter and D. S. Henningson
3.1.1. Algorithm
An important feature of the above algorithm is that it does not require explicit
knowledge of the matrix A: all it requires is a sequence of vectors, as summarized
below.
Consider a sequence {x 0 , . . . , x m } where x j Rn . Dene the empirical Ritz values j
and empirical Ritz vectors v j of this sequence by the following algorithm:
(i) Dene K by (3.3) and nd constants cj such that

m1
r = xm cj x j = x m K c, r span{x 0 , . . . , x m1 }. (3.7)
j =0

(ii) Dene the companion matrix C by (3.6) and nd its eigenvalues and
eigenvectors
C = T1 T, = diag(1 , . . . , m ), (3.8)
1
where eigenvectors are columns of T .
(iii) Dene v j to be the columns of V = K T1 .
If x j = Aj x 0 , then the empirical Ritz values j are the usual Ritz values of A after
m steps of the Arnoldi method, and v j are the corresponding Ritz vectors. These,
then, are (usually) good approximations of the eigenvalues and eigenvectors of A.
However, if we do not have x j = Aj x 0 (for instance, if the sequence is generated by a
nonlinear map), then at this point, it is not clear what the above algorithm produces.
We show below that for a nonlinear system, the algorithm produces approximations
of the Koopman modes and associated eigenvalues.
3.2. Modal decomposition for nonlinear systems
A more general interpretation of the above algorithm is provided by the following
theorem, which will be used below.
Theorem 1. Consider a set of data {x 0 , . . . , x m }, and let j , v j be the empirical Ritz
values and vectors of this sequence. Assume the j are distinct. Then

m
xk = kj v j , k = 0, . . . , m 1, (3.9)
j =1

m
xm = j v j + r,
m r span{x 0 , . . . , x m1 }. (3.10)
j =1

Proof. Note that (3.9) may be written in matrix form as

1 1 21 1m1
1 2 22 2
m1
K := [x 0 x 1 x m1 ] = [v 1 v m ]
... .. .. .. . (3.11)
.
..
. . .
1 m 2m m1
m
The rightmost matrix above is a Vandermonde matrix, which we will denote T. Note
that Vandermonde matrices and companion matrices are closely related, in that T
diagonalizes the companion matrix C dened in (3.6), as long as the eigenvalues
1 , . . . , m are distinct. That is, T is precisely the matrix T in (3.8), and so K = V T,
which veries (3.11), and therefore (3.9). Equation (3.10) then follows from the last
 Spectral analysis of nonlinear ows 7
column of the equality
[x 1 x2 x m ] = K C + r eT = K T1 T + r eT = V T + r eT . (3.12)

To illustrate how this theorem provides a connection with Koopman modes, consider
a vector-valued observable g : M Rp for the dynamical system (2.1) and its
expansion (2.4) in the Koopman modes. Suppose we have a sequence of observations
g(x k ), for k = 0, . . . , m, and let j and v j be the empirical Ritz values and vectors for
this sequence. Then, by Theorem 1, we have

m
g(x k ) = kj v j , k = 0, . . . , m 1, (3.13)
j =1

m
g(x m ) = j v j + r,
m (3.14)
j =1

with r span{g(x 0 ), . . . , g(x m1 )}. Comparing with the expansion (2.5), the empirical
Ritz values j and vectors v j behave in precisely the same manner as the eigenvalues
j and modes v j of U , but for the nite sum in (3.13) instead of the innite sum (2.5).
If r = 0 in (3.14), then as far as the data is concerned, the approximate modes
are indistinguishable from true eigenvalues and Koopman modes of U , with the
expansion (2.4) consisting only of a nite number of terms.
If r = 0, then there is some error, but this is in a sense the best one can do, since
the m + 1 observations cannot in general be spanned by m modes. In fact, by the
projection theorem, the error r in (3.14) is the same as the smallest possible error in
projecting g(x m ) onto any modes v j formed from linear combinations of the rst m
data vectors. In this sense, the values j are then approximations of true eigenvalues j
of U , and the vectors v j are approximations of the Koopman modes v j , scaled by the
constant values j (x 0 ).
Note that an interesting situation occurs if g(x 0 ) = g(x m ). Then in (3.7), c0 = 1 and
cj = 0 for all other values of j . The empirical Ritz values are then the mth roots of
unity, j = e2ij/m , and the Vandermonde matrix T is the discrete Fourier transform
matrix. Thus, in this case, the empirical Ritz vectors are given by the discrete Fourier
transform of the data, at frequencies 2j/m, as illustrated by the example in 2.2.

4. Example: jet in crossow

The jet-in-crossow conguration appears in a variety of applications and is a
common way of mixing a jet uid injected through an orice with a uniform
crossow. Recently, Bagheri et al. (2009b) showed that the jet in crossow exhibits
self-sustained global oscillations that can be associated with vortex shedding in
dierent spatial regions. Using time traces, we extract and quantify here the oscillatory
behaviour of the ow from fully nonlinear direct numerical simulations (DNS) and
show that the computed Koopman modes identify the relevant frequencies and the
corresponding three-dimensional ow structures automatically.
The parameters and the numerical code are the same as in the DNS performed in
Bagheri et al. (2009b); the jet inow ratio is R Vj et /U = 3, the Reynolds number
is Re0 U 0 / = 165 and the ratio between the crossow displacement thickness
and the jet diameter is 0 /D = 1/3. The incompressible NavierStokes equations over
a at plate are solved using a FourierChebyshev spectral method (Chevalier et al.
2007) and the jet with an initially parabolic velocity prole is introduced as an
8 C. W. Rowley, I. Mezic, S. Bagheri, P. Schlatter and D. S. Henningson

(a) (b)

(c) 0 (d)
1.0
0.5
0.95
u1 u2
1.0 0.9
1.5
200 250 300 350 400 450 500 550 600 650 200 250 300 350 400 450 500 550 600 650
t t
(e) 1.0 ( f) 1.0

0.8 0.8

0.6 0.6

0.4 0.4

0.2 0.2

0 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225 0.250 0 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 0.225 0.250
St St

Figure 1. (a) Snapshot of the ow eld at t = 400. Red and grey isocontours represent
2 = 0.1 and u = 0.2 (near the wall) respectively. (b) The same quantities for the time-averaged
ow which also is the rst Koopman mode. (c, d ) Time signal probes located near the wall and
on the jet shear layer respectively (see the text). (e, f ) The spectral content of the corresponding
time signals in (c) and (d ) are shown in black and the magnitudes of the rst seven pairs of
Koopman modes at each frequency are shown in red. The amplitudes are normalized with
their maximum values.

inhomogeneous boundary condition for the wall-normal velocity component at the

wall (y = 0). The grid resolution is 256 201 144 grid points in a computational
box (Lx , Ly , Lz ) = (75, 20, 30)0 . The three-dimensional ow behaviour is triggered by
an asymmetric localized perturbation imposed at t = 0. For the exact form of the
jet-prole and further numerical details see Bagheri et al. (2009b).
The ow physics of the jet-in-crossow has been studied extensively (see e.g. Fric &
Roshko 1994; Kelso, Lim & Perry 1996; Muppidi & Mahesh 2007) and it is shown
that it is mainly characterized by four to ve vortical structures depending on R and
Re. In the present study, we could identify two steady, two unsteady and one quasi-
steady vortex systems. Four of these can be seen in gure 1(a, b) where isocontours
of the 2 criterion (red) (Jeong & Hussain 1995) and the streamwise velocity (grey)
are displayed. The most signicant unsteady feature of the jet is the highly unsteady
shear-layer vortices (see gure 1a): These half-ringed shaped vortices grow along the
jet trajectory, lead to a breakdown of the ordered ow and, eventually, dissipate due
to viscous eects. Connected to the shear layer are the vertically oriented upright
vortices, which are periodically appearing vortices connecting the jet body and the
wall layer in the wake of the jet. These structures are easily identied from the vorticity
eld, but are not visible here. On the other hand, the steady counter-rotating vortex
 Spectral analysis of nonlinear ows 9
(a) (b) 400
1.0

300
0.5

|| j||
Im{j}

200
0

100
0.5

1.0 0 0.03 0.06 0.09 0.12 0.15 0.18 0.21 0.24

1.0 0.5 0 0.5 1.0 St
Re{j}

Figure 2. (a) The empirical Ritz values j . The value corresponding to the rst Koopman
mode is shown with the blue symbol. (b) The magnitudes of the Koopman modes (except the
rst one) at each frequency. In both gures, the colours vary smoothly from red to white,
depending on the magnitude of the corresponding mode.

pair (CVP), characteristic of the jet trajectory, is visible in the time-averaged mean
ow (gure 1b). Similarly visible in the mean, the horse-shoe vortices wrap around the
column of jet very close to the wall and, further downstream, lead to the appearance
of the quasi-steady wall-vortex system as shown by the distortion of the velocity
isocontour in gure 1(a, b). These essentially streamwise-oriented vortices are subject
to low-frequency oscillations of the wall-vortex system originating in a shedding of
the separation zone just downstream of the jet orice, and induce a movement of
the whole jet body. The steady structures were also identied in the steady nonlinear
NavierStokes solution computed numerically by Bagheri et al. (2009b).
As found in Bagheri et al. (2009b), two distinct self-sustained oscillations could
be detected from the DNS. A high-frequency shedding of the shear-layer vortices
and a very low frequency shedding of the wall vortices. Figure 1(c) shows the time
signal of the streamwise velocity u1 (x 1P , t) from a probe located just downstream
of the jet orice and close the wall, x 1P = (x, y, z) = (10.7, 1, 0). In gure 1(e) its
corresponding power spectrum shows the frequency content u1 () of u1 (t). The peak
frequency corresponds to a vortex shedding of wake vortices with the Strouhal number
St f D/Vj et = 0.0174.
In gure 1(d, f ), a second probe located a few jet diameters along the jet trajectory
x 2P = (12, 6, 2), shows a second oscillation that can be identied with the shedding of
the shear-layer vortices. The peak frequency oscillates with St = 0.141 which is nearly
one order of magnitude larger than the low-frequency mode. Note that the peak
frequencies of the power spectra vary slightly depending on the location of the probe.
4.1. Koopman modes and frequencies
In this section we compute the Koopman modes and show that they directly allow
an identication of the various shedding frequencies. The empirical Ritz values
j and the empirical vectors v j of a sequence of ow elds {u0 , u1 , . . . , um1 } =
{u(t = 200), u(t = 202), . . . , u(t = 700)} with m = 251 are computed using the
algorithm described earlier. Thus, the transient time (t < 200) is not sampled and
only the asymptotic motion in phase space is considered.
Figure 2(a) shows that nearly all the Ritz values are on the unit circle |j | = 1
indicating that the sample points ui lie on or near an attracting set. The Koopman
eigenvalue corresponding to the rst Koopman mode is the time-averaged ow and
10 C. W. Rowley, I. Mezic, S. Bagheri, P. Schlatter and D. S. Henningson
(a) (b)

Figure 3. Positive (red) and negative (blue) contour levels of the streamwise velocity
components of two Koopman modes. The wall is shown in grey. (a) Mode 2, with v 2 = 400
and St2 = 0.141. (b) Mode 6, with v 6 = 218 and St6 = 0.0175.

is depicted with blue symbol in gure 2(a). This mode, shown in gure 1(b), captures
the steady ow structures as discussed previously. In gure 2(a), the other (unsteady)
Ritz values vary smoothly in colour from red to white, depending on the magnitude
of the corresponding Koopman mode. The magnitudes dened by the global energy
norm v j , and are shown in gure 2(b) with the same colouring as the spectrum. In
gure 2(b) each mode is displayed with a vertical line scaled with its magnitude at
its corresponding frequency j = Im{log(j )}/t (with t = 2 in our case). Only the
j > 0 are shown, since the eigenvalues come in complex conjugate pairs. Ordering
the modes with respect to their magnitude, the rst (23) and second (45) pair of
modes oscillate with St2 = 0.141 and St4 = 0.136 respectively, whereas the third pair
of modes (67) oscillate with St6 = 0.017. All linear combinations of the frequencies
excite higher modes, for instance, the nonlinear interaction of the rst and third pair
results in the fourth pair, i.e. St8 = 0.157 and so on.
In gures 1(e) and 1(f ) the power spectra of the two DNS time signals (black lines)
are compared to the frequencies obtained directly from the Ritz eigenvalues (red
vertical lines). The shedding frequencies and a number of higher harmonics are in
very good agreement with the frequencies of the Koopman modes. In particular, the
dominant Koopman eigenvalues match the frequencies for the wall mode (St = 0.017)
and the shear-layer mode (St = 0.14). Note that the probe signals are local measures
of the frequencies at one spatial point, whereas the Koopman eigenvalues correspond
to global modes in the ow with time-periodic motion.
The streamwise velocity component u of Koopman modes 2 and 6 are shown
in gure 3. Each mode represents a ow structure that oscillates with one single
frequency, and the superposition of several of these modes results in the quasi-
periodic global system. The high-frequency mode 2 (gure 3a) can be associated with
the shear layer vortices; along the jet trajectory there is rst a formation of ring-like
vortices that eventually dissolve into smaller scales due to viscous dissipation. Also
visible are upright vortices: on the leeward side of the jet, there is a signicant
structure extending towards the wall. This indicates that the shear-layer vortices and
the upright vortices are coupled and oscillate with the same frequency. The spatial
structures of modes 4 and 8 are very similar to those of mode 2, as one expects, since
the frequencies are very close.
On the other hand, the low-frequency mode 6 shown in gure 3(b) features large-
scale positive and negative streamwise velocity near the wall, which can be associated
with shedding of the wall vortices. However, this mode also has structures along
the jet trajectory further away from the wall. This indicates that the shedding of
wall vortices is coupled to the jet body, i.e. the low frequency can be detected nearly
anywhere in the vicinity of the jet since the whole jet is oscillating with that frequency.
 Spectral analysis of nonlinear ows 11

Mode DNS Global POD Koopman

Shear layer 0.141 0.169 0.138, 0.158, 0.121 0.141
Wall 0.017 0.043 0.0188, 0.0094, 0.158, 0.121 0.017
Table 1. Comparison of the frequencies (St = f D/Vjet ) obtained from DNS probes (shown
in gure 1); the global eigenmodes of the linearized NavierStokes; POD modes 1 and 6,
corresponding to mainly shear-layer and wall oscillations, respectively; and Koopman modes.

1.0
0.5
0
0.5
1.0

200 250 300 350 400

Figure 4. Comparison of time coecients: the projection of the ow eld onto the most
energetic POD mode (black), and the coecient of the most energetic Koopman mode (grey).

4.2. Comparison with linear global modes and POD modes

The linear global eigenmodes of the NavierStokes equations linearized about an
unstable steady state solution were computed by Bagheri et al. (2009b) for the
same ow parameters as the current study. They computed 22 complexconjugate
unstable modes using the Arnoldi method combined with a time-stepper approach.
The frequency of the most unstable (antisymmetric) mode associated with the shear-
layer instability was St = 0.169, not far from the value St = 0.14 observed for the
DNS. However, the mode with the lowest frequency associated with the wall vortices
was St = 0.043, far from the observed frequency of St = 0.017. These frequencies
are summarized in table 1. The global eigenmodes capture the dynamics only in
a neighbourhood of the unstable xed point, while the Koopman modes correctly
capture the behaviour on the attractor.
We also compared the Koopman modes with modes determined by POD of the same
dataset. The POD modes themselves are shown in Bagheri, Schlatter & Henningson
(2009a) and capture similar spatial structures to the Koopman modes shown in
gure 3. The most striking distinction is in the time coecients, an example of which
is shown in gure 4: while a single Koopman mode contains, by construction, only a
single frequency component, the POD modes capture the most energetic structures,
resulting in modes that contains several frequencies. The coecient of the rst POD
mode oscillates mainly with frequency St = 0.138, which is close to the shear-layer
oscillation frequency St = 0.141 observed in DNS. However, the signal contains other
frequencies as well, resulting from the interaction of the two fundamental oscillations
(shear-layer and wall), St = 0.138 0.017, which cause the beating shown in gure 4.
The frequencies present in this most energetic POD mode are also summarized in
table 1. Higher POD modes (in this case the sixth) capture the wall oscillations,
but the signal is polluted with other frequencies as well. For situations such as the
jet in crossow where one is interested in studying the dynamics of low-frequency
oscillations (such as wall modes) separate from high-frequency oscillations (such as
12 C. W. Rowley, I. Mezic, S. Bagheri, P. Schlatter and D. S. Henningson
shear-layer modes), the Koopman modes are thus more eective at decoupling and
isolating these dynamics.

5. Conclusions
We have presented a method for studying the dynamical behaviour of nonlinear
systems, and illustrated the method on a jet in crossow. The method involves spectral
analysis of the Koopman operator, an innite-dimensional linear operator dened for
any nonlinear dynamical system. In particular, given a particular observable (e.g.
available measurements from an experiment or simulation), we have dened a set
of Koopman modes associated with this observable, related to eigenfunctions of the
Koopman operator. For the special case of linear systems where the observable is the
full ow state, these modes reduce to the global eigenmodes, and for periodic systems,
the modes can be determined by the discrete temporal Fourier transform. For more
general systems, we have shown that these modes may be computed using a familiar
Arnoldi algorithm, applied to samples of the nonlinear system.
We have used these modes to study the dynamical behaviour of a jet in crossow.
The resulting modes illustrate the dierent spatial structures associated with the
shear layer and the near-wall region. For this example, the Koopman modes capture
the relevant frequencies more accurately than global eigenmodes of the linearized
dynamics, and decouple the dierent frequency components more eectively than
modes determined by POD.

The authors gratefully acknowledge support for this work from the National Science
Foundation (CMS-0347239) and the Air Force Oce of Scientic Research (FA9550-
09-1-0257), and computer-time allocation from the Swedish National Infrastructure
for Computing (SNIC).

REFERENCES
Bagheri, S., Schlatter, P. & Henningson, D. S. 2009a Self-sustained global oscillations in a jet in
crossow. Theor. Comput. Fluid Dyn. Manuscript submitted for publication.
Bagheri, S., Schlatter, P., Schmid, P. J. & Henningson, D. S. 2009b Global stability of a jet in
crossow. J. Fluid Mech. 624, 3344.
Chevalier, M., Schlatter, P., Lundbladh, A. & Henningson, D. S. 2007 A pseudo spectral solver
for incompressible boundary layer ows. Tech. Rep., Trita-Mek 7.
Fric, T. F. & Roshko, A. 1994 Vortical structure in the wake of a transverse jet. J. Fluid Mech. 279,
147.
Holmes, P., Lumley, J. L. & Berkooz, G. 1996 Turbulence, Coherent Structures, Dynamical Systems
and Symmetry. Cambridge University Press.
Jeong, J. & Hussain, F. 1995 On the identication of a vortex. J. Fluid Mech. 285, 6994.
Kelso, R. M., Lim, T. T. & Perry, A. E. 1996 An experimental study of round jets in crossow.
J. Fluid Mech. 306, 111144.
Mezic, I. 2005 Spectral properties of dynamical systems, model reduction and decompositions.
Nonlinear Dyn. 41, 309325.
Mezic, I. & Banaszuk, A. 2004 Comparison of systems with complex behaviour. Physica D 197,
101133.
Muppidi, S. & Mahesh, K. 2007 Direct numerical simulation of round turbulent jets in crossow.
J. Fluid Mech. 574, 5984.
Noack, B. R., Afanasiev, K., Morzynski, M., Tadmor, G. & Thiele, F. 2003 A hierarchy of
low-dimensional models for the transient and post-transient cylinder wake. J. Fluid Mech.
497, 335363.
 Spectral analysis of nonlinear ows 13
Rowley, C. W. 2005 Model reduction for uids using balanced proper orthogonal decomposition.
Intl J. Bifurcation Chaos 15 (3), 9971013.
Ruhe, A. 1984 Rational Krylov sequence methods for eigenvalue computation. Linear Algebr. Appl.
58, 391405.
Saad, Y. 1980 Variations on Arnoldis method for computing eigenelements of large unsymmetric
matrices. Linear Algebr. Appl. 34, 269295.
Schmid, P. & Sesterhenn, J. 2008 Dynamic mode decomposition of numerical and experimental
data. In Sixty-First Annual Meeting of the APS Division of Fluid Dynamics. San Antonio,
Texas, USA.