# Instability of the optimal edge trajectory in the Blasius boundary layer

Miguel Beneitez, Yohann Duguet, Philipp Schlatter, Dan Henningson

► To cite this version:

Miguel Beneitez, Yohann Duguet, Philipp Schlatter, Dan Henningson. Instability of the optimal edge trajectory in the Blasius boundary layer. Journal of Fluid Mechanics, 2023, 971 (24), pp.A42. 10.1017/jfm.2023.545. hal-04284983

HAL Id: hal-04284983

https://hal.science/hal-04284983v1

Submitted on 6 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Instability of the optimal edge trajectory in the Blasius boundary layer

Miguel Beneitez¹,†, Yohann Duguet², Philipp Schlatter¹ and
Dan S. Henningson¹

¹FLOW Centre and Swedish e-Science Research Centre (SeRC), KTH Engineering Mechanics, Royal Institute of Technology, SE-100 44 Stockholm, Sweden

²LISN-CNRS, Campus Universitaire d'Orsay, Université Paris-Saclay, F-91400 Orsay, France

(Received 22 June 2022; revised 14 April 2023; accepted 21 May 2023)

In the context of linear stability analysis, considering unsteady base flows is notoriously
difficult. A generalisation of modal linear stability analysis, allowing for arbitrarily
unsteady base flows over a finite time, is therefore required. The recently developed
optimally time-dependent (OTD) modes form a projection basis for the tangent space.
They capture the leading amplification directions in state space under the constraint that
they form an orthonormal basis at all times. The present numerical study illustrates
the possibility to describe a complex flow case using the leading OTD modes. The
flow under investigation is an unsteady case of the Blasius boundary layer, featuring
streamwise streaks of finite length and relevant to bypass transition. It corresponds to
the state space trajectory initiated by the minimal seed; such a trajectory is unsteady,
free from any spatial symmetry and shadows the laminar-turbulent separatrix for a finite
time only. The finite-time instability of this unsteady base flow is investigated using the
8 leading OTD modes. The analysis includes the computation of finite-time Lyapunov
exponents as well as instantaneous eigenvalues, and of the associated flow structures. The
reconstructed instantaneous eigenmodes are all of outer type. They map unambiguously
the spatial regions of largest instantaneous growth. Other flow structures, previously
reported as secondary, are identified with this method as relevant to streak switching
and to streamwise vortical ejections. The dynamics inside the tangent space features both
modal and non-modal amplification. Non-normality within the reduced tangent subspace,
quantified by the instantaneous numerical abscissa, emerges only as the unsteadiness of
the base flow is reduced.

**Key words:** transition to turbulence, boundary layer stability

† Present address: DAMTP, Centre for Mathematical Sciences, Wilberforce Road,
Cambridge CB3 0WA, UK. Email address for correspondence: mb2467@cam.ac.uk

© The Author(s), 2023. Published by Cambridge University Press. This is an Open Access article,
distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/
licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original
article is properly cited.
# 1. Introduction

Hydrodynamic stability theory aims at characterising the stability of a given base flow to infinitesimal or finite-amplitude disturbances. In most academic cases, the base flow of interest is known analytically and is generally independent of time (Drazin & Reid 2004). There are, however, physical contexts in which the choice of a physically relevant base flow is not obvious. Bypass transition to turbulence in shear flows falls into this category: there is ample experimental and numerical evidence that turbulent fluctuations emerge from the breakdown of laminar streamwise streaks of sufficiently strong amplitude (Morkovin 1969; Henningson, Lundbladh & Johansson 1993; Jacobs & Durbin 2001; Matsubara & Alfredsson 2001; Brandt & Henningson 2002; Brandt, Schlatter & Henningson 2004) rather than from the destabilisation of the steady laminar base flow. Streamwise streaks, originally called Klebanoff modes, are loosely defined as spanwise modulations of the streamwise velocity field (Klebanoff, Tidstrom & Sargent 1962). They are predominantly streamwise-independent structures supporting three-dimensional wiggles convected at different velocities (Brandt et al. 2003). Streaks are not associated mathematically with unstable eigenmodes of the purely laminar base flow, instead they emerge because of the non-normality of the associated linear operator (Schmid & Henningson 2001) via a mechanism called lift-up. This mechanism transfers streamwise vorticity upstream into streaks further downstream (Landahl 1980; Brandt 2014). Careful early experiments have suggested that their breakdown follows an instability mechanism (Bakchinov et al. 1995; Alfredsson & Matsubara 1996; Matsubara & Alfredsson 2001; Asai, Minagawa & Nishioka 2002). The exact temporal dynamics of finite-amplitude streaks is, however, not trivial. In several numerical studies, a frozen (two-dimensional) finite-amplitude streak pattern was considered as a base flow, and its linear stability analysis was carried out by assuming that the perturbations are inviscid (Andersson et al. 2001; Kawahara et al. 2003; Brandt 2007). The unstable eigenfunctions identified break the translational invariance of the initial streaks. The main outcome of the stability analysis of streamwise-invariant streaks is the possibility for two different ways of breaking this streamwise invariance, either by symmetric (varicose) or anti-symmetric (sinuous) eigenmodes. Around that time, Hamilton, Kim & Waleffe (1995) made use of the concept of subcritical streak instability to justify the three-dimensionality of the self-sustaining process in all shear flows (Hamilton et al. 1995; Waleffe 1997). Hœpffner, Brandt & Henningson (2005), following Schoppa & Hussain (2002), showed that streamwise modulations of the streaks observed during transition, although possible as a linear instability of the frozen streaks, can also arise for lower streak amplitudes via non-normal amplification of streak disturbances over a finite time. In a related study, the secondary instability of time-dependent streaks in channel flow was addressed by adopting a finite-time formalism by Cossu, Chevalier & Henningson (2007). Schlatter et al. (2008) studied the secondary instability of streaks via nonlinear impulse response. Linear stability features were later extracted directly from numerical data (Vaughan & Zaki 2011; Hack & Zaki 2014) by considering an instantaneous streamwise-independent base flow. More recently, the stability of streaks in turbulent flows was also considered by focusing on the associated mean flow rather than on instantaneous flow fields (Alizard 2015; Cassinelli, de Giovanetti & Hwang 2017). It remains hence an open question whether there are additional insights for stability analysis by considering fully unsteady three-dimensional base flows. This paper is devoted to a computational exploration of the possibilities offered by this approach.

In the context of initial value problems, an initial condition at time $t = t_0$ is represented by a point in the associated state space. The knowledge of a given initial condition defines uniquely the base flow, i.e. the unsteady state space trajectory initiated by that particular
initial condition. In principle, the arbitrary unsteadiness of the base flow is not an obstacle to modal linear stability analysis (LSA), at least when the base flow corresponds to an attractor defined over unbounded times. The generalisation of eigenvalues is given by (time-independent) Lyapunov exponents (LEs), defined as ergodic averages of the instantaneous divergence rate between trajectories (Viana 2014). The generalisation of the eigenvectors is given by the (time-dependent) covariant Lyapunov vectors (CLVs); see Ginelli et al. (2007), Kuptsov & Parlitz (2012) and Pikovsky & Politi (2016). Eventually, in the present study, an additional theoretical limitation is the requirement that the method be applicable to a base flow defined only over a finite time interval. This requirement is made necessary by the convective nature of the boundary layer and the fact that any spatially localised perturbation to the Blasius flow has to exit a bounded computational domain in a finite time. In this context, most infinite-time concepts such as eigenvalues need to be formally redefined over the finite time interval of interest. While this does not pose any strong mathematical difficulty, it crucially determines the mathematical toolbox relevant for that problem.

We are interested here in a base flow featuring streamwise streaks of finite length and width, with an unsteady dynamics. Since we wish to define the base flow in an unambiguous way, it is initialised at $t=0$ from a well-defined finite-amplitude perturbation to the original laminar Blasius flow. In the present context of identifying the mechanisms allowing for transition from a minimal level of disturbance, the selected initial condition is the laminar base flow, perturbed at $t=0$ by the so-called minimal seed (Kerswell 2018). The minimal seed is defined rigorously as the disturbance of lowest energy capable of triggering turbulence, or equivalently the point on the edge manifold closest to the laminar attractor in energy norm (Kerswell 2018). Its computation is based on a nonlinear optimisation method (Pringle & Kerswell 2010; Cherubini et al. 2011b; Vavaliaris, Beneitez & Henningson 2020) and in practice requires an optimisation time interval (0, $T_{opt}$). The trajectory initiated by this flow field is called optimal edge trajectory. By construction, it is an edge trajectory i.e. it belongs to the invariant set called the laminar-turbulent boundary: some infinitesimal perturbations to such trajectories lead to relaminarisation while others trigger turbulent flow. The concept of edge trajectory was originally introduced in bistable parallel shear flows (Itano & Toh 2001): the asymptotic fate of such edge trajectories form the edge state, a relative attractor in state space, whose stable manifold divides the state space in two disjoint and complementary basins. Its extension to boundary layer flows is trivial for parallel boundary layer flows (Biau 2012; Khapko et al. 2013, 2014, 2016) but less straightforward in spatially developing boundary layer flows like the Blasius boundary layer (Cherubini et al. 2011b; Duguet et al. 2012). In such cases, the concept of a turbulent attractor is not clearly defined, yet edge trajectories can still be identified, at least over finite times. In boundary layers, the edge concept becomes fragile on very long time scales because the laminar Blasius flow can develop instabilities to Tollmien-Schlichting waves over long time horizons (Beneitez et al. 2019; Beneitez, Duguet & Henningson 2020a). In the absence of an asymptotic state, the stability of finite-time edge trajectories cannot be investigated using LEs and CLVs, all based on ergodic infinite-time averages. The generalisations of eigenvalues/LEs on finite times are, trivially, the finite-time LEs (FTLEs). Their large-time limits, when they are defined, coincide indeed with LEs (Haller 2015). The eigenvectors do not, however, admit any simple finite-time generalisation.

We chose for this task the optimally time-dependent (OTD) modes introduced recently by Babaee & Sapsis (2016). The associated formalism has two advantages: it computes physically meaningful directions in the tangent space, and yields accurate numerical
estimates of the FTLEs. The OTD modes approximate the linearised dynamics Blanchard & Sapsis (2019) around the base flow trajectory in an optimal way, yet under the constraint that the modes remain orthogonal at all times. Orthogonality is not a property shared by CLVs. Handling an orthogonal basis is in practice a strong technical advantage over ill-conditioned bases. The trade-off is that the OTD modes do not fulfil the covariance property. Note that, when both are defined, the leading OTD mode still coincides with the leading CLV for sufficiently long times. The reduced linearised operator, obtained by projecting the original operator on the r first OTD modes, can be used to estimate the stability characteristics of the high-dimensional problem, otherwise prohibitively expensive to compute. In particular, the eigenvalues of this reduced-order operator yield an accurate approximation of the FTLEs of the full system (Babaee et al. 2017). Besides, whereas the OTD modes themselves are not interpretable physically, instantaneous eigenmodes can be reconstructed in physical space from the diagonalisation of the reduced-order operator. As shown by Babaee & Sapsis (2016) from specific examples, over shorter time horizons well-initialised OTD modes can capture the non-normality of the underlying dynamics. These properties make OTD modes an interesting tool specifically for transient phenomena. On a technical level, their implementation requires neither solutions of the adjoint system, nor data to be input and no iterative scheme: the OTD modes are computed in real time together with the time-evolving base flow. They, however, need to be initialised at t = 0. There is currently no accepted general way of choosing initial conditions for these modes, although it is expected that past some finite transient time the OTD directions naturally align with the most important directions of the system. The OTD modes have been used recently in several hydrodynamic applications, including the identification of bursting phenomena (Farazmand & Sapsis 2016), the control of linear instabilities (Blanchard, Mowlavi & Sapsis 2019) and the stability of pulsating Poiseuille flow (Kern et al. 2021) as well as for faster edge tracking in high dimension (Beneitez et al. 2020b). The current investigation, motivated by these promising properties, is an opportunity to test a new computational framework for stability calculations considered until now as challenging.

The present study revisits the optimal edge trajectory in the Blasius boundary layer by considering it as the new finite-time base flow, and by determining its stability characteristics using the new finite-time framework offered by OTD modes. In particular, the physical structure of the leading modes will be analysed at different times, with a focus on the influence of the time dependence of the base flow on the results.

The structure of this paper is as follows. The OTD modes are introduced mathematically in a general context in § 2. The computational set-up, the implementation and the details of the reference edge trajectory are described in § 3. Section 4 contains the stability analysis using the proposed methodology. Finally, the conclusions are given and discussed in § 5.

## 2. Theoretical framework

### 2.1. Linearisation around an arbitrary base flow

The context of the current study is very general. Assuming that a spatially discretised flow field can be represented by n independent real-valued degrees of freedom with n ≫ 1 (see e.g. Gibson, Halcrow & Cvitanović 2008), we consider ℝⁿ as the original high-dimensional space of reference. We suppose a non-autonomous dynamical system defined over a time interval [t₀, t₁)

$$ \frac{dQ}{dt} = f(Q, t), \quad Q(t_0) = Q_0, \qquad (2.1a,b) $$
Instability of the optimal edge trajectory

where $Q_0, Q \in \mathbb{R}^n$, and $f: \mathbb{R}^n \to \mathbb{R}^n$ is a diffeomorphism. We suppose both $t_0$ and $t_1$ finite although $t_1 \to +\infty$ is also possible. For a given choice of $Q_0$, we define the solution to (2.1a,b), namely $\tilde{Q}: (t_0: t_1) \to \mathbb{R}^n$ as the base flow whose stability we will now determine.

Let $q(t)$ represent a small perturbation to $\tilde{Q}(t)$, small enough so that the dynamics can be linearised around $\tilde{Q}(t)$ (mathematically $q$ evolves in the tangent space associated with the dynamics). Then $q$ is governed by the linearised equation

$$ \frac{dq}{dt} = L(\tilde{Q}, t)q, \quad L(\tilde{Q}, t) = \nabla_Q f(\tilde{Q}, t), \qquad (2.2a,b) $$

where $\nabla_Q f(\tilde{Q}, t)$ is the $n \times n$ (time-dependent) Jacobian matrix, evaluated along the base flow at time $t$. The OTD modes, to be introduced in § 2.3, form a basis of time-dependent real-valued vectors (a complex-valued definition is also possible, but is not discussed here). They approximate in an optimal way the leading directions of the Jacobian operator. They are better understood after the notion of covariant vectors is discussed in § 2.2.

## 2.2. Covariance property

An ideal basis for the linearised dynamics should allow one to split the whole $n$-dimensional tangent space into a direct sum of subspaces evolving along the flow, each one with its own specific dynamics (Pikovsky & Politi 2016). The associated time-dependent directions spanning these subspaces are referred to as dynamically covariant. If the base flow $\tilde{Q}$ does not depend on time, the covariance property classically defines expanding and contracting eigenspaces, the covariant vectors are the associated eigenvectors and that the temporal rate of change of their norm defines the eigenvalues. In the general case, these vectors are called CLVs, or sometimes simply Lyapunov vectors. By definition, CLVs can be re-interpreted as zeros of the functional

$$ \mathcal{J} = \lim_{\delta t \to 0} \frac{1}{(\delta t)^2} \sum_{i=1}^{n} \| w_i(t + \delta t) - (\nabla F_t^{t+\delta t}) w_i(t) \|^2, \qquad (2.3) $$

where $F_t^{t+\delta t}: \mathbb{R}^n \to \mathbb{R}^n$ is the infinitesimal forward propagator associated with (2.1a,b). The Jacobian matrix $\nabla F_t^{t+\delta t}$ maps a vector of the tangent space $w_i(t)$ at time $t$ to its image at the later time $t + \delta t$ in the corresponding tangent space.

The issues associated with CLVs are twofold. First, they are not necessarily mutually orthogonal at a given time, making the associated basis possibly ill conditioned. Second, the only algorithms known to compute them are proven to be valid only on attractors on which the dynamics is ergodic (Ginelli et al. 2007; Kuptsov & Parlitz 2012); CLVs are thus essentially an inappropriate computational tool for the study of transients.

## 2.3. Optimally time-dependent modes

Babaee & Sapsis (2016) defined the OTD modes $u_i$, $i = 1, \dots, r$ as a computational compromise. These are minimisers of the functional $\mathcal{J}$ in (2.3), under the additional constraint that they form an orthonormal basis at all times. The orthonormality constraint simply reads

$$ \langle u_i(t), u_j(t) \rangle = \delta_{ij}, \quad i, j = 1, \dots, r, \ r \ll n, \qquad (2.4) $$

where $\langle \cdot, \cdot \rangle$ denotes the inner product associated with the $L^2$ norm and $\delta_{ij}$ is the classical Kronecker symbol. The time evolutions of the OTD modes and of a set of initially random
<figure><img src="image_3.png" /><figcaption>Figure 1. Sketch of different vector bases for the tangent space of a given unsteady base flow. Temporal evolution of respectively OTD modes (red, a) and random unit vectors (orange, b) propagated along a base flow trajectory. The leading non-rescaled direction coincides with the leading OTD mode and is shown pointing upwards at all times.</figcaption></figure>

unit vectors are depicted schematically in figure 1 for illustration purposes. Random unit vectors would follow the tangent dynamics and all align rapidly with the most expanding direction, making them poor candidates to describe and analyse the dynamics of the tangent space. The OTD modes follow the tangent dynamics but stay orthonormal at all times, avoiding any alignment issues which might occur using e.g. CLVs. This constraint destroys the interpretability of each OTD direction in terms of the covariant dynamics. However, the orthonormality of the OTD modes is particularly appealing for reduced-order modelling, for instance in the context of control (Blanchard & Sapsis 2019; Blanchard et al. 2019).

As derived in Babaee & Sapsis (2016), the maximisation of $\mathcal{J}$ in (2.3) under the orthogonality constraint (2.4) yields a system of coupled nonlinear evolution equations

$$ \frac{d\mathbf{u}_i}{dt} = L(t)\mathbf{u}_i - \sum_{j=1}^{r} [\langle L(t)\mathbf{u}_i, \mathbf{u}_j \rangle - A_{ij}(t)]\mathbf{u}_j, \quad i = 1, \dots, r. \quad (2.5) $$

The nonlinearity is a direct consequence from the orthonormality constraint. The matrix $A \in \mathbb{R}^{r \times r}$ refers a priori to any skew-symmetric matrix. The discretised equations (2.5) for $i = 1, \dots, r$ form, together with (2.1a,b), a closed $(r+1)$-dimensional system of real-valued ODEs. The $\mathbf{u}_i$ modes depend on the instantaneous vector $\tilde{Q}(t)$, however, $\tilde{Q}$ itself is unaffected by the evolution of the $\mathbf{u}_i$ modes. In Blanchard & Sapsis (2019), a particular choice of $A$ was made

$$ A_{ij} = \begin{cases} -\langle \mathbf{L} \mathbf{u}_j, \mathbf{u}_i \rangle, & j < i \\ 0, & i = j \\ \langle \mathbf{L} \mathbf{u}_i, \mathbf{u}_j \rangle, & j > i, \end{cases} \quad (2.6) $$

which is also considered here. An arbitrary choice of $A$ would lead to a fully coupled system so that each $\mathbf{u}_i$ appears in every equation in (2.5). However, under the current choice the set of $r$ equations in (2.5) has a lower triangular form: the evolution of the $i$th
Instability of the optimal edge trajectory

mode depends only on the modes from 1 to *i*, making the OTD formulation hierarchical. The resulting evolution equation for each OTD mode becomes

$$ \frac{d\mathbf{u}_i}{dt} = \mathbf{L}(t)\mathbf{u}_i - \langle \mathbf{L}(t)\mathbf{u}_i, \mathbf{u}_i \rangle \mathbf{u}_i - \sum_{j=1}^{i-1} [\langle \mathbf{L}(t)\mathbf{u}_i, \mathbf{u}_j \rangle + \langle \mathbf{L}(t)\mathbf{u}_j, \mathbf{u}_i \rangle] \mathbf{u}_j. \quad (2.7) $$

The nonlinear system of *r* (2.7) can be evolved forward in time together with (2.1a,b), from which the matrix *L* can be evaluated at all times. Equation (2.1a,b) remains, however, independent of the evolution of each *u*<sub>*i*</sub>. This results in an (*r* + 1) × *n*-dimensional asymmetrically coupled dynamical system. The OTD modes retain a (short-time) memory of their initial conditions. They are in general not covariant, except for base flows such that all instantaneous eigenvectors remain normal to each other.

## 2.4. The reduced linearised operator

In order to analyse the linearised dynamics within the reduced subspace optimally spanned by the OTD modes, Babaee & Sapsis (2016) introduced the reduced operator *L*<sub>*r*</sub> defined by projecting the high-dimensional operator *L* onto the OTD directions

$$ L_{r_{ij}}(t) = \langle \mathbf{u}_i, \mathbf{L}(t)\mathbf{u}_j \rangle \quad i, j = 1, \dots, r. \quad (2.8) $$

In particular, all the instantaneous stability indicators defined in the next subsection will be derived from algebraic properties of the *r* × *r* matrix *L*<sub>*r*</sub> evaluated at the relevant times.

For a time-independent linearised operator *L* the space spanned by the modes {u<sub>*i*</sub>}<sub>*i* = 1</sub><sup>*r*</sup> converges asymptotically to the most unstable eigenspace of *L*. Moreover, if *L* happens to be also symmetric, the OTD modes coincide with its eigenvectors at all times.

## 2.5. Instantaneous stability indicators

As emphasised in Babaee & Sapsis (2016), although they correspond to divergence-free vector fields, the *u*<sub>*i*</sub> modes do not have a direct physical interpretation as flow fields. More meaningful vector sets can nevertheless be reconstructed instantaneously from the knowledge of the *u*<sub>*i*</sub> modes and the reduced operator. At every time, *L*<sub>*r*</sub>(*t*) can be diagonalised as *L*<sub>*r*</sub> = *E*<sup>*λ*</sup>Δ<sub>*r*</sub>(*E*<sup>*λ*</sup>)<sup>-1</sup> with *E*<sup>*λ*</sup> and Δ<sub>*r*</sub> = diag(*λ*<sub>1</sub>(*t*), ..., *λ*<sub>*r*</sub>(*t*)). The new modes u<sub>*i*</sub><sup>*λ*</sup>, *i* = 1, ..., *r* are defined (using the summation convention) by

$$ \mathbf{u}_i^{\lambda}(t) = E_{ij}^{\lambda}(t) \mathbf{u}_j(t). \quad (2.9) $$

Unlike the *u*<sub>*i*</sub> modes, the u<sub>*i*</sub><sup>*λ*</sup> modes are not necessarily mutually orthogonal. They are interpreted as instantaneous eigenmodes. They are the only velocity fields used for visualisation in this paper. We emphasise that, although the *u*<sub>*i*</sub> modes are real valued, the modes u<sub>*i*</sub><sup>*λ*</sup> are complex valued and come in pairs. Only the real parts, with arbitrary phase, of the associated velocity fields will be represented. The time-dependent numbers *λ*<sub>1</sub>(*t*), ..., *λ*<sub>*r*</sub>(*t*), *i* = 1, ..., *r* are labelled instantaneous eigenvalues and are complex valued.

Another key scalar quantity is the instantaneous numerical abscissa σ(*t*), a positive number, defined as the largest eigenvalue of the symmetrised reduced operator (L + L<sup>*T*</sup>)/2 (Embree & Trefethen 2005). Here, σ corresponds to the largest possible growth rate at a given time, due to both normal and non-normal effects combined together. Whenever *L* is non-normal, σ is strictly larger than the real part of all *λ*<sub>*i*</sub> values.
M. Beneitez, Y. Duguet, P. Schlatter and D.S. Henningson

For a given integer value $r$, it is a natural extension to define $\sigma_r$ as the largest eigenvalue of the symmetrised reduced operator $(L_r + L_r^T)/2$. The gap

$$g_r(t) := \min_i |\sigma_r - \mathrm{Re}(\lambda_i)| = \sigma_r - \mathrm{Re}(\lambda_1), \qquad (2.10)$$

quantifies the non-normality of the reduced linearised operator at each instant. The values of $g_r(t)$ bound from below the value of $g_n(t)$ corresponding to the full high-dimensional system. In the remainder of the paper we will not make a difference between $g_r$ and $g_n$ and will simply use the notation $g(t)$. Note that for arbitrary time-dependent operators and finite $r$, it is possible to have $\sigma_r \approx \max \lambda_i$ even for a non-normal operator.

## 2.6. Finite-time Lyapunov exponents

For a general dynamical system in dimension $n$, characterised by a propagator $F_{t_0}^t$, the Cauchy–Green tensor $C_{t_0}^t$ is defined as

$$C_{t_0}^t(q_0) := [\nabla F_{t_0}^t(q_0)]^T [\nabla F_{t_0}^t(q_0)]. \qquad (2.11)$$

The associated FTLEs are defined directly from the eigenvalues $\gamma_1 > \gamma_2 > \dots > \gamma_n$ of the Cauchy–Green tensor (Haller 2015). These eigenvalues are real and positive by virtue of the positive definiteness of the Cauchy–Green tensor. Each FTLE is defined, for an initial time $t_0$ and a horizon time $T > 0$ (Haller 2015), as

$$(\Lambda_{t_0}^{t_0+T})_i = \frac{1}{T} \log \sqrt{\gamma_i}, \quad i = 1, \dots, n. \qquad (2.12)$$

Babaee et al. (2017) provided an analytical proof that, for any integer $r > 0$, the $r$-dimensional OTD subspace aligns exponentially fast, i.e. for increasing $T$ with the space spanned by the $r$ most dominant left vectors of the Cauchy–Green tensor. The exact rate of convergence depends on the spectrum of the problem at hand. The OTD formulation is hence a robust direct method to estimate the $r$ leading FTLEs $(\Lambda_{t_0}^{t_0+T})_i$, $i = 1, \dots, r$ of the full system at any time $t_0$ (Babaee et al. 2017; Sapsis 2018; Blanchard & Sapsis 2019), provided the time horizon $T$ is large enough. As a consequence the FTLEs $(\Lambda_{t_0}^{t_0+T})_i$, $i = 1, \dots, r$ are evaluated simply by averaging over time the diagonal elements of $L_r$ (Blanchard & Sapsis 2019). They are expressed as

$$(\Lambda_{t_0}^{t_0+T})_i \approx \frac{1}{T} \int_{t_0}^{t_0+T} \langle u_i(\tau), L_r(\tau) u_i(\tau) \rangle d\tau, \quad i = 1, \dots, r. \qquad (2.13)$$

It can be useful to relate the OTD modes to other known vector sets from the literature beyond the CLVs. The Gram–Schmidt vectors are precisely involved in the classical algorithms used for computing FTLEs and hence LEs (see e.g. Shimada & Nagashima 1979). The OTD modes coincide with the so-called Gram–Schmidt vectors, at least in the limit where the Gram–Schmidt vectors are continuously re-orthogonalised (Blanchard & Sapsis 2019). The same modes have also been called sometimes backwards Lyapunov vectors Kuptsov & Parlitz (2012). Unlike the CLVs, the OTD modes depend on the choice of the inner product except for the leading mode.

# 3. Computational set-up

## 3.1. Direct numerical simulation

The Blasius boundary layer is the incompressible flow over a semi-infinite flat plate. It develops at the leading edge of the plate in the absence of a streamwise pressure gradient.
Let x, y, z denote the streamwise, wall-normal and spanwise directions, respectively. Here, **v** is the total velocity field, **v**<sub>B</sub> = (u<sub>B</sub>, v<sub>B</sub>, 0) that of the steady Blasius solution, then **u** := **v** − **v**<sub>B</sub> = (u, v, w) is the perturbation velocity field. All quantities are made non-dimensional using the free-stream velocity U<sub>∞</sub> and the boundary layer thickness δ*(x) := ∫<sub>0</sub><sup>∞</sup>(1 − u<sub>B</sub>(x)/U<sub>∞</sub>) dy of the undisturbed (steady) Blasius flow. A local Reynolds number can be defined as Re<sub>δ*</sub>(x) := U<sub>∞</sub>δ*/ν, with ν the kinematic viscosity of the fluid. The value of Re<sub>δ*</sub> = Re<sub>δ*</sub>(x = 0) is imposed at the upstream end (x = 0) of the computational domain, located at a finite distance downstream of the leading edge.

The boundary conditions for the edge trajectory at the wall (y = 0) are of no-slip and no-penetration type

$$u = v = w = 0, \tag{3.1}$$

and at the upper domain boundary (y = L<sub>y</sub>) of Neumann type to allow for a natural growth of the boundary layer

$$\frac{\partial u}{\partial y} = \frac{\partial v}{\partial y} = \frac{\partial w}{\partial y} = 0. \tag{3.2}$$

A fringe region located at the downstream end of the domain damps outgoing velocity perturbations consistently with the streamwise periodic boundary conditions. The fringe is imposed as a volume force F(t, x, y, z) of the form

$$F = \gamma(x)(\mathcal{U}(x, y, z) - \mathbf{v}(t, x, y, z)), \tag{3.3}$$

where γ(x) is a non-negative fringe function detailed in Chevalier et al. (2007). The streamwise component of U(x, y, z) is defined as

$$\mathcal{U}_x = U(x, y, z) + [U(x + x_L, y, z) - U(x, y, z)]S\left(\frac{x - x_{\text{blend}}}{\Delta_{\text{blend}}}\right), \tag{3.4}$$

where S(x<sub>blend</sub>, Δ<sub>blend</sub>) is a blending function connecting smoothly the outflow to the inflow, and U(x, y, z) solves the boundary layer equations. The wall-normal component of U is obtained via the continuity equation. In the present work the fringe length is Δ<sub>blend</sub> = 600, x<sub>L</sub>=2500 and γ<sub>max</sub> = 0.8.

The present approach has been successfully applied in most works referenced in Chevalier et al. (2007), and in several later publications including Duguet et al. (2012) and Beneitez et al. (2019). The effect of the fringe on outgoing perturbations, allowing for the simulation of spatially developing flows in the presence of periodic boundary conditions was analysed in full mathematical detail in Nordström, Nordin & Henningson (1999).

The temporal integration of the incompressible Navier–Stokes equations is performed using the pseudo-spectral solver SIMSON (Chevalier et al. 2007). This direct numerical simulation (DNS) code solves the equations in the wall-normal velocity-vorticity formulation. The solution is advanced in time using a second-order Crank–Nicholson scheme for the linear terms and a fourth-order low-storage Runge–Kutta scheme for the nonlinear terms. The timestep is fixed to Δt = 0.2 in terms of U<sub>∞</sub> and δ<sub>0</sub>*. The velocity field is expanded along N<sub>x</sub> Fourier modes in the streamwise direction x and N<sub>z</sub> modes in the spanwise direction z, N<sub>y</sub> Chebyshev modes are used in the wall-normal direction y using the Chebyshev-tau method. The evaluation of the nonlinear terms obeys the 3/2-rule for dealiasing.

The additional equations (2.7) ruling the evolution of the OTD modes are advanced in time using the same scheme, based on an explicit evaluation of the inner products at
every collocation point at every timestep. The initial conditions for the modes $i = 1, ..., r$ are spatially localised disturbance velocity fields, consistent with the localised nature of the perturbations to the streaks observed in bypass transition. The boundary conditions for (2.5) are the same as for the original DNS. The choice of boundary conditions is particularly sensitive for the perturbation equations. Further details can be found in Appendix A. The computational requirements for each individual OTD mode are the same as a full DNS.

Although, as described in Beneitez et al. (2019), it would be possible to use a moving box technique to track localised disturbances over long time horizons using limited computational resources, this is not required here because of the limited tracking time. The reference frame is hence understood as the laboratory frame. The computational set-up for the edge tracking is similar to that in Vavaliaris et al. (2020). The computational domain $\Omega$ has dimensions $[L_x, L_y, L_z] = [2500, 60, 100]$ and the velocity field is expanded on $[N_x, N_y, N_z] = [2048, 201, 256]$ modes before dealiasing. This numerical resolution is comparable locally to that used in Duguet et al. (2012) and Beneitez et al. (2019). The computation of the OTD modes starts at initial time $t = 0$ from the (spatially localised) minimal seed computed in Vavaliaris et al. (2020) corresponding to the state shown in Figure 3. It ends at $t = 800$, at which time the localised perturbation has not yet left the computational domain.

The computation of the OTD modes in (2.5) depends on the definition of the inner product, chosen here as

$$ \langle \mathbf{u}, \mathbf{u}' \rangle = \int_{\Omega} (uu' + vv' + ww') \, d\Omega, \qquad (3.5) $$

where $\mathbf{u} = (u, v, w)$ and $\mathbf{u}' = (u', v', w')$ are any two flow fields with finite $L^2$ norm, and $d\Omega = dx dy dz$ is the usual infinitesimal integration element over the numerical domain $\Omega$.

## 3.2. The optimal edge trajectory

The minimal seed $M$ refers to the perturbation closest in kinetic energy to the laminar Blasius boundary layer flow, and able to trigger subcritical transition. This particular optimal condition was selected because it gives rise to a fully nonlinear trajectory relevant for the method tested. Moreover, it is uniquely defined by the parameters for the optimisation algorithm, namely here the Reynolds number value $Re_0^* = 240.458$. That value of $Re_0(x = 0)$ is chosen to match previous works (Cherubini et al. 2011a; Vavaliaris et al. 2020), in particular the original work by Cherubini et al. (2011a) where the non-dimensionalisation differs from the present one. Note that in parallel flows the Reynolds number entirely defines the dynamical system, however, in spatially developing flows the Reynolds number is intrinsically linked to the streamwise coordinate. Consequently, the minimal seed is conditioned by the range of Reynolds numbers (streamwise distances) allowed for in the time evolution of the perturbations. This results in the minimal seed being dependent on the inlet Reynolds number, on the length of the computational domain and on the optimisation time (Vavaliaris et al. 2020; Beneitez et al. 2020a). In Vavaliaris et al. (2020) the chosen optimisation time is $T_{opt} = 400$ and the computational domain length $L_x = 500$. Here, $M$ is computed iteratively using the nonlinear adjoint-based optimisation framework of Rabin, Caulfield & Kerswell (2012), Kerswell (2018). The maximised objective function is the energy gain at a given time $T_{opt}$, $G(T_{opt}) = E(T_{opt})/E(0)$, where $E(t)$ is the perturbation kinetic energy at time $t$. The optimisation framework follows Foures, Caulfield & Schmid (2013) and is based on the
Instability of the optimal edge trajectory

<figure><img src="image_4.png" /><figcaption>Figure 2. State portrait of the optimal edge trajectory in the Blasius boundary layer using the variables $\Omega_x$, $\Omega_y$ and $W$ defined by (3.6)–(3.8). The edge trajectory starts at $t = 0$ from the minimal seed $M$ computed in Vavaliaris et al. (2020). The part of the edge trajectory (green) investigated in this work, $t \in [0, 800]$, is shown using a thicker line. The other trajectories start from the neighbourhood of $M$, they bracket the edge trajectory, and approach either the laminar (blue) or the turbulent state (red). Dots are plotted every 50 time units, highlighting the slowdown already after $t = 100$. Here, $M$ indicates the location of the minimal seed.</figcaption></figure>

implementation into the open-source solver Nek5000 originally implemented by Rinaldi, Canton & Schlatter (2019). The optimal state determined for a near-to-threshold initial energy $E_0$ was bisected using an edge tracking algorithm (Itano & Toh 2001; Skufca, Yorke & Eckhardt 2006), so that the computed trajectory approximates well an edge trajectory for $t \le 800$, the bracketing trajectories differing by less than 2% in the observable used for edge tracking. This property is crucial for the stability study: initialising the base flow for the OTD analysis from outside the edge manifold would possibly result in a different transition scenario, as reported e.g. in Cherubini et al. (2011a). Although the investigation in Beneitez et al. (2019) warned against the possible interference between edge trajectories and unstable Tollmien–Schlichting waves over time scales $O(10^4)$, no such phenomenon will be encountered with the present set-up, since the considered observation time is $O(10^3)$.

A state portrait is shown in figure 2, based on the three global quantities already used in previous studies (Duguet et al. 2012; Beneitez et al. 2019, 2020a)

$$ \Omega_x = (\delta_0^*/\delta)^{1/2} \left( \frac{1}{V} \int_V |\omega_x|^2 \, dv \right)^{1/2}, \qquad (3.6) $$

$$ \Omega_y = (\delta_0^*/\delta)^{1/2} \left( \frac{1}{V} \int_V |\omega_y|^2 \, dv \right)^{1/2}, \qquad (3.7) $$

$$ W = (\delta_0^*/\delta)^{3/2} \left( \frac{1}{V} \int_V |w|^2 \, dv \right)^{1/2}. \qquad (3.8) $$

The quantities $\omega_x$ and $\omega_y$ are the streamwise and wall-normal perturbation vorticity components, respectively, and the integration is carried over the computation domain of volume $V$. The prefactors in powers of $(\delta_0^*/\delta)$ make use of the value of the boundary layer thickness evaluated at the centre of mass, see Duguet et al. (2012). In figure 2, the edge trajectory is highlighted using a thicker (green) line, and the time interval $t \in [0, 800]$
<figure><img src="image_5.png" /><figcaption>Figure 3. Three-dimensional top view of the perturbation velocity field corresponding to the initial condition of the reference trajectory. Contours of $u = -8 \times 10^{-3}$ (blue), $u = 6 \times 10^{-3}$ (red) and global $\lambda_2^* = -4 \times 10^{-4}$, where $\lambda_2^*$ denotes the vortex identification criterion introduced by Jeong & Hussain (1995). The streamwise extension of the minimal seed is ~22 length units.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 4. Three-dimensional top view of the perturbation velocity field reference trajectory at (a) $t = 100$, contours of $u = -8 \times 10^{-3}$ (blue), $u = 6 \times 10^{-3}$ (red) and global $\lambda_2^* = -1 \times 10^{-4}$, where $\lambda_2^*$ denotes the vortex identification criterion introduced by Jeong & Hussain (1995). The distance between the vertical lines is $\Delta x = 50$ (b) $t = 280$, contours of $u = -5 \times 10^{-2}$ (blue), $u = 6 \times 10^{-2}$ (red) and $\lambda_2^* = 1 \times 10^{-4}$. The distance between the vertical lines is $\Delta x = 50$ and (c) $t = 720$, contours of $u = -5 \times 10^{-2}$ (blue), $u = 6 \times 10^{-2}$ (red) and global $\lambda_2^* = 1 \times 10^{-4}$. The distance between the vertical lines is $\Delta x = 100$.</figcaption></figure>

considered in this study is highlighted using a thicker green line (with equispaced dots every 50 time units). The thinner lines in red and blue correspond to trajectories closely bracketing the edge trajectory.

It is useful to recall the main features of the unsteady base flow reported by Vavaliaris et al. (2020). For early times $t \le 60$ the dynamics is dominated by a three-dimensional version of the Orr mechanism (Vavaliaris et al. 2020), where vortical disturbances initially tilted against the mean shear progressively untilt as time increases. For $60 \le t \le 200$ the lift-up mechanism takes over and a pair of streamwise streaks forms. Both mechanisms are known to be non-modal, the stronger energy amplification being associated with the lift-up
Instability of the optimal edge trajectory

<figure><img src="image_7.png" /><figcaption>Figure 5. Real part of the instantaneous eigenvalues Re($\lambda_i$), $i = 1, \dots, 8$ (solid lines) and instantaneous numerical abscissa $\sigma$ (dashed red) of the reduced operator $L_r$, plotted vs time. Here, $\lambda_1$ is highlighted (orange line). Vertical dashed line indicates $t = 50$.</figcaption></figure>

(Schmid & Henningson 2001). For $t \ge 200$ the energy growth slows down. Snapshots of the velocity field along the optimal edge trajectory are shown in figure 4 at times $t = 100, 280$ and $720$. From $t \ge 100$ onwards, the edge trajectory consists of a localised pair of high- and low-speed streaks (Beneitez et al. 2019) with an undulation linked to oblique waves. It experiences a couple of streak-switching events around $t \approx 500$ and $t \approx 700$. The streaks elongate with time but remain always localised in the streamwise and spanwise directions. By construction, typical infinitesimal perturbations of this unstable flow field will make it evolve either towards an incipient turbulent spot or towards the laminar state. It is precisely their state space location on the verge of bypass transition that makes edge trajectories a relevant choice as a base flow (Khapko et al. 2016). Imposing an optimality condition has the advantage of making the current trajectory well defined.

# 4. Results

This section is devoted to the analysis of the stability properties of the optimal trajectory described in § 3 using $r = 8$ OTD modes. The choice of 8 modes aims at producing the largest possible subspace while keeping the simulations computationally feasible. The cost of each OTD mode is comparable to an additional DNS to be run in parallel to the original base flow. Moreover, the choice of number of modes is comparable to previous simulations of similar scale (Babaee & Sapsis 2016). We restrain our study to the time interval $t \in [0, 800]$.

## 4.1. Finite-time stability analysis

### 4.1.1. *Instantaneous growth rates*

We begin by reporting the real part of the instantaneous eigenvalues $\lambda_i(t)$, $i = 1, \dots, r$, computed over the whole trajectory. They are shown together with the instantaneous numerical abscissa $\sigma(t)$ vs time in figure 5. The gap $g(t) = \sigma(t) - \text{Re}(\lambda_1)$, which quantifies the instantaneous non-normality of the reduced operator, is displayed as a black line in figure 6. These quantities have all been defined in § 2.5.

The time series of these instantaneous growth rates can be grossly divided into two phases. In the initial phase for $t \lesssim 100$, the two leading growth rates vary rapidly in
<figure><img src="image_8.png" /><figcaption>Figure 6. Quantifier $g(t)$ of the non-normality in the reduced system (black). Norm of the vector formed by the time derivatives of the observables used in the state portrait of figure 2: $d\Omega_x/dt$, $d\Omega_y/dt$ and $dW/dt$ as defined in (4.1) (red). The figure starts at $t = 50$.</figcaption></figure>

time while the others are all negative. In a second phase starting at $t \approx 100$, Re($\lambda_1$) dominates in the range 0.2–0.3, with a slight decaying trend as time increases. All other eigenvalues remain close to zero in real part, never exceeding 0.1. A quick glance at the state portrait in figure 2 suggests that this second phase corresponds to a clear slowdown of the dynamics of the base flow itself. If the dynamics is quasi-steady, it is expected that the stability properties of the edge trajectory mimic qualitatively the stability properties of steady/travelling edge states reported in other shear flow studies: one large dominating unstable eigenvalue, representing a strong instability in a direction transverse to the edge manifold, associated with many other eigenvalues of lesser magnitude responsible for the slow chaotic fluctuations within the edge manifold (Duguet, Willis & Kerswell 2008). This expectation is largely confirmed by figure 5 for $t \ge 100$.

A finer analysis of the fluctuations of the growth rates is possible both in the initial and the quasi-steady phases. This is achieved by focusing on the gap $g(t)$, interpreted as a measure of instantaneous non-normality within the OTD subspace. For the initial times $t \lesssim 50$, $\sigma = \text{Re}(\lambda_1) = \text{Re}(\lambda_2) > 0$. After $t = 50$, the gap $g$ rises from zero to a maximum of approximately 0.4. It later decreases to smaller values of $\approx 0.1$. As for the other $\lambda_i$ values, they are all negative at $t = 0$ but grow at the same pace and cross zero at $t \approx 100$. At later times, all instantaneous growth rates stabilise, while Re($\lambda_1$) decreases gently in a non-monotonic manner, and $g(t)$ oscillates around low values $\approx 0.05–0.1$. The fact that the peak of $g(t)$ occurs before $t = 100$ is consistent with the reported occurrence of purely non-normal Orr and lift-up mechanisms along the edge trajectory for these times (Vavaliaris et al. 2020). The sensitivity of the edge trajectory appears high where the edge trajectory also experiences strong non-normal amplification. However, the fact that $g \approx 0$, i.e. $\sigma = \lambda_1$ at the earliest times $t \le 50$, may be wrongly attributed to a lack of non-normal potential of $L(t)$. To start with, this is a property of the instantaneous reduced operator $L_r(t)$ computed for a given value of $r$, not necessarily of the full operator $L(t)$. The reverse is yet true: non-normal features of the reduced-order operator $L_r(t)$ carry over to $L(t)$. Moreover, after trying several different initialisations this result was found to depend crucially on the choice of the OTD basis for $t = 0$, at least over early times $t \le 50$. This makes it difficult to draw general conclusions for short enough times,
consistently with the study of Babaee et al. (2017). This is possibly confirmed by the very transient behaviour of the eigenvalues $\lambda_2$ to $\lambda_8$. From $t \approx 50$ on, the non-normal potential within the OTD subspace is high again as expected, judging from the large values of $g(t)$, and transient effects due to the initialisation of the OTD modes can be neglected.

A peak at $t \approx 60$, and a smaller one at $t \approx 100$, are evident in the data for $\sigma(t)$ in figure 5(a). These times are perfectly consistent with the occurrence of both the Orr and the lift-up mechanisms described in Vavaliaris et al. (2020). Two additional bumps for both $g(t)$ and $\sigma(t)$ can also be seen at $t \approx 550$ and $t \approx 720$. According to Vavaliaris et al. (2020), these two times correspond to streak-switching events. This suggests that streak-switching events, themselves an inherent part of the self-sustained mechanism (Khapko et al. 2013; Beneitez et al. 2019), are linked to stronger non-normality than the rest of the edge trajectory.

Eventually, figure 6 also shows the norm of the time derivatives of the three observables $\Omega_x$, $\Omega_y$ and $W$ used in figure 2. This quantity is defined as $\xi$ via

$$ \xi(t) = \sqrt{\left(\frac{d\Omega_x}{dt}\right)^2 + \left(\frac{d\Omega_y}{dt}\right)^2 + C \left(\frac{dW}{dt}\right)^2}, \quad (4.1) $$

where $C$ is a unity-valued constant ensuring the correct dimensionality. It is plotted in figure 6 in connection with the time evolution of $g(t)$. We analyse now these quantities by considering consecutive sub-intervals of the edge trajectory starting from the minimal seed: (i) $t \in [0, 60]$ (Orr mechanism in the base flow) corresponds to a very rapid evolution of the observables $\xi(t)$. The OTD modes, however, take time to catch up with non-normality until $t \approx 80$, as shown by $g(t)$; (ii) $t \in [60, 200]$ corresponds to the lift-up in the base flow effect associated with non-normal growth. Here, $g(t)$ appears largest for $t \approx 80$ and decreases rapidly until $t \approx 130$, where a change in the slope of $g(t)$ can be noticed. The value of $\xi(t)$ mirrors this behaviour, suggesting that non-normality is decreasing as the lift-up of the base flow ends; (iii) the trajectory has reached the relative attractor past $t \ge 200$. In this stage we observe that the slowdown of the dynamics indicated by $\xi(t)$ corresponds to higher values of $g(t)$, and vice versa. This can be seen in the intervals $t \in [300, 400]$, where the dip in $g(t)$ corresponds to a peak in $\xi(t)$, and in $t \in [500, 600]$, where an increase in $g(t)$ corresponds to a dip in $\xi(t)$.

### 4.1.2. Characterisation as an outer mode instability

In the original study on streak breakdown by Vaughan & Zaki (2011), where the base flow consists of a quasi-steady localised streak rather than a time-dependent one, a distinction was made between two types of modes. The main criterion is the wall-normal position of the energy of each mode with respect to the location of the critical layer, the latter being known from inviscid analysis. The modes with a critical layer close to the wall (such as Orr-Sommerfeld modes) are denoted as inner modes, while those with a critical layer in the free stream are denoted as outer modes. Another characterisation of the inner vs outer mode distinction, also suggested by Vaughan & Zaki (2011), relies on the relation between the growth rate of the mode and the streak amplitude. Although the present context differs, notably because of the unsteady aspect of the streaks, such a characterisation can also be applied to the modes determined by our method. Figure 7 shows the values of the two largest instantaneous eigenvalues Re($\lambda_{1,2}$) plotted vs $\Omega_x$, in figure 7(a), and vs the volume average energy in the spanwise direction $||w||^2$. These quantities are used as a proxy for the instantaneous amplitude of the streaky edge state. It can be directly compared with figure 7
<figure><img src="image_9.png" /><figcaption>Figure 7. (a) Green lines: Re($\lambda_{1,2}$) vs $\Omega_x$, red line $\sigma$. The dots denote $t = 50$, the trajectory starts from the values on the left of the figure. (b) Right: *idem* for Re($\lambda_{1,2}$) vs $||w||^2$, the volume average spanwise velocity.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8. Histograms of all $r = 8$ reduced-order FTLEs ($\Lambda_{t_0}^{t_0+T}$)$_i$ ($i = 1, \dots, 8$) sampled over various time window values of $t_0 \in [0, 720]$ along the unsteady base flow trajectory. Horizon times $T = 10, 20, 30, 50$ and $70$. The histograms are normalised so that their integral is equal to 1; (a) $t_0 \in [0, 720]$, (b) $t_0 \in [100, 720]$.</figcaption></figure>

from Vaughan & Zaki (2011), where the growth rate is plotted vs the streak amplitude (called $A_u$). The corresponding figure was used to define a classification of the instability mechanisms: inner mode refers to an instability mode present for arbitrary small values of $A_u$, in contrast with outer modes which are not found for vanishing streak amplitude. In the present case, positive growth rates are only found for non-vanishing values of the observable $\Omega_x \ge 0.05$. Interpreting $\Omega_x$ as an alternative definition of streak amplitude unambiguously indicates that the dominant instability of the edge state should be classified as an outer mode instability.

### 4.1.3. Finite-time Lyapunov exponents

Figure 8(a) shows distributions of FTLEs ($\Lambda_{t_0}^{t_0+T}$)$_i$ ($i = 1, \dots, 8$) computed within the interval $t_0 \in [0, 800]$. Figure 8(b) is similar except that the values of $t_0$ are restrained to
Instability of the optimal edge trajectory

the sub-interval $t_0 \in [100, 800]$. In both plots the time horizon $T$ takes increasing values from 10 to 70. Comparing the different values of $T$ essentially confirms the robustness of the FTLE distributions with respect to the time horizon. The reason why all FTLES from $i = 1$ to 8 are reported together is the frequent change in the ordering of the growth rates, occurring every time an eigenvalue crossing takes place (Babaee et al. 2017). The many negative occurrences in figure 8(a), as well as the largest occurrences ($\ge 0.3$) can be attributed to the choice of initial conditions for the OTD modes, including accidentally co-aligned disturbances. A potential improvement of the initial conditions could be the computation of the eigenvectors associated with the minimal seed, by assuming no time dependency. This would result in eigendirections already within the initial tangent space. Even though there is no guarantee that these directions will remain in the tangent space at later times, they can be expected to be physically relevant at least for the initial times. These occurrences indeed disappear entirely in figure 8(a) after the initial first 100 time units have been discarded, consistently with the results of § 4.1.1. Since the original bisection algorithm is essentially a shooting method (Itano & Toh 2001), we expect one of the FTLEs to be the signature of the instability of the edge manifold. In other words this FTLE is associated with an unstable direction pointing transversally to it. The other additional positive FTLEs have no choice but to be associated with the weak apparently unsteady dynamics taking place within the relative attractor, rather than transversally to it. This conclusion is consistent with the results of § 4.1.1. The two peaks in figure 8(a), close to 0.15 and 0.3, correspond to a higher number of occurrences. They can be related respectively to the slow and fast separation of vortical disturbances, later to be shed from the main edge structure, see Duguet et al. (2012).

4.1.4. Local expansion rates

When dealing with proper attractors defined over unbounded times, it is common to estimate numerically its dimension. Among the different possible definitions, the Kaplan–Yorke dimension $D_{KY}$ is of interest, because it only requires the values of the leading LEs $\lambda_i$, once ranked in descending order $\lambda_1 > \lambda_2 > \dots > \lambda_r$. It is defined as

$$D_{KY} = j + \frac{S_j}{|\lambda_{j+1}|}, \qquad (4.2)$$

where $S_j$ the cumulative sum

$$S_j = \sum_{i=1}^{j} \lambda_i, \qquad (4.3)$$

and $j$ is the only integer such that $S_j > 0$, but $S_{j+1} < 0$. In the present case the long-time LEs $\lambda_1, \dots$ cannot be computed since the dynamics takes place over finite times. The above definition can, however, be generalised to finite-time problems by considering either the instantaneous or the finite-time exponents (Kuptsov & Kuznetsov 2018). The current analysis is based on the sum $S_j$, rather than on the effective dimension $D_{KY}$ which can be constructed from $S_j$ in (4.3) only if $j$ is large enough. Indeed, with the present value of $r = 8$, there are not enough negative exponents to define $D_{KY}$ according to (4.2). Geometrically, $S_j(t)$ is understood as the instantaneous rate of change of the volume of an infinitesimal state space element defined in the corresponding $j$-dimensional subspace (Kuptsov & Kuznetsov 2018). In figure 9 we show the cumulative sum $S_j(t)$ as a function of time, computed in two different ways. Figure 9(a) has $S_j(t)$ based on the instantaneous
<figure><img src="image_11.png" /><figcaption>Figure 9. (a) Cumulative sum $S_j(t)$ of Re($\lambda_i$), $i = 1, \dots, j$ at different times and for $j = 1, \dots, 8$. (b) Cumulative sum of $\Lambda_{t_0}^{t_0+T}$ at different times.</figcaption></figure>

growth rates Re($\lambda_i$), $i = 1, \dots, j$. Figure 9(b) has $S_{t_0}^{t_0+T}$ based on the FTLEs $\Lambda_{t_0}^{t_0+T}$, which are computed over an entire time interval.

It is observed that, for $j \le 8$ and $t \le 800$, both cumulative sums never become negative. This confirms that the instantaneous and the finite-time Kaplan–Yorke dimensions of the underlying relative attractor are both strictly larger than 8. Interestingly, $S_j$ decreases with $t_0$, up to $t = 550$, for all $j$ values. For the last values of $t_0$ plotted, $S_j$ even eventually decreases with $j$, which suggests that instantaneous eigenvalues with negative growth rate start to contribute to the instantaneous/FTLE spectrum at later times. From a geometric point of view, the fact that $S_j$ stays always positive suggests that the volume of infinitesimal state space elements of the reduced $r$-dimensional space grows with time. This is in contrast with the full $n$-dimensional space where such a volume has to decrease, since the original dynamical system (2.1a,b) is dissipative. In other words, the present reduction, with the choice of $r = 8$, does not incorporate enough dissipative modes, only active modes. Conducting a similar numerical experiment with much larger $r$ is as of today too demanding in terms of memory requirements, at least for the Blasius flow.

### 4.1.5. Summary

The main learnings from the OTD stability analysis restrained to $r = 8$ modes are the following: the dominant edge instability qualifies an outer mode mechanism linked with the wall-normal vorticity of the localised streak. Past the initial 50 time units where the analysis depends on the initialisation of the modes, several mechanisms can be identified from the three peaks in the FTLE spectrum. The dominant instability corresponds to an instability transverse to the edge manifold, while the others correspond to the slow variability of the edge trajectory itself: the dynamics of the perturbations mimics the dynamics inherent to the base flow itself, including the streak phenomenon. The local dimension of the tangent space exceeds the value of $r = 8$. Finally, we observe that the non-normal amplification of disturbances increases when the change of the base flow in time becomes slower and vice versa.
# Instability of the optimal edge trajectory

## 4.2. Modal structures

Beyond global indicators characterising the tangent dynamics, a description of the modal structures in physical space is required. We recall (see § 2.5) that the flow fields visualised correspond to the real part of the vectors $\mathbf{u}_t^1$ defined in (2.9). Two-dimensional visualisations are shown for two different times, namely $t = 280$ and $720$. The modes come in complex conjugate pairs for the considered times, therefore we only display here every other mode among the computed ones. The velocity field of the base flow at these two times, selected along the edge trajectory after the initial transient, is shown in figure 4. It consists of a wiggly finite-length streak flanked with shorter streamwise vortices. At these two times, both snapshots are comparable, the main differences being the longer streamwise extent together (of approximately $500\delta_0^*$) with a spanwise narrower structure of extent $40\delta_0^*$ at the later time. Taking into account the dynamics of the base flow near these two times enriches the description. Near $t = 280$ the formation of streaks by the lift-up mechanism is almost mature (Vavaliaris et al. 2020) and the dynamics relaxes towards quasi-steady motion. By contrast, in the time units following $t = 720$, low and high-speed streak are on the verge of exchanging their spanwise position.

The instantaneous eigenmodes for $t = 280$ are first shown in figures 10–13. The representation, inspired by the experimental figures of Balamurugan & Mandal (2017), is based on a pseudocolour plot of the streamwise velocity perturbation for the reference trajectory, overlapped with lines indicating 40%–100% of the maximum range of the vorticity normal to the planes at $z = -4$, $y = 2.5$ and $x = 325$. The planes are selected to intersect relevant regions of the main structure. We describe now the observed flow structures. The present method as well as the underlying modal decomposition are new in fluid mechanics apart from Babaee & Sapsis (2016). Therefore for pedagogic reasons we chose to display the flow fields of every computed instantaneous eigenmode, omitting the redundant conjugate modes.

For $t = 280$, the spatial structure of each of the 8 leading OTD modes superimposes well with the active part of the main structure, which consists of a sinuous streak of finite length. As a consequence the OTD modes inherit this sinuous structure. Importantly, no spatial symmetry has been imposed neither on the base flow nor on the disturbances modes. This differs from the classical study of Andersson *et al.* (2001) where the base flow has no streamwise dependence. The long-standing question about the symmetries of the leading eigenmodes, namely whether they are symmetric with respect to the plane $z = 0$ (sinuous) or antisymmetric (varicose), becomes irrelevant here. In particular the varicose symmetry, which is consistent with the formation of hairpin vortices, is not characteristic of any of the modes investigated. The classical conclusion of Andersson *et al.* (2001), namely that the sinuous instability of streaks is the most unstable mechanism of paramount importance for streak breakdown, remains valid. Further visualisation of the modes at $t = 280$ highlights the shear layers in the flow, visible in the *xy* plane. The *xz*-plane shows that most of the activity of the mode is located within the active core of the streak and its upstream tail. The *yz*-plane confirms the localisation of the mode on the top shear layer. For all modes, energy is located mostly within the active core or upstream of it. This is in line with the former observation that secondary structures shed downstream of the edge state are not key ingredients of the self-sustained cycle (Duguet *et al.* 2012). Streamwise velocity profiles for the instantaneous eigenmodes are shown in figure 14. They suggest robust localisation close to the edge of the boundary layer. In all panels in figure 14, the *y*-location for the largest amplitude of the streaks is displaced towards larger values with increasing *x*: the head of the streaks characterising the edge trajectory appears tilted upwards. This is again consistent with the description of outer mode instability in Vaughan & Zaki (2011).
M. Beneitez, Y. Duguet, P. Schlatter and D.S. Henningson

<figure><img src="image_12.png" /><figcaption>Figure 10. Snapshots at $t = 280$, $u_1^d$. Two-dimensional grey scale planes of the perturbation velocity of the edge trajectory with arbitrary amplitude, contour lines of $\omega_k k = \{x, y, z\}$, normal to the corresponding plane for each $\{u_1^d\}$, corresponding to 40% to 100% of its maximum value. For each mode, from top to bottom the planes displayed are $z = -4$, $y = 2.5$ and $x = 325$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11. Snapshots at $t = 280$. Mode $u_3^d$. Same as figure 10.</figcaption></figure>
Instability of the optimal edge trajectory

<figure><img src="image_14.png" /><figcaption>Figure 12. Snapshots at $t = 280$. Mode $u_5^l$. Same as figure 10.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13. Snapshots at $t = 280$. Mode $u_8^l$. Same as figure 10.</figcaption></figure>

The relevance of this region is furthermore consistent with the interpretation in Hack & Zaki (2014), where streak instability proceeds via outer modes localised near the edge of the boundary layer. As for the differences between the different modes $u_1^l, \dots, u_8^l$, at $t = 280$ they are not very pronounced yet. Only $u_1^l$ stands out through a less pronounced
<figure><img src="image_16.png" /><figcaption>Figure 14. Streamwise velocity profiles of $u_1^\lambda$ at $t = 280$ corresponding to positions $x = 320$ (solid), $x = 340$ (dashed), $x = 360$ (dot-dashed) for (a) $z = -4$, (b) $z = -2$, (c) $z = 2$, (d) $z = 4$.</figcaption></figure>

tail of streamwise vorticity at the upstream edge. It was checked that perturbing the edge trajectory at $t = 280$ by $u_1^\lambda$, with amplitude $\pm 10^{-4}$, leads either to a turbulent flow or to relaminarisation. This confirms that this eigendirection is transverse to the edge manifold at the considered time.

Most features discussed above are also attested at a later time $t = 720$ as shown in figures 15 and 16, just before streak switching takes place. There are, however, noticeable differences. At $t = 720$, all the instantaneous eigenvalues are still positive, with $\lambda_1$ strictly larger than the other eigenvalues and $\lambda_8$ closer to zero. The leading OTD mode $u_1^\lambda$ is still similar in shape to $u_3^\lambda$ and $u_5^\lambda$, while $u_8^\lambda$ clearly displays a different structure; $u_1^\lambda$ displays strong activity at the edge of the boundary layer, upstream of the active core, strictly above the corresponding shear layer of the base flow (it is most visible on the streamwise velocity component). More noticeable is the fact that the modal structures are lifted towards the edge of the boundary layer, see e.g. the $xy$ plane of figure 16(a) for $u_1^\lambda$. The vortical structures associated with this mode form a larger angle with the wall than the base flow itself. It was again checked that the eigendirection $u_1^\lambda$ is transverse to the edge manifold at the considered time. The structures highlighted in $u_8^\lambda$ are of particular interest. They correspond to the region where a new high-speed streak is in the process of being spanned (see the supplementary material in Vavaliaris et al. (2020) for further evidence). The corresponding OTD mode(s) should hence not only be understood as the manifestation
<figure><img src="image_17.png" /><figcaption>Figure 15. Streamwise velocity profiles of $u_1^d$ at $t = 720$ corresponding to positions $x = 480$ (solid), $x = 500$ (dashed), $x = 520$ (dot-dashed) for (a) $z = -4$, (b) $z = -2$, (c) $z = 2$, (d) $z = 4$.</figcaption></figure>

of an instability of a simple instability-free base flow, instead it can be interpreted as precursor(s) of events that will anyway occur along the edge trajectory. The positive FTLEs associated with the corresponding instantaneous eigenmode are a signature of short-term unpredictability, they quantify the temporal volatility of the streak-switching phenomenon.

Further strengthening the discussion above, figure 20 shows the same snapshots as in figure 4 superimposed now with contours of $\lambda_2$ for the leading OTD mode. It can be seen in both figures 20(a) and 20(b) that the instability mode is mostly localised within the edge structure. The localisation within the active core is even clearer in figure 20(b). Furthermore, figure 21 shows greater localisation in the side where a new high-speed streak is to be generated.

Some elements of this analysis could have been anticipated. The OTD framework, in line with the whole concept of Lyapunov analysis, is a generalisation of modal stability analysis to arbitrarily unsteady base flows. Non-normal features can be captured provided an insightful initialisation of the OTD modes, yet these features are not expected to persist over longer time horizons, e.g. those involved in the evaluation of FTLEs. However, the Orr as well as the lift-up mechanism, which dominate the dynamics at early times, are intrinsically non-normal mechanisms of finite duration. In principle, a large number of eigenvectors is needed to capture transient growth accurately. This explains why so many
<figure><img src="image_18.png" /><figcaption>Figure 16. Snapshots at $t = 720$, $u_1^l$. Two-dimensional grey scale planes of the perturbation velocity of the edge trajectory with arbitrary amplitude, contour lines of $\omega_k$ ($k = \{x, y, z\}$), normal to the corresponding plane for each $\{u_1^l\}$, corresponding to 40% to 100% of its maximum value. For each mode, from top to bottom the planes displayed are $z = -1.6$, $y = 2.5$ and $x = 610$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 17. Snapshots at $t = 720$. Mode $u_3^l$. Same as figure 16.</figcaption></figure>
Instability of the optimal edge trajectory

<figure><img src="image_20.png" /><figcaption>Figure 18. Snapshots at $t = 720$. Mode $u_5^l$. Same as figure 16.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 19. Snapshots at $t = 720$. Mode $u_8^l$. Same as figure 16. Note that the large red structure in the $xy$-plane is located in the region where a new streak will be spanned (Vavaliaris et al. 2020) and the supplementary material therein.</figcaption></figure>

modes possess a similar structure. This trend is aggravated by the fact that for small $r$, the captured non-normality is an estimate of the non-normality of the whole system.

If the description in terms of few OTD modes can seem irrelevant at the earliest times when non-normality dominates, the situation becomes tractable again with small $r$ as soon as the growth of the streaks slows down. The corresponding visualisations for $t = 280$ and
<figure><img src="image_22.png" /><figcaption>Figure 20. Three-dimensional top view of the leading instantaneous eigenmode $\mathbf{u}_1^\lambda$ with arbitrary amplitude, superimposed on the edge trajectory from figure 4; (a) $t = 280$, contours of $u = -5 \times 10^{-2}$ (white), $u = 6 \times 10^{-2}$ (black) for the same snapshot as in figure 4 and $\lambda_2 = -0.29\%$ of its maximum value (green) for the leading OTD mode, (b) $t = 720$, contours of $u = -5 \times 10^{-2}$ (white), $u = 6 \times 10^{-2}$ (black) for the same snapshot as in figure 4 and $\lambda_2 = -0.47\%$ (green) of its maximum value for the leading OTD mode.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 21. Three-dimensional top view of the eighth instantaneous eigenmode $\mathbf{u}_8^\lambda$, with arbitrary amplitude, superimposed on the edge trajectory from figure 4 $t = 720$, contours of $u = -5 \times 10^{-2}$ (white), $u = 6 \times 10^{-2}$ (black) for the same snapshot as in figure 4 and $\lambda_2 = -2.3\%$ (green) of its maximum value.</figcaption></figure>

$t = 720$ are displayed in figures 10–13 and figures 16–19. At this stage the instantaneous eigenvalue distribution as well as FTLE distribution is more comparable to the usual spectrum of edge state solutions, see figure 8: one dominant unstable eigenvalue marking a direction locally transversal to the edge manifold, several weakly positive eigenvalues expressing the chaotic nature of the edge state fluctuations and (not appreciable here because of the small value of $r$) a large set of stable eigenvalues expressing the attraction of the edge state within the edge manifold.

One clear feature from physical space visualisations, regardless of the quantity plotted, is how the localised support of all OTD modes, except here for $\mathbf{u}_8^\lambda$, superimposes exactly with the location of the edge state. This suggests that the present modes, if they contribute to an instability of the edge state, would not make the main coherent edge state spread spatially, at least at the level of the linearised dynamics. As far as the unsteady dynamics restricted to the edge manifold is concerned, this suggests that shift sideways are excluded near $t \approx 280$ whereas they are likely to occur at $t \approx 720$. Such sideways shifts have been reported in most edge states of boundary layer flows (Khapko et al. 2013, 2016; Beneitez et al. 2019), ASBL (Khapko et al. 2016) as well as channel flow (Toh & Itano 2003). As in the present case, the shift phases are usually short and alternate with long shift-free phases. Another robust feature of all localised edge states concerns transition from the edge state to the turbulent state: the transition process consists of two consecutive steps:
first a local intensification of the disturbances within the active core, followed by spatial
spreading (Mellibovsky et al. 2009; Duguet, Willis & Kerswell 2010; Duguet et al. 2013).

The consecutive nature of these two events would suggest that the spreading phase
is nonlinear, while the intensification phase can be understood partially from the linear
instability of the edge state. The fact that the spreading is reflected in the spatial structure
of at least one instantaneous eigenmode $u_8^l$ at the later time $t = 720$, suggests, however,
that spanwise spreading can be partially predicted and described at this time by linear
mechanisms. These new results suggest further study.

**5. Conclusion and outlooks**

We have used the recently developed framework of the OTD modes to study the linearised
dynamics about a segment from a well-defined unsteady base flow. The methodology was
applied to a complex hydrodynamic case at the limit of our computational capabilities, and
yields results in line with the expected physics. It even performs beyond expectations by
revealing new physical phenomena. The physical system under investigation is the Blasius
boundary layer flow. The original trajectory under scrutiny belongs by construction to the
edge manifold delimiting bypass from natural transition. However, the study is restricted
to timespans short enough such that Tollmien–Schlichting waves do not have time to
affect the transition process. This unsteady trajectory is re-interpreted as an unsteady base
flow, whose linear (modal) stability analysis is expected to contain information about the
stability of localised streaks, as observed in instances of bypass transition. This choice
of base flow, due to its three-dimensionality and its unsteady dynamics, represents an
excellent test case for a new stability approach.

Limiting ourselves, as a computational compromise, to a projection basis consisting
of only 8 OTD modes, we have computed the instantaneous eigenvalues along the
unsteady trajectory. The streaky base flow displays a couple of unstable complex conjugate
eigenvalues which dominate the finite-time stability of the trajectory. The remaining
eigenvalues investigated have a positive real part as well, yet with a smaller magnitude.
This is consistent with the expectations for chaotic dynamics within the edge manifold,
although the notion of chaos is usually kept for the infinite-time frameworks.

Numerical evidence suggests that the leading instability mechanism(s) in this study
correspond to an outer mode as described by Vaughan & Zaki (2011) and Hack &
Zaki (2014), even if the corresponding perturbations lack the long wavelength structure
characteristic of streak eigenmodes reported so far (Andersson et al. 2001; Brandt 2014).
We have also analysed the FTLEs along the trajectory by considering several time
horizons. The results confirm the presence of one fast unstable direction vs many slower
state space directions. Moreover, we could confirm that the underlying invariant set has a
finite-time fractal dimension strictly larger than 8.

The leading modal structures obtained from the OTD modes are not trivial to describe,
mainly due to the lack of spatial symmetry of the base flow. The main property exploited
in this study regards the spatial localisation of the modes. Most of the modes computed for
$r = 8$ display, in an instantaneous fashion, the same localisation properties as the original
base flow. The most unstable perturbations display a positive instantaneous growth rate,
and their vortical activity is classically located in the region adjacent to the streaks, where
the total shear is highest (Schmid & Henningson 2001). In particular, the perturbations in
the $xy$-plane are tilted from the wall by an angle larger than the base flow, particularly at
larger times. Some of the modal perturbations extracted also display vortical fluctuations
upstream of the base flow, while one identified mode even displays localisation on the
spanwise side of the base flow (at a later time only). It is suggested that the latter
eigenmode plays an active role as precursor in streak-switching events, the same events that lead the localised edge state to propagate sideways. Downstream fluctuations are, however, absent from the leading instantaneous eigenmodes, suggesting that they are not fundamental to the temporal sustainment of the edge state (Duguet *et al.* 2012).

Although the method originally targets a modal description of the relevant finite-time instabilities, it can also capture non-normal amplification mechanisms (Babaee & Sapsis 2016). In practice, the exact amount of non-normality predicted, as well as the associated energy amplification, are constantly underestimated for finite $r$ compared with the full-dimensional problem, mainly because a larger number of instantaneous eigenmodes would be required to faithfully capture non-normal effects. Nevertheless these results confirm that non-normal effects also play a role in the streak breakdown phenomenon (Schoppa & Hussain 2002; Høpffner *et al.* 2005).

This study highlights the relatively large sensitivity, on short times, of the instantaneous eigenvalues to the initialisation of the OTD modes. It is expected from theoretical arguments (Babaee *et al.* 2017) that FTLEs can be safely computed from the eigenvalues only past a transient time, which is a priori unknown and case dependent. A detailed comparison between two different arbitrary initialisations suggests that, in the present case, only the early times prior to $t \approx 50$ are highly dependent on the choice made for $t=0$ (cf. figure 23. Although this transient can be considered as short relative to the complete transition process, it still represents a clear limitation of the method as far as early times are concerned. At times larger than 50, the instantaneous eigenvalues $\lambda_1, \dots$ evolve qualitatively similarly with time although instantaneously eigenvalues may differ between the two simulations. The corresponding trend is also valid for the numerical abscissa $\sigma$. If the dynamics belonged to an attractor, the time-averaged FTLEs would converge to the LEs, known to be independent of the initialisation Pikovsky & Politi (2016). Although the present case does not revolve around a genuine attractor in state space, the results in figure 23 clearly suggest that the late-time dynamics can be considered as temporally converged. Note that for $r$ large enough the discrepancy between different initialisations is expected to vanish even at finite times, for instantaneous eigenvalues as well as for the numerical abscissa. However, additional modes (and thus larger $r$) also imply a significant increase in computational time.

On the technical level, several points require further discussion and study:

(i) the size of the OTD subspace cannot be determined a priori (Babaee & Sapsis 2016; Kern *et al.* 2021). This is in particular relevant to capture the non-normality along the reference trajectory, where a large number of modes are required. Note that in cases with extensive systems, or ‘weak turbulence’, such as Kuramoto–Sivashinsky (Cvitanović, Davidchack & Siminos 2010) just a few modes are required to entirely describe the most unstable subspace, whereas in pulsating Poiseuille flow more than 70 modes are required to fully describe the non-normal behaviour (Kern *et al.* 2021). However, it has been shown that much lower number of modes $r \approx 6$ could already bring relevant physical insight (Kern *et al.* 2021);

(ii) eigenvalue crossing can make the OTD basis readapt multiple times;

(iii) the modal structures arising from the OTD framework are not associated with a single mode in the sense of classical LSA. The projected OTD modes contain information about several different mechanisms taking place at the same time, in particular for very complex reference trajectories.

To further clarify the potential of the OTD modes in the present complex flow case, we gather our main results in the following list:
Instability of the optimal edge trajectory

(i) first demonstration that the stability analysis of unsteady trajectories is technically possible for a complex large-scale system, without resorting to average LEs or CLVs, even in the case where those might not be available;

(ii) evidence that one unstable mode dominates over all the others at all times, a feature not at all obvious for an aperiodic flow;

(iii) quantification of FTLEs along the edge trajectory, including the early times;

(iv) quantification of the growth of state space volumes as time progresses, showing that at later times the requested number of modes is reduced compared with earlier times;

(v) first quantitative evidence for non-normal effects in an aperiodic flow;

(vi) occurrence of spanwise shifts detected in the higher-order modes at late times;

(vii) evidence that the sinuous symmetry prevails over the streamwise-independent structures throughout all the study. In particular varicose perturbations, known as an alternative way to break streamwise independence and popularised by hairpin vortex studies, appear absent from our study;

(viii) evidence that the new modes found in the present analysis can also be described as outer modes.

Looking ahead, although for intermediate times the OTD modes capture the non-normal features of the underlying linear dynamics, for large times the proposed methodology (edge tracking together with LSA using OTD modes) is essentially a generalisation of modal stability analysis to unsteady cases. Persistent consequences of the non-normality include for instance the finite-time instabilities likely to occur during the Orr mechanism (for $t < 60$) and the lift-up at later times. Both require further extensions of this methodology for a quantitative prediction. The optimal framework proposed by Schmid (2007) is intrinsically non-modal, and it is well suited to the identification of the disturbance most amplified in finite time over an unsteady base flow. The corresponding adjoint-looping algorithm was used successfully by Cossu et al. (2007) in channel flow, except that the reference trajectory chosen was not an edge trajectory but a linear transient. It would be interesting to apply the same methodology on an unsteady edge trajectory and compare the results with the present ones, to see whether one of the methods can predict the finite-time growth of coherent structures not captured by the other technique. Moreover, the possibility of combining these several techniques together is interesting for future developments in stability analysis.

**Acknowledgements.** M.B. and Y.D. would like to thank H. Babaee and S. Kern for discussions about the OTD modes.

**Funding.** Financial support by the Swedish Research Council (VR) grant no. 2016-03541 is gratefully acknowledged. The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) partially funded by the Swedish Research Council through grant agreement no. 2018-05973.

**Declaration of interests.** The authors report no conflict of interest.

Author ORCIDs.

Miguel Beneitez https://orcid.org/0000-0002-4045-7262;

Yohann Duguet https://orcid.org/0000-0001-6258-0475;

Philipp Schlatter https://orcid.org/0000-0001-9627-5903;

Dan S. Henningson https://orcid.org/0000-0001-7864-3071.

Appendix A. Further computational details

This appendix provides further details about the computation of the OTD modes using a pseudospectral approach, and in particular using the SIMSON code (Chevalier *et al.* 2007).
Consider the equations for the evolution of the OTD modes about a trajectory evolved with the Navier–Stokes equations

$$ \dot{u}_i = L_{NS}(u_i) - \langle L_{NS}(u_i), u_i \rangle u_i - \sum_{j=1}^{i-1} [\langle L_{NS}(u_j), u_j \rangle + \langle L_{NS}(u_j), u_i \rangle] u_j, \quad (A1) $$

where $L_{NS}$ denotes the linearised Navier–Stokes operator. Bold letters denote quantities in physical space, which are discretised into $\mathbb{R}^n$ degrees of freedom. The linearised Navier–Stokes functional without external forcing applied to a field $u_i$ reads

$$ L_{NS}(u_i) = -(U_b \cdot \nabla)u_i - (u_i \cdot \nabla)U_b - \nabla p_i + \frac{1}{\text{Re}}\nabla^2 u_i. \quad (A2) $$

The additional constraint to the linearised Navier–Stokes equations is introduced in SIMSON in the form of an explicit forcing at each timestep. Boundary conditions are implemented into the linear part of the solver, while the nonlinear terms are evaluated explicitly. In the present case, evaluating explicitly the inner products involving the linearised Navier–Stokes operator can produce erroneous results if the boundary conditions on the additional forcing term are not applied properly. The term $L_{NS}(u_i)$ needs to contain the boundary conditions corresponding to the linearised operator to provide a correct forcing term and $L_r$ in the computations. In particular, it is necessary to apply Neumann boundary conditions on the free stream and Dirichlet boundary conditions at the wall

$$ u_i(y=0) = \frac{\partial u_i}{\partial y}(y=L_y) = 0, \quad (A3) $$

when recovering the wall-normal velocity from the fourth-order equation arising from the velocity–vorticity formulation (Chevalier et al. 2007). This differs from the main body of the implementation in SIMSON since the correct boundary conditions need to be applied in the explicit term involving $L_{NS}(u_i)$ as well as the implicit part of the solver.

The OTD modes converge exponentially fast to the most unstable directions of the Cauchy–Green tensor (Babaee et al. 2017), and after a long time only depend on the point of the trajectory where they are computed (Blanchard & Sapsis 2019). However, the OTD modes depend on their initialisation (Babaee & Sapsis 2016; Kern et al. 2021). It has been observed that there is not a universal time for which the OTD subspace is converged to the most unstable directions. Nevertheless, relevant physical features may be observed from early times (Babaee & Sapsis 2016).

## Appendix B. Initial conditions for the OTD modes

An additional point to consider is that, if one of the directions not part of the basis becomes unstable enough, the basis will need to re-adapt. This is due to (2.5) being evolved continuously, whereas the introduction of a different vector in the most unstable subspace occurs discontinuously (Babaee & Sapsis 2016).

The choice of the initial condition for the OTD modes plays therefore a crucial role for the OTD framework. Since our reference trajectory consists of several finite-time events of interest, our goal is to choose initial conditions which adapt as quickly as possible to the most unstable dynamics. We therefore chose initial conditions which are physically relevant to excite instability mechanisms on the edge trajectory. The initial condition for the first mode is exactly the perturbation to the Blasius boundary layer associated with the edge state. This represents infinitesimal perturbations of the same
Instability of the optimal edge trajectory

<figure><img src="image_24.png" /><figcaption>Figure 22. Snapshots at $t = 120$. Contours of the wall-normal velocity $V$ for $u_1^t$ obtained from two different sets of initial conditions. (a) Subset of initial conditions used in the current work. (b) Alternative initial conditions described in Appendix B. Note that according to (2.9), the {$u^t$} are not normalised.</figcaption></figure>

shape as the edge trajectory. The initial conditions for modes 2–8 correspond to pairs of counter-rotating vortices with different spatial extensions. The counter-rotating vortices are also rotated about the y axis to remove any symmetric constraint. This set of initial conditions is not orthogonal by construction and therefore a Gram-Schmidt algorithm is performed before initialising the OTD computations. The results reported in the body of the paper correspond to these initial conditions.

To check the robustness of the results, alternative sets of initial conditions have been tested using $r = 4$: (i) the 4 leading modes from the results in the main body of the paper and (ii) random noise. Using $r = 4$ modes only appeared sufficient to illustrate the main aspects of the subsequent checks.

A comparison for the wall-normal component of the leading projected OTD mode at $t = 120$, obtained using the two different sets of initial conditions can be seen in figure 22. The figure shows an agreement about the general physical features of the perturbation, i.e. the high-speed streak flanked by two low-speed streaks is present in both cases. However, no exact match is observed. The convergence to a unique set of OTD modes is expected to be exponentially fast (Babaee & Sapsis 2016; Babaee et al. 2017; Blanchard & Sapsis 2019), but the explicit times are strongly case dependent.
<figure><img src="image_25.png" /><figcaption>Figure 23. Instantaneous eigenvalues for $r = 4$ corresponding to two different initialisations: (a) edge state perturbation over the Blasius boundary layer and counter rotating vortices, (b) random noise.</figcaption></figure>

The most unstable instantaneous eigenvalues are shown in figure 23. It can be observed that, in the case of the random noise, the initial peak is lost. It is reasonable to assume that the unsteady base flow changes too fast during the initial times while the OTD subspace has not had enough time to adapt. On the other hand, the second peak identified at $t \sim 80$ is well identified with both sets of initial conditions.

We should consider random noise as the worst choice of initial conditions, since it is entirely agnostic to the underlying reference trajectory. The results presented above further strengthen the importance of the choice of initial conditions. They indicate that, although the OTD approach is robust at large enough times, it remains dependent on the initialisation for times earlier than $t \approx 100$.

REFERENCES

ALFREDSSON, P.H. & MATSUBARA, M. 1996 Streaky structures in transition. *Proc. Trans. Boundary Layer Aeronaut.* **46**, 373–386.

ALIZARD, F. 2015 Linear stability of optimal streaks in the log-layer of turbulent channel flows. *Phys. Fluids* **27** (10), 105103.

ANDERSSON, P., BRANDT, L., BOTTARO, A. & HENNINGSON, D.S. 2001 On the breakdown of boundary layer streaks. *J. Fluid Mech.* **428**, 29.

ASAI, M., MINAGAWA, M. & NISHIOKA, M. 2002 The instability and breakdown of a near-wall low-speed streak. *J. Fluid Mech.* **455**, 289.

BABAEE, H., FARAZMAND, M., HALLER, G. & SAPSIS, T.P. 2017 Reduced-order description of transient instabilities and computation of finite-time Lyapunov exponents. *Chaos* **27** (6), 063103.

BABAEE, H. & SAPSIS, T.P. 2016 A minimization principle for the description of modes associated with finite-time instabilities. *Proc. R. Soc. Lond. A* **472** (2186), 20150779.

BAKCHINOV, A.A., GREK, G.R., KLINGMANN, B.G.B. & KOZLOV, V.V. 1995 Transition experiments in a boundary layer with embedded streamwise vortices. *Phys. Fluids* **7** (4), 820–832.

BALAMURUGAN, G. & MANDAL, A.C. 2017 Experiments on localized secondary instability in bypass boundary layer transition. *J. Fluid Mech.* **817**, 217.

BENEITEZ, M., DUGUET, Y. & HENNINGSON, D.S. 2020a Modeling the collapse of the edge when two transition routes compete. *Phys. Rev. E* **102** (5), 053108.

BENEITEZ, M., DUGUET, Y., SCHLATTER, P. & HENNINGSON, D.S. 2020b Edge manifold as a lagrangian coherent structure in a high-dimensional state space. *Phys. Rev. Res.* **2** (3), 033258.

BENEITEZ, M., DUGUET, Y., SCHLATTER, P. & HENNINGSON, D.S. 2019 Edge tracking in spatially developing boundary layer flows. *J. Fluid Mech.* **881**, 164–181.

BIAU, D. 2012 Laminar-turbulent separatrix in a boundary layer flow. *Phys. Fluids* **24** (3), 034107.
BLANCHARD, A., MOWLAVI, S. & SAPSIS, T.P. 2019 Control of linear instabilities by dynamically consistent order reduction on optimally time-dependent modes. *Nonlinear Dyn.* **95** (4), 2745–2764.

BLANCHARD, A. & SAPSIS, T.P. 2019 Analytical description of optimally time-dependent modes for reduced-order modeling of transient instabilities. *SIAM J. Appl. Dyn. Syst.* **18** (2), 1143–1162.

BRANDT, L. 2007 Numerical studies of the instability and breakdown of a boundary-layer low-speed streak. *Eur. J. Mech. (B/Fluids)* **26** (1), 64–82.

BRANDT, L. 2014 The lift-up effect: the linear mechanism behind transition and turbulence in shear flows. *Eur. J. Mech. (B/Fluids)* **47**, 80–96.

BRANDT, L., COSSU, C., CHOMAZ, J-M., HUERRE, P. & HENNINGSON, D.S. 2003 On the convectively unstable nature of optimal streaks in boundary layers. *J. Fluid Mech.* **485**, 221–242.

BRANDT, L. & HENNINGSON, D.S. 2002 Transition of streamwise streaks in zero-pressure-gradient boundary layers. *J. Fluid Mech.* **472**, 229–261.

BRANDT, L., SCHLATTER, P. & HENNINGSON, D.S. 2004 Transition in boundary layers subject to free-stream turbulence. *J. Fluid Mech.* **517**, 167.

CASSINELLI, A., DE GIOVANETTI, M. & HWANG, Y. 2017 Streak instability in near-wall turbulence revisited. *J. Turbul.* **18** (5), 443–464.

CHERUBINI, S., DE PALMA, P., ROBINET, J.-CH. & BOTTARO, A. 2011a The minimal seed of turbulent transition in the boundary layer. *J. Fluid Mech.* **689**, 221–253.

CHERUBINI, S., DE PALMA, P., ROBINET, J.-CH. & BOTTARO, A. 2011b Edge states in a boundary layer. *Phys. Fluids* **23** (5), 051705.

CHEVALIER, M., SCHLATTER, P., LUNDBLADH, A. & HENNINGSON, D.S. 2007 Simson – a pseudo-spectral solver for incompressible boundary layer flow. *Tech. Rep. TRITA-MEK* 2007:07.

COSSU, C., CHEVALIER, M. & HENNINGSON, D.S. 2007 Optimal secondary energy growth in a plane channel flow. *Phys. Fluids* **19** (5), 058107.

CVITANOVIĆ, P., DAVIDCHACK, R.L. & SIMINOS, E. 2010 On the state space geometry of the Kuramoto–Sivashinsky flow in a periodic domain. *SIAM J. Appl. Dyn. Syst.* **9** (1), 1–33.

DRAZIN, P.G. & REID, W.H. 2004 *Hydrodynamic Stability*, 2nd edn. Cambridge University Press.

DUGUET, Y., MONOKROUSOS, A., BRANDT, L. & HENNINGSON, D.S. 2013 Minimal transition thresholds in plane Couette flow. *Phys. Fluids* **25** (8), 084103.

DUGUET, Y., SCHLATTER, P., HENNINGSON, D.S. & ECKHARDT, B. 2012 Self-sustained localized structures in a boundary-layer flow. *Phys. Rev. Lett.* **108** (4), 044501.

DUGUET, Y., WILLIS, A.P. & KERSWELL, R.R. 2010 Slug genesis in cylindrical pipe flow. *J. Fluid Mech.* **663**, 180.

DUGUET, Y., WILLIS, A.P. & KERSWELL, R.R. 2008 Transition in pipe flow: the saddle structure on the boundary of turbulence. *J. Fluid Mech.* **613**, 255–274.

EMBREE, M. & TREFETHEN, L.N. 2005 Spectra and pseudospectra: the behavior of nonnormal matrices and operators.

FARAZMAND, M. & SAPSIS, T.P. 2016 Dynamical indicators for the prediction of bursting phenomena in high-dimensional systems. *Phys. Rev. E* **94** (3), 032212.

FOURES, D.P.G., CAULFIELD, C.P. & SCHMID, P.J. 2013 Localization of flow structures using infty-norm optimization. *J. Fluid Mech.* **729**, 672–701.

GIBSON, J.F., HALCROW, J. & CVITANOVIĆ, P. 2008 Visualizing the geometry of state space in plane Couette flow. *J. Fluid Mech.* **611**, 107–130.

GINELLI, F., POGGI, P., TURCHI, A., CHATÉ, H., LIVI, R. & POLITI, A. 2007 Characterizing dynamics with covariant Lyapunov vectors. *Phys. Rev. Lett.* **99** (13), 130601.

HACK, M.J.P. & ZAKI, T.A. 2014 Streak instabilities in boundary layers beneath free-stream turbulence. *J. Fluid Mech.* **741**, 280.

HALLER, G. 2015 Lagrangian coherent structures. *Annu. Rev. Fluid Mech.* **47** (1), 137–162.

HAMILTON, J.M., KIM, J. & WALEFFE, F. 1995 Regeneration mechanisms of near-wall turbulence structures. *J. Fluid Mech.* **287** (1), 317–348.

HENNINGSON, D.S., LUNDBLADH, A. & JOHANSSON, A.V. 1993 A mechanism for bypass transition from localized disturbances in wall-bounded shear flows. *J. Fluid Mech.* **250**, 169–207.

HÖPFNER, J., BRANDT, L. & HENNINGSON, D.S. 2005 Transient growth on boundary layer streaks. *J. Fluid Mech.* **537** (1), 91–100.

ITANO, T. & TOH, S. 2001 The dynamics of bursting process in wall turbulence. *J. Phys. Soc. Japan* **70** (3), 703–716.

JACOBS, R.G. & DURBIN, P.A. 2001 Simulations of bypass transition. *J. Fluid Mech.* **428**, 185–212.

JEONG, J. & HUSSAIN, F. 1995 On the identification of a vortex. *J. Fluid Mech.* **285**, 69–94.
KAWAHARA, G., JIMÉNEZ, J., UHLMANN, M. & PINELLI, A. 2003 Linear instability of a corrugated vortex sheet-a model for streak instability. *J. Fluid Mech.* **483**, 315.

KERN, J.S., BENEITEZ, M., HANIFI, A. & HENNINGSON, D.S. 2021 Transient linear stability of pulsating poiseuille flow using optimally time-dependent modes. *J. Fluid Mech.* **927**, A6.

KERSWELL, R.R. 2018 Nonlinear nonmodal stability theory. *Annu. Rev. Fluid Mech.* **50** (1), 319–345.

KHAPKO, T., DUGUET, Y., KREILOS, T., SCHLATTER, P., ECKHARDT, B. & HENNINGSON, D.S. 2014 Complexity of localised coherent structures in a boundary-layer flow. *Eur. Phys. J. E* **37** (4), 32.

KHAPKO, T., KREILOS, T., SCHLATTER, P., DUGUET, Y., ECKHARDT, B. & HENNINGSON, D.S. 2013 Localized edge states in the asymptotic suction boundary layer. *J. Fluid Mech.* **717**.

KHAPKO, T., KREILOS, T., SCHLATTER, P., DUGUET, Y., ECKHARDT, B. & HENNINGSON, D.S. 2016 Edge states as mediators of bypass transition in boundary-layer flows. *J. Fluid Mech.* **801**.

KLEBANOFF, P.S., TIDSTROM, K.D. & SARGENT, L.M. 1962 The three-dimensional nature of boundary-layer instability. *J. Fluid Mech.* **12** (1), 1–34.

KUPTSOV, P.V. & PARLITZ, U. 2012 Theory and computation of covariant Lyapunov vectors. *J. Nonlinear Sci.* **22** (5), 727–762.

KUPTSOV, P.V. & KUZNETSOV, S.P. 2018 Lyapunov analysis of strange pseudohyperbolic attractors: angles between tangent subspaces, local volume expansion and contraction. *Regular Chaotic Dyn.* **23** (7), 908–932.

LANDAHL, M.T. 1980 A note on an algebraic instability of inviscid parallel shear flows. *J. Fluid Mech.* **98** (2), 243–251.

MATSUBARA, M. & ALFREDSSON, P.H. 2001 Disturbance growth in boundary layers subjected to free-stream turbulence. *J. Fluid Mech.* **430**, 149.

MELLIBOVSKY, F., MESEGUER, A., SCHNEIDER, T.M. & ECKHARDT, B. 2009 Transition in localized pipe flow turbulence. *Phys. Rev. Lett.* **103** (5), 054502.

MORKOVIN, M.V. 1969 On the many faces of transition. In *Viscous Drag Reduction*, pp. 1–31. Springer.

NORDSTRÖM, J., NORDIN, N. & HENNINGSON, D. 1999 The fringe region technique and the fourier method used in the direct numerical simulation of spatially evolving viscous flows. *SIAM J. Sci. Comput.* **20** (4), 1365–1393.

PIKOVSKY, A. & POLITI, A. 2016 *Lyapunov Exponents: A Tool to Explore Complex Dynamics*. Cambridge University Press.

PRINGLE, C.C.T. & KERSWELL, R.R. 2010 Using nonlinear transient growth to construct the minimal seed for shear flow turbulence. *Phys. Rev. Lett.* **105** (15), 154502.

RABIN, S.M.E., CAULFIELD, C.P. & KERSWELL, R.R. 2012 Triggering turbulence efficiently in plane Couette flow. *J. Fluid Mech.* **712**, 244.

RINALDI, E., CANTON, J. & SCHLATTER, P. 2019 The vanishing of strong turbulent fronts in bent pipes. *J. Fluid Mech.* **866**, 487–502.

SAPSIS, T.P. 2018 New perspectives for the prediction and statistical quantification of extreme events in high-dimensional dynamical systems. *Phil. Trans. R. Soc. Lond.* A **376** (2127), 20170133.

SCHLATTER, P., BRANDT, L., DE LANGE, H.C. & HENNINGSON, D.S. 2008 On streak breakdown in bypass transition. *Phys. Fluids* **20** (10), 101505.

SCHMID, P.J. 2007 Nonmodal stability theory. *Annu. Rev. Fluid Mech.* **39**, 129–162.

SCHMID, P.J. & HENNINGSON, D.S. 2001 *Stability and Transition in Shear Flows*. Springer.

SCHOPPA, W. & HUSSAIN, F. 2002 Coherent structure generation in near-wall turbulence. *J. Fluid Mech.* **453**, 57.

SHIMADA, I. & NAGASHIMA, T. 1979 A numerical approach to ergodic problem of dissipative dynamical systems. *Prog. Theor. Phys.* **61** (6), 1605–1616.

SKUFCA, J.D., YORKE, J.Y. & ECKHARDT, B. 2006 Edge of chaos in a parallel shear flow. *Phys. Rev. Lett.* **96** (17).

TOH, S. & ITANO, T. 2003 A periodic-like solution in channel flow. *J. Fluid Mech.* **481**, 67–76.

VAUGHAN, N.J. & ZAKI, T.A. 2011 Stability of zero-pressure-gradient boundary layer distorted by unsteady Klebanoff streaks. *J. Fluid Mech.* **681**, 116.

VAVALIARIS, C., BENEITEZ, M. & HENNINGSON, D.S. 2020 Optimal perturbations and transition energy thresholds in boundary layer shear flows. *Phys. Rev. Fluids* **5** (6), 062401.

VIANA, M. 2014 *Lectures on Lyapunov Exponents*, vol. 145. Cambridge University Press.

WALEFFE, F. 1997 On a self-sustaining process in shear flows. *Phys. Fluids* **9** (4), 883–900.