# A minimization principle for the description of time-dependent modes associated with transient instabilities

###### Abstract

We introduce a minimization formulation for the determination of a finite-dimensional, time-dependent, orthonormal basis that captures directions of the phase space associated with transient instabilities. While these instabilities have finite lifetime they can play a crucial role by either altering the system dynamics through the activation of other instabilities, or by creating sudden nonlinear energy transfers that lead to extreme responses. However, their essentially transient character makes their description a particularly challenging task. We develop a minimization framework that focuses on the optimal approximation of the system dynamics in the neighborhood of the system state. This minimization formulation results in differential equations that evolve a time-dependent basis so that it optimally approximates the most unstable directions. We demonstrate the capability of the method for two families of problems: i) linear systems including the advection-diffusion operator in a strongly non-normal regime as well as the Orr-Sommerfeld/Squire operator, and ii) nonlinear problems including a low-dimensional system with transient instabilities and the vertical jet in crossflow. We demonstrate that the time-dependent subspace captures the strongly transient non-normal energy growth (in the short time regime), while for longer times the modes capture the expected asymptotic behavior.

## 1 Introduction

A broad range of complex systems in nature and technology are characterized by the presence of strongly transient dynamical features associated with finite-time instabilities. Examples include turbulent flows in engineering systems (e.g. Kolmogorov [15] and unstable plane Couette flow [26], reactive flows in combustion [52, 41]), turbulent flows in geophysical systems (e.g. climate dynamics [31, 30], cloud process in tropical atmospheric convection [24, 23]), nonlinear waves (e.g. in optics [3, 4] or water waves [38, 37, 20, 19]), and mechanical systems [34, 5, 29, 35, 25].

These systems are characterized by very high dimensional attractors, intense nonlinear energy transfers between modes and broad spectra. Despite their complexity, the transient features of these dynamical systems are often associated with low-dimensional structures, i.e. a small number of modes, whose strongly time-dependent character, however, makes it particularly challenging to describe with the classical notion of time-independent modes. This is because these modes, despite their connection with intense energy transfers and transient dynamics, often have low energy and hence they are “buried” in the complex background of modes that are not associated with intense growth or decay but only with important energy. These transient modes often act as “triggers” or “precursors” of higher energy phenomena or instabilities and a thorough analysis of their properties can have important impact for i) the understanding of the system dynamics and in particular the mechanisms associated with transient features (see e.g. [50, 56, 54]), ii) the prediction and quantification of upcoming instabilities that are triggered through these low-energy dynamical processes (see e.g. [53, 18, 19]), iii) the control and suppression of these instabilities by suitably focusing the control efforts in the low energy regime of these transient phenomena (see e.g. [40, 55, 17]).

Transient dynamics is central in understanding a wide range of fluid mechanics problems. In the context of hydrodynamic stability, non-normality of the linearized Navier-Stokes operator can cause significant *transient energy growth* [10, 14]. It has been well-established now that the eigenvalue analysis fails to predict the short-time evolution of perturbations for convective flows [59]. Instead,
the transient energy growth can be better understood by analyzing the pseudospectra of the linearized operator [58, 42]. Transition from laminar flow to turbulence, is another
active area of research in fluid mechanics, to which, understanding the transient
dynamical features is crucial. For a finite-amplitude disturbance, *bypass transition* has been observed in wall-bounded shear flows [49, 43],
in which case the non-normal growth of localized disturbances leads to small turbulent spots, bypassing the secondary instability process [27]. Recent computational and experimental studies also demonstrate the *sudden transition* from laminar to turbulent motion in pipe flows, where turbulence forms from localized patches called *puffs* [36, 6]. On the other hand,
intermittent behavior is the hallmark of turbulent fluid flows.
Turbulent
*chaotic bursts* appearing in spatially and temporally
localized events can dramatically change the dynamics of the system [21]. Prototype systems that mimic these properties were introduced and analyzed in [32, 13, 33, 34].
Transient dynamics have a fundamental role in the intermittent behavior of passive tracers as well, where even elementary models without positive Lyapunov exponents [11, 12] have been able to reproduce intermittent behavior observed in complex models. For such systems, the fundamental role of the random resonance between Fourier modes of the turbulent velocity field and the passive tracer has been recently illustrated [57].

In the context of uncertainty quantification and prediction a new family of stochastic methods relying on the so-called Dynamical Orthogonality condition was recently developed to deal with the strongly transient features of stochastic systems. The Dynamically Orthogonal (DO) Field Equations [46, 47] and Dynamically Bi-orthogonal (BO) equations [16] evolve a subspace according to the system stochastic PDE and the current statistical state of the system. Despite their success in resolving low-dimensional stochastic attractors for PDES [44, 45, 48], these methodologies are often too expensive to implement for high dimensional systems (e.g. DO or BO require the simultaneous solution of many PDEs). In addition, in order to obtain an accurate description of the time-dependent dynamics many modes should be included in the analysis and the computational cost increases very rapidly (especially for systems with high complexity).

Our aim in this work is to develop a method that will generate adaptively a time-dependent basis that will capture strongly transient phenomena. This approach will rely on system observables obtained either through high fidelity numerical solvers or measurements, as well as the linearized equations of the system. The core of our approach is a minimization principle that will seek to minimize the distance between the local vector field of the system, constrained over the direction of the time-dependent modes, and the rate of change of the time-dependent modes. A direct minimization of the defined functional will result in evolution equations for the optimally time-dependent (OTD) basis elements. For systems characterized by transient responses these modes will adapt according to the (independently) computed or measured system history in a continuous way, capturing at each time the transiently most unstable directions of the system. For sufficiently long times where the system reaches an equilibrium we prove that the developed equations provide the most unstable directions of the system in the asymptotic limit.

We demonstrate the developed approach over a series of applications, including linear and nonlinear systems. As a first example we consider the advection-diffusion operator where we show how the OTD basis captures the directions associated with the non-normal behavior. The second example involves the Orr-Sommerfeld/Squire operator that governs the evolution of infinitesimal disturbances in parallel viscous flows. Our goal here is also the computation of time-dependent modes that explain the transient growth of energy due to non-normal dynamics. The third problem involves low-dimensional dynamical system as well as the nonlinear transient dynamics of a jet in cross flow. Using the developed framework we compute the modes associated with the transient but also asymptotically unstable directions of the phase space and we assess their time-dependent stability properties.

## 2 Optimally time-dependent (OTD) basis for transient instabilities

Let the dynamical system

defined on a state space We denote by the position of the trajectory at time that is initiated on . Also, let

(1) |

denote the linearized dynamical system around the trajectory and let the inner product between two elements and be denoted as . The linear time-dependent dynamical system represented by equation (1) has the solution

(2) |

where is the propagator that maps the state of the system at time to . The propagator can be represented as the ordered product of infinitesimal propagators

(3) |

where lies in with . Our aim is to evolve a basis , i.e. a set of time-dependent, orthonormal modes, so that optimally follows for all times. To achieve this goal we formulate the following quantity, which measures the distance between the action of infinitesimal propagator on an orthonormal basis and . We have:

(4) |

where is an arbitrary and time-dependent orthonormal basis, i.e. for every time instant it satisfies the orthonormality condition:

(5) |

We observe that in the functional given by equation (4):

where is the identity matrix. Replacing the above equations into the functional given by equation (4) results in:

(6) |

In Figure 1 we illustrate the distance function . For each direction we have, due to the normalization property , the rate of change lying in the orthogonal complement of . Under this constraint we choose as the vector that is closest to the image of under the effect of the operate .

We emphasize that the minimization of the function (6) is considered only with respect to the time-derivative (rate of change) of the basis, instead of the basis itself. This is because we do not want to optimize the subspace that the operator is acting on, but rather find an optimal set of vectors, that best approximates the linearized dynamics in the subspace . We then solve the resulted equations and compute . We will refer to these modes as the optimally time-dependent (OTD) modes, and the space that these modes span as the OTD subspace.

Note that a direct minimization of the function (6), over finite-time intervals, with respect to the modes would result the well-known Euler-Lagrange equations. However, in this case the emphasis is put on finding the optimal ‘input’ to the operator in order for the function (6) to be minimized, while here our aim is different, i.e. we seek to find an optimal set of vectors that approximate in the best possible way, when is acting on .

To summarize the above discussion, we are considering the following minimization problem:

In the reminder of this paper we will demonstrate that this minimization principle allows to capture any arbitrary and transient kind of growth caused by the linear operator, e.g. non-normal or exponential. In what follows we obtain a differential formulation for the evolution equations that correspond to the proposed minimization problem.

### 2.1 Evolution equations for the time-dependent modes

Before we proceed to the minimization of the function given by equation (6), we express the orthonormality constraints for the basis elements in terms of their time derivatives. By differentiating over time, we have

(7) |

The above condition is satisfied if the following condition is valid:

(8) |

where is any arbitrary function for which . Clearly, for any choice of the above condition guarantees that if is initially orthonormal, it will remain orthonormal for all times. As we will see the choice of does not change the evolved subspace. However, it allows for different formulations of the evolution equations. Using the constraint (8) we have the following Theorem.

###### Theorem 2.1

Proof: We first formulate the minimization problem that also takes into account the appropriate number of Lagrange multipliers, with . In this way we obtain:

(10) |

We consider the first derivative with respect to each to obtain the set of equations

To obtain an extremum we need the right hand side of the last equation to vanish:

(11) |

which should be solved together with condition (8). We take the inner product of equation (11) with mode and obtain

Using the last equation and substituting in (11) will result in the evolution equations (9). This completes the proof.

In a more compact form the evolution equation for a finite-dimensional operator can be obtained, where we express the OTD subspace in a matrix whose column is . The function is correspondingly expressed in the matrix notation as with . Thus, the evolution equation for the OTD modes can be expressed as:

(12) |

where denotes the transpose of a matrix.

Now we define the *reduced operator*
that is obtained by projecting the original
operator onto the subspace . Thus,

(13) |

Therefore the OTD equation could be equivalently expressed as:

(14) |

In what follows we will need to the definition of equivalence between two subspaces.

###### Definition 2.1

The two OTD subspaces and are equivalent at time if there exists a transformation matrix such that , where is an orthogonal rotation matrix, i.e. .

In the following we show that the evolution of two OTD subspaces and under different choice of , which are initially equivalent (at ), will remain equivalent for every time , *i.e.* where is an orthogonal rotation matrix governed by the matrix differential equation:

(15) | ||||

with and being the two different choices of for the evolution of and , respectively, while in the initial orthogonal rotation matrix, *i.e.* .
We first prove the following Lemma.

###### Lemma 2.1

The solution to the matrix differential equation given by equation (15), remains an orthogonal rotation matrix for every time given that the initial condition is an orthogonal rotation, *i.e.* , and and are skew-symmetric matrices.

Proof: We show that for every :

Clearly is a fixed point for the above equation.

Next we prove that the for a given dynamical system, the OTD subspaces that are initially equivalent, remain equivalent for all times.

###### Theorem 2.2

Suppose that and satisfy the evolution equation (12) with different
choices of functions denoted by and respectively.
We also assume that the two bases are initially equivalent, *i.e.
* , where is an orthogonal rotation matrix. Then the subspaces and
are equivalent for with a rotation matrix governed by the matrix differential
equation (15).

Proof: We plug into the OTD equation for :

Multiplying both sides from the right with and using the identity results in:

Now we substitute from equation (15) into the above equation to obtain:

which is the evolution equation for the OTD basis . This completes the proof.

Theorem 2.2 implies that the difference between the two bases will evolve along directions already contained in the initially common subspace. To this end, both bases will continue to span the same subspace and the variation between the two is only an internal rotation. Therefore, the two family of equations will result in the same time-dependent subspace. There are multiple choices for the function and we now examine a special one.

#### 2.1.1 The Dynamically Orthogonal formulation

The simplest choice for the function in (8) is for all . The resulted evolution equations in this case will have the form

(16) |

where is the orthogonal projection operator onto the subspace . Note that corresponds to the dynamical orthogonality (DO) condition [46, 44] that has been employed to derive closed equations for the solution of stochastic PDEs. In this case uncertainty is resolved only along specific modes that evolve with time by projecting the original equation of the system over these directions. The evolution of these modes (stochastic subspace) is done according to equations derived using the DO condition and they have the general form of system (16). The equivalence of system (16) with the minimization problem (6) provides a clear interpretation for the evolution of the DO modes.

### 2.2 Steady linearized dynamics

Here we consider the special case where is a time-independent operator. We prove that the basis defined through the introduced minimization principle will asymptotically span the eigenvectors of associated with the most intense instabilities (i.e. eigenvalues with largest real part). In particular we have the following theorem:

###### Theorem 2.3

Let be a steady and diagonalizable operator that represents the linearization of an autonomous dynamical system. Then

i) Equation (10) has equilibrium states that consist of all the -dimensional subspaces in the span of distinct eigenvectors of .

ii) From all the equilibrium states there is only one that is a stable solution for equation (10). This is given by the subspace spanned by the eigenvectors of associated with the eigenvalues having the largest real part.

Proof: Let where is the matrix of eigenvectors of the operator with the column of being the eigenvectors: , and is the diagonal eigenvalue matrix whose entries are: .

i) First, we show that a subspace that is in the span of
precisely eigenvectors of the operator , is an equilibrium state for equation
(16). Without loss of generality, we consider the first eigenvectors to span such space, *i.e.* associated with , and therefore can be expressed in eigenvector coordinates as:
where is the projection coefficients.
Therefore .
As a result, is in the null space of the orthogonal projector , and thus .
This result is independent of the choice of eigenvectors. To this end, we note that, for an operator with distinct eigenvalues, there exists a number of of such equilibrium states.

ii) Next, we show that from all the equilibrium states , the subspace that is spanned by the eigenvectors associated with the largest eigenvalues, is the only stable equilibrium. First, we investigate the stability of . We denote the complement of the space with associated with the corresponding eigenvalues of . We consider a perturbation that belongs to the orthogonal complement of , i.e. :

We note that the orthonormality condition for , *i.e.* , is satisfied for :

In the above equation, we use the orthonormality condition of that implies: . Moreover, since we immediately obtain:

(17) |

which requires the evolution of the perturbation, *i.e.* , to remain
orthogonal to for all times.
Now, linearizing the evolution equation stated in equation (16) around the equilibrium state yields:

In the above equation, the term , since and . Therefore the evolution equation for the perturbation equation becomes:

(18) |

Next, we transform the evolution equation (18) into the eigenvector coordinates. The perturbation can be expressed in eigenvector coordinates as: , where . Replacing and into equation (18) results in:

(19) |

Having assumed that is non-deficient, we multiply both sides of equation (19) by :

(20) | ||||

where we introduced , and used the orthonormality condition of . We then multiply both sides by to obtain:

(21) |

where . The perturbation can be decomposed as:

(22) |

where and , with . The matrix can be written as:

(23) | ||||

In the above equation, we note that:

(24) | ||||

(25) | ||||

In the last equation, we used the orthonormality condition of and introduced . Replacing equations (24) and (25) into equation (23) results in:

(26) |

By replacing equation (26) into equation (21), the evolution equation becomes:

(27) |

The evolution equation for the perturbation must satisfy the orthogonality constraint expressed by equation (17). The orthogonality constraint requires that:

Since , the above orthogonality condition, in fact, imposes constraints on the evolution of perturbation . In the following, we simplify these constraints. In the above equation:

Therefore, the orthogonality constraint is equivalent to:

(28) |

It follows that

(29) |

or equivalently,

(30) |

Using the orthogonality constraint given by equation (30), and using equation (27), the evolution equation for becomes:

(31) |

The above equation shows that the evolution of can be expressed solely based on the evolution of . Thus, the stability of the solution only depends on the stability of . The evolution of in the index notation is given by:

(32) |

Therefore, the asymptotic stability of requires that:

(33) |

The inequality (33) shows that the subspace with eigenvalues having the largest real part is the only stable solution to equation (16). This completes the proof.

###
2.3 Time-dependent linearized dynamics

Consider the evolution of an arbitrary (not orthonormal) subspace under the time-dependent linearized dynamics which is governed by:

(34) | ||||

and consider the corresponding OTD evolution:

(35) | ||||

We choose the initial state of the OTD basis such that and span the same subspace. In the following Theorem we show that the OTD modes exactly span the subspace . More specifically, we show that the two subspaces and are related via a time-dependent transformation matrix.

###### Theorem 2.4

Let be the subspace evolved under the time-dependent linearized dynamics, eq. (34); and be the OTD basis evolved with eq. (35). Then, assuming that initially the two subspaces are equivalent, i.e. there is a matrix such that , there exists a linear transformation such that where is the solution of the matrix differential equation:

(36) | ||||

Proof: We substitute the transformation into the evolution equation for :

Substituting from equation (36) after rearrangement results in:

Multiplying both sides of the above equation by from the right results in:

which is the OTD equation for . The initial condition is also . This completes the proof.

For dynamical systems with persistent instabilities the evolution under the time-dependent linearized equation (34) is unstable, and grows rapidly. Even for stable dynamical systems, as we move away from , almost any vector approaches to the least stable direction. However the evolution of the OTD modes, due to their built-in orthonormality, is always stable, and as we will demonstrate in our results, the OTD evolution leads to a well-conditioned numerical algorithm that peels off the most unstable directions of the dynamics.

### 2.4 Mode ranking within the subspace

Having derived the subspace that corresponds to the most unstable directions, the next step is to rank these directions internally, i.e. within the subspace. As we describe below, there are two ways to rank the OTD basis based on the growth rate. Both of these approaches amount to an internal rotation within the OTD subspace.

(1) The *instantaneous growth rate* [22, 25] is obtained by computing
the eigenvalues of the symmetric part of , *i.e* :

(37) |

where , and is the matrix of eigenvectors.
We rank these values from the least stable to the most stable directions,
*i.e.*:

We note that is the *numerical abscissa* of the operator .
Therefore represents the maximum instantaneous growth rate within the OTD subspace.

(2) The *instantaneous eigenvectors* of the reduced operator can also be obtained through the equation:

(38) |

where the eigendirections are represented by . The instantaneous eigenvalues are denoted by , where are ranked from the eigenvalue with the largest real part , to the eigenvalue with the smallest real part .

Based on the above two rankings we define the rotated OTD basis:

(39) |

where is the ranked representation of the OTD modes defining the space , and is the rotation matrix obtained from either of the two described strategies for mode ranking (eq. (37) or eq. (38)). For a non-Hermitian operator , the ranking based on the instantaneous eigenvectors results in modes which many not be mutually orthogonal. For this case, the orthogonal representation of the least stable directions can be obtained by performing a Gram-Schmidt orthonormalization.

## 3 Linear dynamics

Linear dynamics and in particular non-normal behaviour has a critical role in determining the short-time fate of a disturbance in both linear and nonlinear dynamical systems. To study the evolution of an initial condition under the effect of
linear or linearized dynamics, reduction in eigenfunction coordinates is often utilized. However for highly non-normal operators, a large number of eigenfunctions are required to correctly capture the non-normal behavior, since the eigenfunctions are far from orthogonal and, in fact they constitute a highly ill-conditioned basis (*i.e.* In what follows we demonstrate the computational efficiency of the OTD modes in capturing non-normal behavior and contrast the OTD basis with the eigenfunction basis for linear systems.

### 3.1 Advection-diffusion operator

As the first case, we consider the advection-diffusion operator which has a wide range of applications in fluid mechanics, financial mathematics, and many other fields. Particularly we are interested to study the effect of non-normality on the reduced operator, both in the short-time and the long-term asymptotic behavior. The operator with zero boundary condition is given by:

(40) | ||||

where is the diffusion coefficient and is the advection speed. The evolution equation for the OTD basis is expressed by:

(41) | ||||

where the inner product and the induced norm are:

To solve the evolution equation we use Chebyshev collocation discretization
implemented in
*chebfun* [1]. The OTD basis is initialized with orthonormalized
modes:

for all the cases considered in this section.

In the case of a nearly normal operator, *i.e* large , an optimal
basis must be close to the dominant eigenfunctions for all times. On the
other hand,
for a set of parameters that corresponds to a non-normal operator, an optimal basis should differ from the eigenfunctions
for short-time dynamics, and only for should it converge to the operator
eigenfunctions. In the remaining of this section, we demonstrate that the
OTD subspace captures short- and long-time dynamics for both normal and
non-normal operators.

To analyze the behavior of the OTD basis we choose . The instantaneous eigenvalues of the reduced operator, , for are shown with colored solid lines in Figure 2(a). The real part of the least stable eigenvalues of the operator are shown with the dashed lines. It is clear that the instantaneous eigenvalues quickly approach the least stable eigenvalues of . In Figure 3(a), the first two elements of the OTD basis, internally rotated in the eigendirection, and , are shown at and . These are superimposed with the two most unstable eigenfunctions of the operator. At , the modes are very close to the eigenfunctions of and at the modes have essentially converged to the eigenfunctions. This demonstrates that due to the strongly normal behavior of the operator at , the eigenfunctions explain the dynamics accurately for both short time and long time, and the OTD basis quickly converges to the space spanned by the eigenfunctions.

At , however, the instantaneous eigenvalues converge with a much
slower rate and much later, *i.e.* , to the operator eigenvalues,
as it is shown
in Figure 2(b). Accordingly, as it can be seen in Figure
3(b), the OTD basis elements are different from the eigenfunctions
of
at , and only until later the basis approaches the eigenfunctions.
Another important observation is related with the advection direction of the OTD basis, which is left-ward. For instance at , the basis has advected an approximate distance of to the left. As a result the OTD basis
has a strong presence in the region . This is not the case for the eigenfunctions
of which are primarily concentrated near .

*i.e.*and : (a) ; (b) .

### 3.2 Orr-Sommerfeld/Squire operator

As the second case we consider the Orr-Sommerfeld/Squire (OS/SQ) equation that governs the evolution of infinitesimal disturbances in parallel viscous flows. The eigenvalues of OS/SQ operator are highly sensitive to perturbations, and its eigenfunctions are linearly dependent, resulting in a highly ill-conditioned linear dynamical system. To this end, the OS/SQ equation is considered only for demonstration purposes, i.e. to illustrate that the OTD modes capture correctly the short-time evolution of the infinitesimal disturbances, as well as their asymptotic (long-term) behavior.

We consider the plane Poiseuille flow in which the base-flow velocity is uni-directional given by , with . The disturbance is taken to be:

with

where and are the wall-normal velocity and the vorticity, respectively, and and denote the streamwise and spanwise wave numbers, respectively.

The Orr-Sommerfeld equation in velocity-vorticity formulation is given by:

(42) |

with boundary conditions:

(43) |

where is a linear operator:

(44) |

with:

and is the modulus of the wave vector and . For the complete derivation of the OS/SQ equations, we refer to [51].

For the choice of the inner product we use the energy measure, which provides a physically motivated norm that arises naturally from the OS/SQ equation [50]. The inner product is given by:

(45) |

where:

(46) |

and denotes complex conjugate. In the following we consider a discrete representation of the operator . Assuming a solution of the form: