Instability of the optimal edge trajectory in the Blasius boundary layer

Abstract

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.

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\gg 1$ (see e.g. Gibson, Halcrow & Cvitanović 2008), we consider $\mathbb {R}^{n}$ as the original high-dimensional space of reference. We suppose a non-autonomous dynamical system defined over a time interval $[t_0,t_1)$

(2.1a,b)\begin{equation} \frac{{\rm d}{\boldsymbol{Q}}}{{\rm d}t}={\boldsymbol{f}}({\boldsymbol{Q}},t),\quad {\boldsymbol{Q}}(t_0)={\boldsymbol{Q}}_0, \end{equation}

where ${\boldsymbol {Q}}_0,{\boldsymbol {Q}}\in \mathbb {R}^{n}$, and ${\boldsymbol {f}}:\mathbb {R}^{n} \rightarrow \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 ${\boldsymbol {Q}}_0$, we define the solution to (2.1a,b), namely $\bar {{\boldsymbol {Q}}}:(t_0:t_1)\rightarrow \mathbb {R}^{n}$ as the base flow whose stability we will now determine.

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

(2.2a,b)\begin{equation} \frac{{\rm d}{\boldsymbol{q}}}{{\rm d}t}={\boldsymbol{L}}(\bar{{\boldsymbol{Q}}},t){\boldsymbol{q}}, \quad {\boldsymbol{L}}(\bar{{\boldsymbol{Q}}},t)=\boldsymbol{\nabla}_{\boldsymbol{Q}} {\boldsymbol{f}}(\bar{{\boldsymbol{Q}}},t), \end{equation}

where ${\boldsymbol {\nabla }}_{\boldsymbol {Q}} {\boldsymbol {f}}(\bar {{\boldsymbol {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 $\bar {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

(2.3)\begin{equation} \mathcal{J}=\lim_{\delta t\to 0}\frac{1}{(\delta t)^2}\sum_{i=1}^{n}||{\boldsymbol{w}}_i(t+\delta t)-({\boldsymbol{\nabla}} {\boldsymbol{F}}_{t}^{t+\delta t}){\boldsymbol{w}}_i(t)||^2, \end{equation}

where ${\boldsymbol {F}}_{t}^{t+\delta t}:\mathbb {R}^{n} \rightarrow \mathbb {R}^{n}$ is the infinitesimal forward propagator associated with (2.1a,b). The Jacobian matrix ${\boldsymbol {\nabla }} {\boldsymbol {F}}_{t}^{t+\delta t}$ maps a vector of the tangent space $\boldsymbol {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 ${\boldsymbol {u}}_i$, $i=1,\ldots,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

(2.4)\begin{equation} \langle {\boldsymbol{u}}_i(t), {\boldsymbol{u}}_j(t)\rangle=\delta_{ij},\quad i,j=1,\ldots,r, \ r\ll n, \end{equation}

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 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).

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.

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

(2.5)\begin{equation} \frac{{{\rm d}\boldsymbol{u}}_i}{{\rm d}t} = {\boldsymbol{L}}(t) {\boldsymbol{u}}_i-\sum_{j=1}^r [\langle {\boldsymbol{L}}(t) {\boldsymbol{u}}_i, {\boldsymbol{u}}_j\rangle-{\boldsymbol{A}}_{ij}(t)]{\boldsymbol{u}}_j,\quad i=1,\ldots,r . \end{equation}

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

(2.6)\begin{equation} {\boldsymbol{A}}_{ij} = \left\{\!\!\begin{array}{ll} -\langle {\boldsymbol{L}} {\boldsymbol{u}}_j, {\boldsymbol{u}}_i \rangle, & j< i\\ 0, & i=j \\ \langle \boldsymbol{L} {\boldsymbol{u}}_i, {\boldsymbol{u}}_j \rangle, & j>i, \end{array}\right. \end{equation}

which is also considered here. An arbitrary choice of $\boldsymbol {A}$ would lead to a fully coupled system so that each ${\boldsymbol {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{\rm th}$ mode depends only on the modes from $1$ to $i$, making the OTD formulation hierarchical. The resulting evolution equation for each OTD mode becomes

(2.7)\begin{equation} \frac{{\rm d}{\boldsymbol{u}}_i}{{\rm d}t} = {\boldsymbol{L}}(t) {\boldsymbol{u}}_i-\langle {\boldsymbol{L}}(t) {\boldsymbol{u}}_i, {\boldsymbol{u}}_i \rangle {\boldsymbol{u}}_i-\sum_{j=1}^{i-1} [\langle {\boldsymbol{L}}(t) {\boldsymbol{u}}_i, {\boldsymbol{u}}_j\rangle + \langle {\boldsymbol{L}}(t) {\boldsymbol{u}}_j, {\boldsymbol{u}}_i\rangle ]{\boldsymbol{u}}_j. \end{equation}

The nonlinear system of $r$ (2.7) can be evolved forward in time together with (2.1a,b), from which the matrix ${\boldsymbol {L}}$ can be evaluated at all times. Equation (2.1a,b) remains, however, independent of the evolution of each ${\boldsymbol {u}}_i$. This results in an $(r+1)\times 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 ${\boldsymbol {L}}_r$ defined by projecting the high-dimensional operator ${\boldsymbol {L}}$ onto the OTD directions

(2.8)\begin{equation} \boldsymbol{L}_{r_{ij}}(t) = \langle {\boldsymbol{u}}_i,{\boldsymbol{L}}(t){\boldsymbol{u}}_j\rangle\quad i,j=1,\dots,r. \end{equation}

In particular, all the instantaneous stability indicators defined in the next subsection will be derived from algebraic properties of the $r \times r$ matrix ${\boldsymbol {L}}_r$ evaluated at the relevant times.

For a time-independent linearised operator $L$ the space spanned by the modes $\{u_i\}_{i=1}^{r}$ 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 ${\boldsymbol {u}}_i$ 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 ${\boldsymbol {u}}_i$ modes and the reduced operator. At every time, $L_r(t)$ can be diagonalised as ${\boldsymbol {L}}_r={\boldsymbol {E}}^{\lambda }\boldsymbol {\varLambda }_r(\boldsymbol {E}^{\lambda })^{-1}$ with ${\boldsymbol {E}}^{\lambda }$ and ${\boldsymbol {\varLambda }}_r=\text {diag}(\lambda _1(t),\ldots,\lambda _r(t))$. The new modes ${\boldsymbol {u}}_i^{\lambda }, i=1,\ldots,r$ are defined (using the summation convention) by

(2.9)\begin{equation} {\boldsymbol{u}}_i^{\lambda}(t)=\boldsymbol{E}^{\lambda}_{ij}(t){\boldsymbol{u}}_j(t). \end{equation}

Unlike the ${\boldsymbol {u}}_i$ modes, the ${\boldsymbol {u}}_i^{\lambda }$ 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 ${\boldsymbol {u}}_i$ modes are real valued, the modes ${\boldsymbol {u}}_i^{\lambda }$ 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 $\lambda _1(t),\ldots,\lambda _r(t)$, $i=1,\ldots,r$ are labelled instantaneous eigenvalues and are complex valued.

Another key scalar quantity is the instantaneous numerical abscissa $\sigma (t)$, a positive number, defined as the largest eigenvalue of the symmetrised reduced operator $({\boldsymbol {L}}+{\boldsymbol {L}}^T)/2$ (Embree & Trefethen 2005). Here, $\sigma$ corresponds to the largest possible growth rate at a given time, due to both normal and non-normal effects combined together. Whenever ${\boldsymbol {L}}$ is non-normal, $\sigma$ is strictly larger than the real part of all $\lambda _i$ values.

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

(2.10)\begin{equation} g_r(t):=\min\limits_i|\sigma_r-{\rm Re}(\lambda_i)|=\sigma_r-{\rm Re}(\lambda_1), \end{equation}

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 ${\boldsymbol {F}}_{t_0}^{t}$, the Cauchy–Green tensor $\boldsymbol {C}_{t_0}^{t}$ is defined as

(2.11)\begin{equation} {\boldsymbol{C}}_{t_0}^{t}({\boldsymbol{q}}_0) :=[{\boldsymbol{\nabla}} {\boldsymbol{F}}_{t_0}^{t}({\boldsymbol{q}}_0)]^{\rm T} [{\boldsymbol{\nabla}} {\boldsymbol{F}}_{t_0}^{t}({\boldsymbol{q}}_0)]. \end{equation}

The associated FTLEs are defined directly from the eigenvalues $\gamma _1>\gamma _2>\cdots >\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

(2.12)\begin{equation} (\varLambda_{t_0}^{t_0+T})_i=\frac{1}{T}\log{\sqrt{\gamma_i}},\quad i=1,\ldots,n. \end{equation}

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 $(\varLambda _{t_0}^{t_0+T})_i,~i=1,\ldots,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 $(\varLambda _{t_{0}}^{t_{0}+T})_i$, $i=1,\ldots,r$ are evaluated simply by averaging over time the diagonal elements of ${\boldsymbol {L}}_r$ (Blanchard & Sapsis 2019). They are expressed as

(2.13)\begin{equation} (\varLambda_{t_{0}}^{t_{0}+T})_i \approx \frac{1}{T}\int_{t_{0}}^{t_{0}+T}\langle {\boldsymbol{u}}_i(\tau),{\boldsymbol{L}}_r(\tau){\boldsymbol{u}}_i(\tau)\rangle \,{\rm d}\tau,\quad i=1,\ldots,r. \end{equation}

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, $\boldsymbol {v}$ is the total velocity field, $\boldsymbol {v}_{B}=(u_B,v_B,0)$ that of the steady Blasius solution, then $\boldsymbol {u}:=\boldsymbol {v}-\boldsymbol {v}_{B}=(u,v,w)$ is the perturbation velocity field. All quantities are made non-dimensional using the free-stream velocity $U_ {\infty }$ and the boundary layer thickness $\delta ^{*}(x):=\int _{0}^{\infty } (1-u_B(x)/U_{\infty })\,{{\rm d}y}$ of the undisturbed (steady) Blasius flow. A local Reynolds number can be defined as $Re_{\delta ^{*}}(x):=U_{\infty }\delta ^{*}/\nu$, with $\nu$ the kinematic viscosity of the fluid. The value of $Re_{\delta ^{*}_{0}}=Re_{\delta ^{*}}(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

(3.1)\begin{equation} u=v=w=0 , \end{equation}

and at the upper domain boundary ($y=L_y$) of Neumann type to allow for a natural growth of the boundary layer

(3.2)\begin{equation} \frac{\partial u}{\partial y}=\frac{\partial v}{\partial y}=\frac{\partial w}{\partial y}=0. \end{equation}

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 $\boldsymbol {F}(t,x,y,z)$ of the form

(3.3)\begin{equation} \boldsymbol{F}=\gamma (x)(\mathcal{U}(x,y,z)-\boldsymbol{v}(t,x,y,z)), \end{equation}

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

(3.4)\begin{equation} \mathcal{U}_x=U(x,y,z)+[U(x+x_L,y,z)-U(x,y,z)]S\left(\frac{x-x_{blend}}{\varDelta_{blend}} \right), \end{equation}

where $S(x_{blend},\varDelta _{blend})$ 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 $\mathcal{U}$ is obtained via the continuity equation. In the present work the fringe length is $\varDelta _{blend}=600$, $x_L$=2500 and $\gamma _{max}=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 $\Delta t=0.2$ in terms of $U_{\infty }$ and $\delta _0^{*}$. The velocity field is expanded along $N_x$ Fourier modes in the streamwise direction $x$ and $N_z$ modes in the spanwise direction $z$, $N_y$ 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,\ldots,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 $\varOmega$ 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

(3.5)\begin{equation} \langle {\boldsymbol{u}},{\boldsymbol{u'}} \rangle = \int_{\varOmega} (uu'+vv'+ww')\,{\rm d}\varOmega, \end{equation}

where $\boldsymbol {u}=(u,v,w)$ and $\boldsymbol {u'}=(u',v',w')$ are any two flow fields with finite $L^2$ norm, and ${\rm d}\varOmega ={\rm d}\kern0.06em x\,{\rm d}y\,{\rm d}z$ is the usual infinitesimal integration element over the numerical domain $\varOmega$.

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_{\delta _0^{*}}=240.458$. That value of $Re_{\delta }(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 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\leqslant 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)

(3.6)$$\begin{gather} \varOmega_x=(\delta_0^*/\delta)^{{1}/{2}}\left(\frac{1}{V}\int_V|\omega_x|^2\,\mathrm{d}v\right)^{{1}/{2}}, \end{gather}$$

(3.7)$$\begin{gather}\varOmega_y=(\delta_0^*/\delta)^{{1}/{2}}\left(\frac{1}{V}\int_V|\omega_y|^2\,\mathrm{d}v\right)^{{1}/{2}}, \end{gather}$$

(3.8)$$\begin{gather}W=(\delta_0^*/\delta)^{{3}/{2}}\left(\frac{1}{V}\int_V|w|^2\,\mathrm{d}v\right)^{{1}/{2}}. \end{gather}$$

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]$ 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.

Figure 2. State portrait of the optimal edge trajectory in the Blasius boundary layer using the variables $\varOmega _x$, $\varOmega _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.

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 ${\sim }22$ length units.

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$.

It is useful to recall the main features of the unsteady base flow reported by Vavaliaris et al. (2020). For early times $t \leqslant 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 \leqslant t \leqslant 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 (Schmid & Henningson 2001). For $t \geqslant 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 \geqslant 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,\ldots,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) - {\rm 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.

Figure 5. Real part of the instantaneous eigenvalues ${\rm Re}(\lambda _i)$, $i=1,\ldots,8$ (solid lines) and instantaneous numerical abscissa $\sigma$ (dashed red) of the reduced operator ${\boldsymbol {L}}_r$, plotted vs time. Here, $\lambda _1$ is highlighted (orange line). Vertical dashed line indicates $t=50$.

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: ${\rm d}\varOmega _x/{\rm d}t$, ${\rm d}\varOmega _y/{\rm d}t$ and ${\rm d}W/{\rm d}t$ as defined in (4.1) (red). The figure starts at $t=50$.

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 time while the others are all negative. In a second phase starting at $t \approx$ 100, ${\rm 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 \geqslant 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 ={\rm Re}(\lambda _1)={\rm 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 ${\rm 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 \leqslant 50$, may be wrongly attributed to a lack of non-normal potential of $\boldsymbol {L}(t)$. To start with, this is a property of the instantaneous reduced operator $\boldsymbol {L}_r(t)$ computed for a given value of $r$, not necessarily of the full operator $\boldsymbol {L}(t)$. The reverse is yet true: non-normal features of the reduced-order operator $\boldsymbol {L}_r(t)$ carry over to $\boldsymbol {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\leqslant 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 $\varOmega _x$, $\varOmega _y$ and $W$ used in figure 2. This quantity is defined as $\xi$ via

(4.1)\begin{equation} \xi(t) = \sqrt{\left(\frac{{\rm d}\varOmega_x}{{\rm d}t}\right)^2+\left(\frac{{\rm d}\varOmega_y}{{\rm d}t}\right)^2+C\left(\frac{{\rm d}W}{{\rm d}t}\right)^2}, \end{equation}

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\geqslant 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 ${\rm Re}(\lambda _{1,2})$ plotted vs $\varOmega _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 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 $\varOmega _x \geqslant 0.05$. Interpreting $\varOmega _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.

Figure 7. (a) Green lines: ${\rm Re}(\lambda _{1,2})$ vs $\varOmega _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 ${\rm Re}(\lambda _{1,2})$ vs $||w||^2$, the volume average spanwise velocity.

Figure 8. Histograms of all $r=8$ reduced-order FTLEs $(\varLambda _{t_{0}}^{t_{0}+T})_i$ ($i=1,\ldots,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]$.

4.1.3. Finite-time Lyapunov exponents

Figure 8(a) shows distributions of FTLEs $(\varLambda _{t_{0}}^{t_{0}+T})_i$ ($i=1,\ldots,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 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 ($\geqslant$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>\cdots >\lambda _r$. It is defined as

(4.2)\begin{equation} D_{KY}=j+\frac{S_j}{|\lambda_{j+1}|}, \end{equation}

where $S_j$ the cumulative sum

(4.3)\begin{equation} S_j=\sum_{i=1}^{j}\lambda_{i}, \end{equation}

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,\ldots$ 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 growth rates ${\rm Re}(\lambda _i)$, $i=1,\ldots,j$. Figure 9(b) has $S_{t_0}^{t_0+T}$ based on the FTLEs $\varLambda _{t_{0}}^{t_{0}+T}$, which are computed over an entire time interval.

Figure 9. (a) Cumulative sum $S_j(t)$ of ${\rm Re}(\lambda _i)$, $i=1,\ldots,j$ at different times and for $j=1,\ldots,8$. (b) Cumulative sum of $\varLambda _{t_{0}}^{t_{0}+T}$ at different times.

It is observed that, for $j\leqslant 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.

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 ${\boldsymbol {u}}_i^{\lambda }$ 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.

Figure 10. Snapshots at $t=280$, $u^{\lambda }_1$. 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_i^{\lambda }\}$, 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$.

Figure 11. Snapshots at $t=280$. Mode $u^{\lambda }_3$. Same as figure 10.

Figure 12. Snapshots at $t=280$. Mode $u^{\lambda }_5$. Same as figure 10.

Figure 13. Snapshots at $t=280$. Mode $u^{\lambda }_8$. Same as figure 10.

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$.

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). 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^{\lambda },\ldots,u_8^{\lambda }$, at $t=280$ they are not very pronounced yet. Only $u_1^{\lambda }$ stands out through a less pronounced 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^{\lambda }_1$ 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^{\lambda }_1$. 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^{\lambda }_1$ is transverse to the edge manifold at the considered time. The structures highlighted in $u^{\lambda }_8$ 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 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.

Figure 15. Streamwise velocity profiles of $u_1^{\lambda }$ 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$.

Figure 16. Snapshots at $t=720$, $u^{\lambda }_1$. 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_i^{\lambda }\}$, 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$.

Figure 17. Snapshots at $t=720$. Mode $u^{\lambda }_3$. Same as figure 16.

Figure 18. Snapshots at $t=720$. Mode $u^{\lambda }_5$. Same as figure 16.

Figure 19. Snapshots at $t=720$. Mode $u^{\lambda }_8$. 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.

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 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 $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.

Figure 20. Three-dimensional top view of the leading instantaneous eigenmode ${\boldsymbol {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.

Figure 21. Three-dimensional top view of the eighth instantaneous eigenmode ${\boldsymbol {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.

One clear feature from physical space visualisations, regardless of the quantity plotted, is how the localised support of all OTD modes, except here for $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^{\lambda }$ 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, \ldots$ 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:

1. (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);

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

3. (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:

1. (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;

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

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

4. (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;

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

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

7. (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;

8. (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.