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 Reference Drazin and Reid2004). 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 Reference Morkovin1969; Henningson, Lundbladh & Johansson Reference Henningson, Lundbladh and Johansson1993; Jacobs & Durbin Reference Jacobs and Durbin2001; Matsubara & Alfredsson Reference Matsubara and Alfredsson2001; Brandt & Henningson Reference Brandt and Henningson2002; Brandt, Schlatter & Henningson Reference Brandt, Schlatter and Henningson2004) 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 Reference Klebanoff, Tidstrom and Sargent1962). They are predominantly streamwise-independent structures supporting three-dimensional wiggles convected at different velocities (Brandt et al. Reference Brandt, Cossu, Chomaz, Huerre and Henningson2003). 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 Reference Schmid and Henningson2001) via a mechanism called lift-up. This mechanism transfers streamwise vorticity upstream into streaks further downstream (Landahl Reference Landahl1980; Brandt Reference Brandt2014). Careful early experiments have suggested that their breakdown follows an instability mechanism (Bakchinov et al. Reference Bakchinov, Grek, Klingmann and Kozlov1995; Alfredsson & Matsubara Reference Alfredsson and Matsubara1996; Matsubara & Alfredsson Reference Matsubara and Alfredsson2001; Asai, Minagawa & Nishioka Reference Asai, Minagawa and Nishioka2002). 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. Reference Andersson, Brandt, Bottaro and Henningson2001; Kawahara et al. Reference Kawahara, Jiménez, Uhlmann and Pinelli2003; Brandt Reference Brandt2007). 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 (Reference Hamilton, Kim and Waleffe1995) 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. Reference Hamilton, Kim and Waleffe1995; Waleffe Reference Waleffe1997). Hœpffner, Brandt & Henningson (Reference Hœpffner, Brandt and Henningson2005), following Schoppa & Hussain (Reference Schoppa and Hussain2002), 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 (Reference Cossu, Chevalier and Henningson2007). Schlatter et al. (Reference Schlatter, Brandt, De Lange and Henningson2008) studied the secondary instability of streaks via nonlinear impulse response. Linear stability features were later extracted directly from numerical data (Vaughan & Zaki Reference Vaughan and Zaki2011; Hack & Zaki Reference Hack and Zaki2014) 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 Reference Alizard2015; Cassinelli, de Giovanetti & Hwang Reference Cassinelli, de Giovanetti and Hwang2017). 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 Reference Viana2014). The generalisation of the eigenvectors is given by the (time-dependent) covariant Lyapunov vectors (CLVs); see Ginelli et al. (Reference Ginelli, Poggi, Turchi, Chaté, Livi and Politi2007), Kuptsov & Parlitz (Reference Kuptsov and Parlitz2012) and Pikovsky & Politi (Reference Pikovsky and Politi2016). 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 Reference Kerswell2018). 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 Reference Kerswell2018). Its computation is based on a nonlinear optimisation method (Pringle & Kerswell Reference Pringle and Kerswell2010; Cherubini et al. Reference Cherubini, De Palma, Robinet and Bottaro2011b; Vavaliaris, Beneitez & Henningson Reference Vavaliaris, Beneitez and Henningson2020) 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 Reference Itano and Toh2001): 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 Reference Biau2012; Khapko et al. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013, Reference Khapko, Duguet, Kreilos, Schlatter, Eckhardt and Henningson2014, Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016) but less straightforward in spatially developing boundary layer flows like the Blasius boundary layer (Cherubini et al. Reference Cherubini, De Palma, Robinet and Bottaro2011b; Duguet et al. Reference Duguet, Schlatter, Henningson and Eckhardt2012). 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. Reference Beneitez, Duguet, Schlatter and Henningson2019; Beneitez, Duguet & Henningson Reference Beneitez, Duguet and Henningson2020a). 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 Reference Haller2015). 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 (Reference Babaee and Sapsis2016). 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 (Reference Blanchard and Sapsis2019) 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. Reference Babaee, Farazmand, Haller and Sapsis2017). 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 (Reference Babaee and Sapsis2016) 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 Reference Farazmand and Sapsis2016), the control of linear instabilities (Blanchard, Mowlavi & Sapsis Reference Blanchard, Mowlavi and Sapsis2019) and the stability of pulsating Poiseuille flow (Kern et al. Reference Kern, Beneitez, Hanifi and Henningson2021) as well as for faster edge tracking in high dimension (Beneitez et al. Reference Beneitez, Duguet, Schlatter and Henningson2020b). 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ć Reference Gibson, Halcrow and 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)$
\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
\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 Reference Pikovsky and Politi2016). 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
\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. Reference Ginelli, Poggi, Turchi, Chaté, Livi and Politi2007; Kuptsov & Parlitz Reference Kuptsov and Parlitz2012); CLVs are thus essentially an inappropriate computational tool for the study of transients.
2.3. Optimally time-dependent modes
Babaee & Sapsis (Reference Babaee and Sapsis2016) 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
\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 Reference Blanchard and Sapsis2019; Blanchard et al. Reference Blanchard, Mowlavi and Sapsis2019).
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 (Reference Babaee and Sapsis2016), the maximisation of 
$\mathcal {J}$ in (2.3) under the orthogonality constraint (2.4) yields a system of coupled nonlinear evolution equations
\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 (Reference Blanchard and Sapsis2019), a particular choice of 
$\boldsymbol {A}$ was made
\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
\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 (Reference Babaee and Sapsis2016) introduced the reduced operator 
${\boldsymbol {L}}_r$ defined by projecting the high-dimensional operator 
${\boldsymbol {L}}$ onto the OTD directions
\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 (Reference Babaee and Sapsis2016), 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
\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 Reference Embree and Trefethen2005). 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
\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
\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 Reference Haller2015). 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 Reference Haller2015), as
\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. (Reference Babaee, Farazmand, Haller and Sapsis2017) 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. Reference Babaee, Farazmand, Haller and Sapsis2017; Sapsis Reference Sapsis2018; Blanchard & Sapsis Reference Blanchard and Sapsis2019), 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 Reference Blanchard and Sapsis2019). They are expressed as
\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 Reference Shimada and Nagashima1979). 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 Reference Blanchard and Sapsis2019). The same modes have also been called sometimes backwards Lyapunov vectors Kuptsov & Parlitz (Reference Kuptsov and Parlitz2012). 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
\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
\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
\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. (Reference Chevalier, Schlatter, Lundbladh and Henningson2007). The streamwise component of 
$\mathcal {U}(x,y,z)$ is defined as
\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. (Reference Chevalier, Schlatter, Lundbladh and Henningson2007), and in several later publications including Duguet et al. (Reference Duguet, Schlatter, Henningson and Eckhardt2012) and Beneitez et al. (Reference Beneitez, Duguet, Schlatter and Henningson2019). 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 (Reference Nordström, Nordin and Henningson1999).
The temporal integration of the incompressible Navier–Stokes equations is performed using the pseudo-spectral solver SIMSON (Chevalier et al. Reference Chevalier, Schlatter, Lundbladh and Henningson2007). 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. (Reference Beneitez, Duguet, Schlatter and Henningson2019), 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. (Reference Vavaliaris, Beneitez and Henningson2020). 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. (Reference Duguet, Schlatter, Henningson and Eckhardt2012) and Beneitez et al. (Reference Beneitez, Duguet, Schlatter and Henningson2019). The computation of the OTD modes starts at initial time 
$t=0$ from the (spatially localised) minimal seed computed in Vavaliaris et al. (Reference Vavaliaris, Beneitez and Henningson2020) 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
\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. Reference Cherubini, De Palma, Robinet and Bottaro2011a; Vavaliaris et al. Reference Vavaliaris, Beneitez and Henningson2020), in particular the original work by Cherubini et al. (Reference Cherubini, De Palma, Robinet and Bottaro2011a) 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. Reference Vavaliaris, Beneitez and Henningson2020; Beneitez et al. Reference Beneitez, Duguet and Henningson2020a). In Vavaliaris et al. (Reference Vavaliaris, Beneitez and Henningson2020) 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 (Reference Rabin, Caulfield and Kerswell2012), Kerswell (Reference Kerswell2018). 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 (Reference Foures, Caulfield and Schmid2013) and is based on the implementation into the open-source solver Nek5000 originally implemented by Rinaldi, Canton & Schlatter (Reference Rinaldi, Canton and Schlatter2019). The optimal state determined for a near-to-threshold initial energy 
$E_0$ was bisected using an edge tracking algorithm (Itano & Toh Reference Itano and Toh2001; Skufca, Yorke & Eckhardt Reference Skufca, Yorke and Eckhardt2006), 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. (Reference Cherubini, De Palma, Robinet and Bottaro2011a). Although the investigation in Beneitez et al. (Reference Beneitez, Duguet, Schlatter and Henningson2019) 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. Reference Duguet, Schlatter, Henningson and Eckhardt2012; Beneitez et al. Reference Beneitez, Duguet, Schlatter and Henningson2019, Reference Beneitez, Duguet and Henningson2020a)
$$\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}$$
$$\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}$$
$$\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. (Reference Duguet, Schlatter, Henningson and Eckhardt2012). 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. (Reference Vavaliaris, Beneitez and Henningson2020). 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 (Reference Jeong and Hussain1995). 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 (Reference Jeong and Hussain1995). 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. (Reference Vavaliaris, Beneitez and Henningson2020). For early times 
$t \leqslant 60$ the dynamics is dominated by a three-dimensional version of the Orr mechanism (Vavaliaris et al. Reference Vavaliaris, Beneitez and Henningson2020), 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 Reference Schmid and Henningson2001). 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. Reference Beneitez, Duguet, Schlatter and Henningson2019) 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. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016). 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 Reference Babaee and Sapsis2016). 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 Reference Duguet, Willis and Kerswell2008). 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. Reference Vavaliaris, Beneitez and Henningson2020). 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. (Reference Babaee, Farazmand, Haller and Sapsis2017). 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. (Reference Vavaliaris, Beneitez and Henningson2020). 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. (Reference Vavaliaris, Beneitez and Henningson2020), 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. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013; Beneitez et al. Reference Beneitez, Duguet, Schlatter and Henningson2019), 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
\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 (Reference Vaughan and Zaki2011), 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 (Reference Vaughan and Zaki2011), 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 (Reference Vaughan and Zaki2011), 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. Reference Babaee, Farazmand, Haller and Sapsis2017). 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 Reference Itano and Toh2001), 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. (Reference Duguet, Schlatter, Henningson and Eckhardt2012).
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
\begin{equation} D_{KY}=j+\frac{S_j}{|\lambda_{j+1}|}, \end{equation}
where 
$S_j$ the cumulative sum
\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 Reference Kuptsov and Kuznetsov2018). 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 Reference Kuptsov and Kuznetsov2018). 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. Reference Vavaliaris, Beneitez and Henningson2020) 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 (Reference Balamurugan and Mandal2017), 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 (Reference Babaee and Sapsis2016). 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. (Reference Andersson, Brandt, Bottaro and Henningson2001) 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. (Reference Andersson, Brandt, Bottaro and Henningson2001), 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. Reference Duguet, Schlatter, Henningson and Eckhardt2012). 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 (Reference Vaughan and Zaki2011). The relevance of this region is furthermore consistent with the interpretation in Hack & Zaki (Reference Hack and Zaki2014), 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. (Reference Vavaliaris, Beneitez and Henningson2020) 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. Reference Vavaliaris, Beneitez and Henningson2020) 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. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2013, Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016; Beneitez et al. Reference Beneitez, Duguet, Schlatter and Henningson2019), ASBL (Khapko et al. Reference Khapko, Kreilos, Schlatter, Duguet, Eckhardt and Henningson2016) as well as channel flow (Toh & Itano Reference Toh and Itano2003). 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. Reference Mellibovsky, Meseguer, Schneider and Eckhardt2009; Duguet, Willis & Kerswell Reference Duguet, Willis and Kerswell2010; Duguet et al. Reference Duguet, Monokrousos, Brandt and Henningson2013).
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 (Reference Vaughan and Zaki2011) and Hack & Zaki (Reference Hack and Zaki2014), even if the corresponding perturbations lack the long wavelength structure characteristic of streak eigenmodes reported so far (Andersson et al. Reference Andersson, Brandt, Bottaro and Henningson2001; Brandt Reference Brandt2014). 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 Reference Schmid and Henningson2001). 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. Reference Duguet, Schlatter, Henningson and Eckhardt2012).
Although the method originally targets a modal description of the relevant finite-time instabilities, it can also capture non-normal amplification mechanisms (Babaee & Sapsis Reference Babaee and Sapsis2016). 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 Reference Schoppa and Hussain2002; Hœpffner et al. Reference Hœpffner, Brandt and Henningson2005).
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. Reference Babaee, Farazmand, Haller and Sapsis2017) 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 (Reference Pikovsky and Politi2016). 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 Reference Babaee and Sapsis2016; Kern et al. Reference Kern, Beneitez, Hanifi and Henningson2021). 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 Reference Cvitanović, Davidchack and Siminos2010) 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. Reference Kern, Beneitez, Hanifi and Henningson2021). However, it has been shown that much lower number of modes
$r\approx 6$ could already bring relevant physical insight (Kern et al. Reference Kern, Beneitez, Hanifi and Henningson2021);(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:
(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 (Reference Schmid2007) 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. (Reference Cossu, Chevalier and Henningson2007) 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.
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. Reference Chevalier, Schlatter, Lundbladh and Henningson2007). Consider the equations for the evolution of the OTD modes about a trajectory evolved with the Navier–Stokes equations
\begin{equation} \boldsymbol{\dot{u}}_i = \boldsymbol{L}_{NS} (\boldsymbol{u}_i)-\langle \boldsymbol{L}_{NS} (\boldsymbol{u}_i), \boldsymbol{u}_i \rangle \boldsymbol{u}_i-\sum_{j=1}^{i-1} [\langle \boldsymbol{L}_{NS} (\boldsymbol{u}_i), \boldsymbol{u}_j\rangle + \langle \boldsymbol{L}_{NS} (\boldsymbol{u}_j), \boldsymbol{u}_i\rangle ]\boldsymbol{u}_j, \end{equation}
where 
$\boldsymbol {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 
$\boldsymbol {u}_i$ reads
\begin{equation} \boldsymbol{L}_{NS} (\boldsymbol{u}_i) =- ({\boldsymbol{U}}_b \boldsymbol{\cdot}\boldsymbol{\nabla}) {\boldsymbol{u}}_i - ({\boldsymbol{u}}_i \boldsymbol{\cdot}\boldsymbol{\nabla}) {\boldsymbol{U}}_b - \boldsymbol{\nabla} p_i + \frac{1}{\mbox{Re}} \nabla^2 {\boldsymbol{u}}_i. \end{equation} 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 
$\boldsymbol {L}_{NS}(\boldsymbol {u}_i)$ needs to contain the boundary conditions corresponding to the linearised operator to provide a correct forcing term and 
${\boldsymbol {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
\begin{equation} {\boldsymbol{u}}_i(y=0)=\frac{\partial {\boldsymbol{u}}_i}{\partial y}(y=L_y) = 0, \end{equation}
when recovering the wall-normal velocity from the fourth-order equation arising from the velocity–vorticity formulation (Chevalier et al. Reference Chevalier, Schlatter, Lundbladh and Henningson2007). 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 
$\boldsymbol {L}_{NS}(\boldsymbol {u}_i)$ as well as the implicit part of the solver.
Figure 22. Snapshots at 
$t=120$. Contours of the wall-normal velocity 
$V$ for 
$u^{\lambda }_1$ 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^{\lambda }\}$ are not normalised.
The OTD modes converge exponentially fast to the most unstable directions of the Cauchy–Green tensor (Babaee et al. Reference Babaee, Farazmand, Haller and Sapsis2017), and after a long time only depend on the point of the trajectory where they are computed (Blanchard & Sapsis Reference Blanchard and Sapsis2019). However, the OTD modes depend on their initialisation (Babaee & Sapsis Reference Babaee and Sapsis2016; Kern et al. Reference Kern, Beneitez, Hanifi and Henningson2021). 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 Reference Babaee and Sapsis2016).
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 Reference Babaee and Sapsis2016).
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 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.
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.
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 Reference Babaee and Sapsis2016; Babaee et al. Reference Babaee, Farazmand, Haller and Sapsis2017; Blanchard & Sapsis Reference Blanchard and Sapsis2019), but the explicit times are strongly case dependent.
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$.
Open access
