A Ginzburg–Landau model for linear global modes in open shear flows

Abstract The Ginzburg–Landau equation (GLE) can phenomenologically model several key features of non-equilibrium systems including those in fluid mechanics. Its validity in real flows, however, remains questionable. Here, we show that the linear GLE can be formulated such that it has the same Wentzel–Kramers–Brillouin (WKB) approximation as for the linear global stability problem in open shear flows. We use the GLE to model the linear global modes of three different wakes and find that it can accurately capture the linear growth rate and frequency to first order in the WKB approximation. Furthermore, we find that it can also provide the shapes of the direct and adjoint eigenvectors and the regions of maximal structural sensitivity. The proposed model requires only the basic flow as input, but gives robust predictions and is computationally inexpensive. As well as opening up new possibilities for GLE-based control strategies, the proposed model makes accurate stability calculations possible, even for some computationally intractable open shear flows.


Introduction
Transition to turbulence is a classic problem in fluid mechanics (Reynolds 1883) and its control is important in many practical applications. For example, one may seek to delay it for drag reduction on an aircraft wing or to promote it for enhanced mixing in a combustor. Although a multi-step process, transition usually begins with the development of hydrodynamic instability on the steady-state basic flow, as shown by the theoretical work of Ruelle & Takens (1971) and by the experiments of Gollub & Swinney (1975). In the present paper, we study one such instability, namely linear global instability, which occurs in a variety of open shear flows (e.g. wakes, jets and mixing layers) and can eventually lead to turbulence via a short sequence of further instabilities.
We study this problem for two reasons. First, linear global instability is known to cause the entire flow to act as a hydrodynamic oscillator with an underlying spatio-temporal structure called a global mode; the Bénard-von Kármán vortex street in a cylinder wake is a classic example (Provansal, Mathis & Boyer 1987). Such global modes can be modelled as pattern forming fronts via the celebrated Ginzburg-Landau equation (GLE) (Newell & Whitehead 1969;Stewartson & Stuart 1971;Moon, Huerre & Redekopp 1983;Cross & Hohenberg 1993;Bohr et al. 1998;Pier & Huerre 2001;Aranson & Kramer 2002;van Saarloos 2003). Moreover, Wentzel-Kramers-Brillouin (WKB) theory shows how the GLE coefficients can be related to local stability analysis (Huerre & Monkewitz 1990;Chomaz, Huerre & Redekopp 1991;Le Dizes et al. 1996), resulting in various GLE-based control strategies (Park, Ladd & Hendricks 1993;Roussopoulos & Monkewitz 1996;Lauga & Bewley 2004;Bagheri et al. 2009;Chen & Rowley 2011;Oehler & Illingworth 2018). However, while qualitatively (and to an extent also quantitatively) successful, the use of the GLE to model global modes remains elusive, even for canonical flows (Pier et al. 1998;Couairon & Chomaz 1999;Juniper & Pier 2015). This limits the direct application of the GLE to real flows (van Saarloos & Hohenberg 1990). The second reason for our studying this problem is that calculations of the WKB approximation, which enable the prediction of linear global instability from the stationary phase argument (Cross & Hohenberg 1993), although mathematically possible (Chomaz et al. 1991;Monkewitz, Huerre & Chomaz 1993), have yet to be achieved owing to their sensitivity to noise in the input data (Le Dizes et al. 1996;Siconolfi et al. 2017). Here, we show that subtle adjustments to the GLE can enable it to accurately model the linear global modes of real flows, providing an innovative way of calculating the WKB approximation for instability predictions.

Global stability problem and its WKB approximation
The linear global instability of a flow is determined by the evolution of infinitesimal perturbations on the steady-state basic flow velocity U. In an incompressible flow, it is governed by the linearized Navier-Stokes equations, where (u, p) (x, y, z, t) are the velocity and pressure perturbations, respectively, is the Laplace operator and Re is the Reynolds number. The global stability problem can be formulated as a generalized eigenvalue problem by decomposing (u, p) = (û,p)(x, y, z) exp(−iωt), ωBq = Lq, (2.2) where the eigenvalue ω gives the temporal evolution, the eigenvectorq = (û,p) gives the spatial structure of the respective velocity and pressure perturbations and, for a given basic flow, the operators B and L are constants (A 1). A flow is then linearly globally unstable if the fastest growing eigenvalue (ω G ) has a positive imaginary part. The validity of the GLE in modelling a global mode (q exp(−iω G t)) is based on the application of WKB theory, which relies on the local quasi-parallel flow assumption (Huerre & Monkewitz 1990;Monkewitz et al. 1993). In parallel flows, where one direction, say x, is homogeneous, the eigenvectorq(x, y, z) becomesq( y, z) exp(ikx), where k = k r + ik i is the streamwise wavenumber, with wave propagation occurring in the x-direction at the group velocity ∂ k ω. Equation (2.2) then becomes ωBq = (L 0 + kL 1 + k 2 L 2 )q, where the operators L 0 , L 1 and L 2 are constants for a given basic flow (A 2). The eigenvectorq can be uniquely determined by the solution pair (ω, k) and is expressed aŝ q (ω, k, y, z). Consequently, the solution around a chosen pair (ω 1 , k 1 ) can be written as a dispersion relation between ω and k, where ω k1 = ∂ k ω(ω 1 , k 1 ), ω kk1 = ∂ 2 k ω(ω 1 , k 1 ) and so on. This is a Taylor series expansion around (ω 1 , k 1 ), which is valid in a small but finite region near (ω 1 , k 1 ) (Schmid & Henningson 2001). A parallel flow is then globally unstable if the disturbances grow in time without advecting away at the group velocity. This is the condition for local absolute instability where the mode with zero group velocity, represented by the solution pair (ω 0 , k 0 ), grows in time (ω 0i > 0) (Huerre & Monkewitz 1985).
WKB theory provides a mathematical framework to extend this simplification to non-parallel flows, where a localized but finite region of local absolute instability can give rise to global instability (Chomaz, Huerre & Redekopp 1988). The global mode can then be formulated as the WKB approximation (Chomaz et al. 1991;Monkewitz et al. 1993), in terms of a series of simple exponentials (Bender & Orszag 1978) where ε is a small parameter representing the slow scale, X = εx, at which the global mode varies in the x-direction. Under the WKB approximation, the wavenumber and the mode shape vary in X, givingq 0 exp iε −1 X k dX − iω g 0 t as the leading-order approximation (zeroth-order solution). The amplitude A and the frequency correction ω ε comprise the second exponential term (first-order solution). On inserting this approximation into the global stability problem (2.2), we get at zeroth order (ε 0 ) where B, L 1 and L 2 are the same as in a parallel flow (A 2), while L n 0 replaces L 0 (A 3). Equation (2.5) is then solved separately at each streamwise station to obtain k at the selected frequency ω g 0 . This selection of the frequency is not trivial and its first mathematically consistent formalism was given by Chomaz et al. (1991). Their criterion gives the leading-order solution ω g 0 as corresponding to the saddle point of ω 0 (X), i.e. ω g 0 = ω 0 (X t ) such that ∂ X ω 0 (X t ) = 0, where, as for a parallel flow, ω 0 (X) corresponds to the local zero-group-velocity mode. The location X t is called the turning point (this term is borrowed from quantum physics) and, in general, it lies in the complex X-plane.
The equation at first order (ε 1 ) is then given as where all the coefficients are functions of X calculated using the equation at zeroth order (2.5), d kk = −∂ 2 k D/∂ ω D and d kω = −∂ 2 kω D/∂ ω D, where D is the local dispersion relation based on (2.5). For the exact derivation of d kk and d kω , please see Monkewitz et al. (1993). We will use the Taylor series expansion of D for the GLE in § 3 (below (3.5)) that will clarify the terms d kk and d kω . Lastly, the term δω represents the contribution from the spatial variations in the local eigenvectorq 0 (ω, k, y, z; X) that do not originate from variations in k. It may also represent the contribution of the terms in the operator L n 0 that are considered too small to be included in the zeroth-order equation (see appendix A).
For the WKB approximation to hold, variations in the amplitude should be slow as compared to its wavenumber. This leads to the condition This condition cannot be satisfied at X t where ω k = 0 and (2.6) is singular. This singularity is of second order (because ∂ X ω 0 = 0) and thus the WKB approximation is invalid in the region ofX = ε −0.5 (X − X t ). The global stability problem must therefore be solved separately in the neighbourhood of X t for the scaled amplitudeĀ and the frequency correctionω ε . The equation in the neighbourhood of X t is then given as This equation is identical to (45) in Huerre & Monkewitz (1990) and (4.13) in Monkewitz et al. (1993). Here k 0X = ∂ X k 0 and ω 0XX = ∂ 2 X ω 0 are calculated at X t . The WKB approximation is finally constructed by solving (2.8) exactly (in terms of Hermite polynomials) and then asymptotically matching the solution in the neighbourhood of X t with the WKB approximations obtained on both sides of, but far away from, X t . Correct branches of the local solutions on either side of X t must be selected such that the boundary conditions ((û,p)(x, y, z) → 0 as |x| → ∞) are satisfied (Monkewitz et al. 1993). The problem, however, is in finding X t (which lies in the complex plane) because it is very sensitive to data on the real axis (Le Dizes et al. 1996).

Ginzburg-Landau equation and its WKB approximation
The GLE has been in use for several decades to model waves in spatially developing flows (Newell & Whitehead 1969;Stewartson & Stuart 1971). Huerre & Monkewitz (1990) used it to illustrate the nature of local convective and absolute instabilities by determining the coefficients of the GLE from local stability analysis, where Ψ represents a global mode with amplitude variation in the x-direction and satisfying the boundary conditions Ψ (x) → 0 as |x| → ∞. In this formulation, the local shape in ( y, z) is determined by the local saddle points (ω 0 , k 0 ) and is obtained as q 0 via (2.5). The same GLE was used by Chomaz et al. (1991) to obtain the frequency selection criterion that is also valid for the Navier-Stokes equations (Monkewitz et al. 1993). Here, we follow Gupta & Wan (2019) by starting with a more general formulation of the GLE where the coefficients are not fixed based on the local saddle points but are left to be determined later As for the global stability problem above, the WKB approximation of Ψ can be constructed as On inserting this approximation into the GLE, we get at leading order (O(ε 0 )) Similar to (2.3), this is a Taylor series expansion of the local dispersion relation (based on (2.5)) around the solution pair (ω 1 , k 1 ). The leading-order solution is governed by the frequency selection criterion of Chomaz et al. (1991), i.e. ω g 0 is determined as the saddle point of ω 0 (X).
At first order (O(ε 1 )), the solution for A and ω ε is determined as where ω k and ω kk are calculated at (ω g 0 , k) from (3.4). Equation (3.5) is equivalent to the corresponding equation for the global stability problem (2.6) such that the terms d kk and d kω are calculated based on a Taylor series expansion of the local dispersion relation (3.4). This gives d kk = ω kk and d kω = 0.
Similar to that for the global stability problem above, the WKB approximation is invalid in the neighbourhood ofX = ε −0.5 (X − X t ). The GLE must therefore be solved separately in this region for the scaled slow amplitudeĀ and the frequency correctionω ε . To this end, the dispersion relation in the neighbourhood of X t is first expanded as where k 0X = ∂ X k 0 , ω kX = ∂ X ω k and ω 0XX = ∂ 2 X ω 0 . The equation in the neighbourhood of X t is then given as This equation is identical to that for the global stability problem (2.8). The boundary conditions (Ψ (x) → 0 as |x| → ∞) are ensured by selecting the correct branches of k 1 on either side of X c (Juniper & Pier 2015). Here X c is the point on the real axis where the branch cut (emanating from X t ) crosses the real X-axis (please see Le Dizes et al. (1996) and Juniper & Pier (2015) for details). The WKB approximation of the global stability problem is based on the exact local stability results (see § 2), whereas the WKB approximation of the GLE is based on a Taylor series expansion (3.4) truncated to the second term. Because the Taylor series expansion is valid only in a small region near the chosen (ω, k), the solution of the GLE (Ψ ) can correctly represent the global mode only if ω 1 is taken to be close to ω g 0 . This explains why the conventional choice of ω 1 = ω 0 (X) (Chomaz et al. 1991;Le Dizes et al. 1996) may not be optimal, because ω 0 can be very different from ω g 0 in most of the flow domain (see second row of figure 1).
The problem with setting ω 1 = ω g 0 , however, is that ω g 0 is not known a priori. In principle, it can be calculated using analytic continuation of ω 0 (X) to the complex X plane (Chomaz et al. 1991). This process, however, is highly sensitive to the data on the real axis (i.e. ill conditioned) (Le Dizes et al. 1996), leading to errors higher than ω ε itself (Siconolfi et al. 2017). In fact, the only known convergent calculations of the WKB approximation require ad hoc smoothing before analytic continuation can be performed (Siconolfi et al. 2017). We propose to sidestep analytic continuation by solving an additional equation that can be considered a third-order GLE, with ξ → 0 as |x| → ∞ and ξ x → 0 as x → ∞. At first order, the solution of the WKB approximation of ξ is identical to that of Ψ but with the underlying Taylor series expansion containing up to the third term We solve (3.2) and (3.8) simultaneously while iterating ω 1 until both equations give identical eigenvalues. Because ω kkk1 is not small, such values of ω 1 must be close to ω g 0 , apart from differences of the order of ∂ 3 X A and ∂ 2 X k, which are expected to be small for slowly varying global modes.

Two-dimensional flows
We first consider two canonical open shear flows: (figure 1a,c,e,g,i) a confined wake and (figure 1b,d, f,h,j) a cylinder wake. Both are two-dimensional with two-dimensional instabilities and have been previously studied by Juniper & Pier (2015) via local stability analysis. The first flow is weakly non-parallel (due to free-slip confinement walls) and consists of three streams injected from the left boundary: a slow inner stream at velocity U 1 (width 2h 1 ) sandwiched between two fast outer streams at velocity U 2 (each of width h 2 ). This flow is defined by the confinement ratio (h ≡ h 2 /h 1 = 1), the inverse shear ratio (Λ −1 ≡ (U 1 − U 2 )/(U 1 + U 2 ) = −1.2), and the Reynolds number (Re ≡ U 2 h 1 /ν, with ν as the kinematic viscosity). The second flow is strongly non-parallel with Re ≡ U ∞ D/ν as the only parameter, where U ∞ is the free-stream velocity and D is the cylinder diameter. The basic flows (top row of 1) are obtained using a spectral element code (Nektar++ Cantwell et al. 2015) with 1600 and 1868 macro elements of polynomial order 7 for the two cases, respectively. For the unstable cases, i.e. above the critical Reynolds number, selective frequency damping is used to obtain the solutions of the steady-state basic flow (Åkervik et al. 2006). We find that the critical Reynolds number for (i) the confined wake is ≈80, which matches well with the global stability results of Juniper, , and for (ii) the cylinder wake is ≈47, which matches well with the experimental results of Provansal et al. (1987).
The local stability code used to obtain the dispersion relations (2.3) andq 0 relies on Chebyshev spectral collocation, in which y is discretized on a Gauss-Lobatto-Chebyshev grid mapped as y = r/(1 − r 2 + (1/y max )), where r is the Chebyshev grid from 0 to 1 (64 points) and y max is 2 for the confined wake but is 20 (an arbitrarily large number) for the cylinder wake. The boundary conditions (BCs) at y max are set as no slip for the confined wake and free slip for the cylinder wake. For the confined wake, no-slip BCs are used for the local stability calculations even though free-slip BCs are used for the base-flow calculations; this is to ensure that the no-penetration condition is satisfied at the confinement wall, consistent with the analysis by Juniper et al. (2011). The second row in figure 1 shows the absolute growth rate (ω 0i ). Both flows contain regions of local absolute instability that grow in length with Re, leading to global instability at a critical Re.
The GLE (3.2) assumes that the domain is infinite (from −∞ to ∞), but neither of the two cases here is truly infinite. This is acceptable so long as the boundaries do not significantly affect the flow instability (Monkewitz et al. 1993). To obtain the GLE, we replace the upstream regions containing abrupt changes (x < 0.1 for the confined wake and x < 0.5 for the cylinder wake) with constant inflow profiles, as per Triantafyllou & Karniadakis (1990). These profiles are extrapolated to an upstream location (x ≈ −3) where the Dirichlet BC (A = 0) is imposed. At a sufficiently far downstream location (x ≥ 8 and 6 for the confined and cylinder wakes, respectively), Neumann BCs (∂ x A = 0, ∂ 2 x A = 0 for (3.8)) are imposed. Lastly, x in (3.2) and (3.8) is discretized with a fourth-order accurate central differencing scheme (dx = 0.1). The resultant problem is solved using the 'eigs' command in MATLAB.
The third and fourth rows in figure 1 show the global stability results obtained from the GLE in terms of k r and k i , respectively. Large negative values of k i near the inlet show that perturbations decay rapidly in the upstream direction, justifying the Dirichlet BC.
Nearly constant values of k r and k i in the downstream region justify the Neumann BCs. The square markers denote the point (X c ) where the branch cut, emanating from X t , crosses the real X-axis. Suitable branches of the local solutions are selected on either side of X c for determining the GLE coefficients (Juniper & Pier 2015). The last row in figure 1 shows the estimated ratio of the second to first exponents of the WKB approximation. This ratio has large values near X c because this is close to X t . For the cylinder wake, this ratio is also large near x = 0.5, indicating that the flow becomes more non-parallel as Re increases. Next, we compare the global modes obtained from the GLE with the exact results obtained by directly solving (2.2) via the Arnoldi method, with a Krylov subspace of size 32, implemented in Nektar++ (He et al. 2017). Figure 2 shows the linear (a) growth rate and (b) frequency. In both wakes, the model results match the exact results well, particularly near the critical Re. The errors, calculated as the difference between the model and the exact ω G (normalized by the exact ω G ), are small near the critical Re, at less than 2.0 % and 2.5 % in the confined and cylinder wakes, respectively. These errors are lower than the 3.1 % error reported by Siconolfi et al. (2017) for a cylinder wake at Re = 50. Crucially, our method does not require any ad hoc smoothing. The errors increase with Re: at approximately twice the critical Re, they reach 2.1 % and 5.5 % in the confined and cylinder wakes, respectively. In the confined wake, the small increase in error with increasing Re is solely due to the difference between the GLE (3.2) and (3.8), which increases as the global mode becomes more unstable. In the cylinder wake, the error increases mainly because the flow becomes more non-parallel at higher Re (see the last row of figure 1).

Direct and adjoint eigenvectors
The top row of figure 3 shows the spatial structure of the global modes in terms of the y-velocity perturbations (the direct eigenvectors), representing the dominant linear vortex shedding modes in the cylinder wake. The upstream boundary in these figures is fixed at x = 0.5, exactly at the downstream end of the cylinder. The eigenvectors from the model match the exact results well, albeit the shapes from the model are marginally wider in the y-direction and their amplitudes are slightly shifted in the downstream direction. This downstream shift of the global modes from the model is consistent with the under-prediction of frequency ω by the model. The phase speed (ω/k) of an instability wave is determined by the basic flow velocity, which is fixed. Under-prediction of ω is therefore expected to lead to an over-prediction of the wavelength (inverse of k). We also note here that, although the linear global-mode frequency, growth rate and amplitude variation are calculated to first order in the WKB approximation, the mode shape is calculated to only zeroth order. A further correction to the mode shape, equivalent to first order in the WKB approximation, can be calculated by solving an extended eigenvalue problem (termed Model + EEV, shown in the second row of figure 3) as in Huang & Wu (2015). This calculation procedure, however, is not always robust and is only applicable in the downstream region where ∂ x k, and hence the correction, is small.
The direct eigenvectors indicate where the instability is most observable and hence where sensors could be placed. For controller design, we also need the adjoint  eigenvectors, which indicate where actuators could be placed. The third row of figure 3 shows the adjoint eigenvectors (also in terms of the y-component of velocity) from the model, calculated to leading order as in Juniper & Pier (2015), and its comparison with the exact adjoint eigenvector, calculated by solving the adjoint of (2.2) (Giannetti & Luchini 2007). Finally, the last row of figure 3 shows the overlap between the direct and adjoint eigenvectors. This overlap, called the structural sensitivity region, indicates where the instability is generated, i.e. the wavemaker region (Giannetti & Luchini 2007). The square markers in the figure denote regions of maximal structural sensitivity, which is where a passive control device, such as a small control cylinder (Strykowski & Sreenivasan 1990), would have the largest effect on the global instability (Giannetti & Luchini 2007); these results could be further improved by including the sensitivity to base-flow modifications, as shown by Marquet, Sipp & Jacquin (2008) and Pralits, Brandt & Giannetti (2010). The location of maximal structural sensitivity is predicted very accurately at Re = 40 and reasonably well at Re = 80. As is the case with the eigenvalue results, the errors increase as the flow becomes more non-parallel at higher Re.

Application to mean flow stability analysis
Linear global modes obtained from basic flow stability analysis are useful because they describe the evolution of infinitesimal perturbations in a flow . Such modes, however, often fail to accurately characterize the nonlinearly saturated finite-amplitude perturbations. This has led to interest in the application of stability analysis on the time-averaged (mean) flow; please see the early work of Malkus (1956) and Stuart (1958) as well as many subsequent studies focusing on cylinder wakes (Zielinska et al. 1997;Pier 2002;Noack et al. 2003;Barkley 2006). In seminal work, Stuart (1958) observed that the growth of perturbations is initially governed by linear theory, but once the perturbations become sufficiently large, they distort the mean flow such that an equilibrium state is reached where the energy transfer from the mean flow to the perturbations is balanced precisely by the energy dissipated via viscous effects. The idea, therefore, is that mean flow stability analysis can account for this energy transfer and can thus characterize the nonlinearly saturated finite-amplitude perturbations.
In figure 4(a), we present the basic flow (top half) and the mean flow (bottom half) for the cylinder wake at Re = 80. Vortex shedding in the cylinder wake distorts the basic flow such that the recirculation bubble in the mean flow becomes shorter. Increasing Re lengthens the recirculation bubble in the basic flow, but shortens it in the mean flow. This opposing trend was attributed by Zielinska et al. (1997) to a mean flow correction caused by the fundamental unstable mode. Consistent with this observation, figure 4(b) shows that while the region of local absolute instability grows in length with Re for the basic flow, it shrinks in length with Re for the mean flow. This highly nonlinear effect can also be seen in figure 4(e), where the deviation of the nonlinear vortex-shedding frequency from the linear global-mode frequency grows rapidly with Re. Again, the local stability results shown in figure 4(c) capture this frequency shift, which was also reported by Pier (2002). Finally, Barkley (2006) performed a global stability analysis on the mean flow of a cylinder wake and found that (i) the obtained eigenfrequency exactly tracks the Strouhal number of nonlinear vortex shedding and (ii) the mean flow is marginally stable (ω Gi ≈ 0).
An open question is whether the present framework can be extended to obtain mean flow global stability results similar to those reported by Barkley (2006). To answer this, we apply the same procedure as explained above but with the basic flow replaced by the mean flow to obtain ω G . The results are presented as symbols in figures 4(d) and 4(e). We find that the present framework slightly under-predicts the nonlinear frequency (figure 4e) but over-predicts the stabilizing effect of the mean flow distortion (figure 4d). Similar to the basic flow results, the mean flow results are better near the critical Re but the errors increase as the flow becomes more non-parallel at higher Re.
From these results, we can conclude that the present framework is also applicable to mean flow stability analysis provided that the constraints for slow evolution in the streamwise direction are satisfied, much as they were for basic flow stability analysis. We must, however, caution that the mean flow is not generally a solution of the Navier-Stokes equations. The results from mean flow stability analysis, therefore, should be carefully interpreted. Sipp & Lebedev (2007) showed that the mean flow is only marginally stable when the higher harmonics are not present; otherwise the energy transfer between different harmonics can affect the evolution of nonlinearly saturated perturbations, an effect that cannot be accounted for by mean flow stability analysis. We also note that another field of application of mean flow stability analysis is turbulent flows. The validity conditions for turbulent mean flow stability analysis have been provided by Beneddine et al. (2016) while the use of the present framework to model convective instabilities in a turbulent mean flow has been demonstrated by Gupta & Wan (2019).

A three-dimensional flow: the flow past a sphere
Lastly, we demonstrate the applicability of the GLE for modelling the leading linear global mode and calculating the linear global stability results in the wake of a sphere at Re ≡ U ∞ D/ν = 300, where D is the sphere diameter and U ∞ is the free-stream velocity. This three-dimensional (3-D) flow undergoes a supercritical pitchfork bifurcation  at Re ≈ 212, after which it loses axisymmetry but gains planar symmetry (Fabre, Auguste & Magnaudet 2008). It then undergoes a supercritical Hopf bifurcation at Re ≈ 275, transitioning to a limit cycle characterized by periodic vortex shedding (Tomboulides & Orszag 2000). Direct global stability analysis of three-dimensional flows is still considered computationally expensive (Theofilis 2011). Only recently have researchers reported direct global stability results for non-axisymmetric sphere wakes (Citro et al. 2016) and the corresponding WKB approximation (Siconolfi et al. 2017). Sphere wakes, therefore, represent a challenging yet verifiable flow on which the present GLE-based model can demonstrate its computational efficiency over the direct method as well as its advantage over the conventional WKB approximation involving analytic continuation of local stability results.
The basic flow for the sphere wake (figure 5a) is obtained using Nektar++ with 4700 hexahedral elements of polynomial order 5. The corresponding direct stability results (figure 5b) are obtained by discretizing the spatial operator using 587 500 points (see He et al. (2019) for details). The present method, by contrast, involves performing biglobal stability analysis at each x-location, requiring spatial discretization in only the y-and z-directions. Here we use 36 points in y from 0 to 3 and 72 points in z from −3 to 3, all discretized on a Gauss-Lobatto-Chebyshev grid, as in § 4.1. Global stability results are then obtained by solving the GLE, which requires discretization in only the x-direction. The present method, therefore, effectively reduces a problem of size Linear global mode ( f ) k i and (g) k r obtained from the model (square markers denote X c ).
GLE remains the same as for the two-dimensional flows of § 4.1. In the implementation, however, we do not directly solve the third-order GLE (3.8), which requires higher spatial resolution for the biglobal stability calculations; instead we insert the solution of the GLE into (3.8) and minimize the residual in the neighbourhood of X c . The global mode calculated from the GLE (bottom half of figure 5(b) and its 3-D view in figure 5h) is in good agreement with the directly calculated global mode (top half of figure 5(b) and its 3-D view in figure 5c), particularly in the region of local absolute instability (i.e. upstream of x ≈ 2). It is important to note that the linear global mode here varies more rapidly than the linear global modes in the two-dimensional wakes of § 4.1, particularly downstream of x ≈ 2 (compare k r in figure 5(g) with those in the third row of figure 1). The first-order corrections to the mode shape (the second row of figure 3) are therefore expected to be large, explaining the mismatch between the model and exact linear global mode shapes in the downstream region. The eigenvalue calculated from the GLE (ω G = 0.79 + 0.08i) is in good agreement with the direct results (ω G = 0.81 + 0.10i) with less than 4 % error. Such agreement is reassuring considering that this sphere wake is much more non-parallel than the two-dimensional wakes considered in § 4.1. However, the global mode shapes calculated from analytic continuation of the local stability results of Siconolfi et al. (2017) (figures 9(a) and 9(b) of their paper) match less well with the exact solutions, even in the region of local absolute instability (figures 9(e) and 9( f ) of their paper).
Analytic continuation of local stability results to the complex plane is an ill-posed and ill-conditioned problem (Le Dizes et al. 1996). The WKB approximation of the stability results calculated by Siconolfi et al. (2017) via analytic continuation is therefore impressive, particularly for the cylinder wake. For the sphere wake, however, analytic continuation becomes even more ill conditioned because this flow has a much smaller pocket of local absolute instability and is very non-parallel. Figure 5(d) shows that the pocket of local absolute instability in this flow extends to only x ≈ 2, downstream of which ω 0i quickly decreases to negative values. This implies that the local absolute instability results downstream of x = 2 quickly become irrelevant to the global mode. This can also be seen from figure 5( f ) where k i is already very positive at x = 2 and consequently the global-mode amplitude quickly decreases downstream of x ≈ 2. The functions ω 0 , k 0 , ω kk and their derivatives, which are required for analytic continuation, are therefore only available on a thin strip from x = 0.5 to x ≈ 2.0-3.0. When analytic continuation is performed away from the strip in the complex plane, the errors are known to increase linearly with distance (normalized by the strip length) from the strip (Trefethen 2020), thus limiting the use of analytic continuation for calculating the WKB approximation in such highly non-parallel flows.
We would like to take this opportunity to clarify the suitability of the expression used by Siconolfi et al. (2017) to calculate the corrected global-mode shapes (figures 9(c) and 9(d) of their paper). In that expression, the spatial mode shape et al. 2017). This expression, however, is inconsistent with the WKB approximation in which the first-order corrections come from the variation in local stability results, such as the contribution of ∂k(X, ω g 0 )/∂ x. The first-order correction term A(X) is given by Monkewitz et al. (1993), which explicitly includes ∂k(X, ω g 0 )/∂ x among other terms (see (2.6)). Giannetti & Luchini (2006) and Viola, Arratia & Gallaire (2016) have presented solvability conditions for calculating such first-order correction terms for the forced response of locally convectively unstable flows. In the present framework, the correction term (2.6) is simplified as (3.5) and is included in the global mode obtained by solving the GLE. The caveat, however, is that there is no unique choice as to which physical variable A(X) relates to. For example, in the sphere wake, it is related to the norm of the z-velocity perturbations. Because all the velocity perturbations follow nearly the same trend (i.e. reaching a maximum at around x ≈ 2), the global-mode shape is not significantly affected by this choice.

Discussion
The present framework predicts the WKB approximation to global stability results for open shear flows. The implementation is straightforward as it involves determining the coefficients of the GLE from local stability analysis and avoids analytic continuation, which is an ill-posed and ill-conditioned problem. However, care must be taken to ensure that a smooth transition of the local wavenumber k occurs across X c (where the branch cut crosses the real X-axis), the boundaries do not significantly affect the results, and the imaginary part of ω kk1 remains mostly negative, satisfying causality (Chomaz et al. 1991). We find in § 4 that the error in the global stability predictions remains within approximately 5 % even for highly non-parallel flows at Reynolds numbers twice the critical values. Such levels of error in global stability calculations are acceptable because they can also arise from subtle differences in the BCs and fluid properties between experiments and simulations. The application of the present framework, therefore, also depends on its ability to handle such uncertainties. It can easily be demonstrated that the GLE-based model provides robust predictions so long as the large-scale features of the obtained coefficients remain unchanged. We have performed a quick test of this for the two-dimensional basic flow cases by adding to the coefficients k 1 , ω k1 and ω kk1 random noise with a normal distribution and a standard deviation of 5.0 % (normalized by the spatial average of the coefficient). The resultant change in the ω G predictions has a normalized standard deviation of under 0.5 % in all the cases studied.
The present framework also models the linear global modes in terms of the GLE, which have been the subject of many studies on flow control, e.g. Park et al. (1993), Roussopoulos & Monkewitz (1996), Lauga & Bewley (2004), Bagheri et al. (2009), Chen & Rowley (2011 and Oehler & Illingworth (2018). In particular, Lauga & Bewley (2004) and Bagheri et al. (2009) studied the control of Ginzburg-Landau systems within the framework of linear time-invariant control theory. Later, Chen & Rowley (2011) and Oehler & Illingworth (2018) used a similar framework to explore optimal sensor and actuator placement in spatially developing flows. Such interest in the GLE arises from the fact that it can model the instabilities in spatially developing flows (Chomaz 2005;Bagheri et al. 2009) at a cost low enough to make feedback control viable (Oehler & Illingworth 2018). Its applications, however, have until now remained only qualitative.
Here, we show that the GLE can quantitatively model flow instabilities in spatially developing flows, thus paving the way for the application of GLE-based control strategies to increasingly realistic flows. To realize this, however, one must first answer three questions: (i) How can the actuator inputs or the sensor outputs be transformed into the GLE variable (i.e. Ψ in (3.2)) (Bagheri et al. 2009)? (ii) Can the GLE work on non-modal (Chen & Rowley 2011), quasi-periodic (Leclercq et al. 2019) or multi-mode systems (Carini, Auteri & Giannetti 2015)? (iii) Can nonlinear effects be ignored (Leclercq et al. 2019)? These questions are generally problem specific but have been addressed in great detail by Bagheri et al. (2009). Furthermore, the issue of projecting the full system onto a reduced-order system has been addressed by Sipp & Schmid (2016).

Conclusions
Figures 2 and 3 capture the main achievements of this paper with regard to the prediction of global instability and the modelling of linear global modes in two-dimensional flows, while figures 4 and 5 extend the applicability to time-averaged and highly non-parallel three-dimensional flows, respectively. The two main strengths of the GLE-based model are that (i) it requires as input the basic flow in only a limited part of the domain in order to be able to robustly predict the linear global mode, and (ii) it is computationally efficient, reducing a triglobal stability problem to a series of simpler biglobal stability problems. The proposed model therefore opens up new possibilities for GLE-based linear control strategies to be directly applied to realistic flows. Crucially, the model makes accurate instability predictions possible even for flows that would otherwise be computationally intractable.
It is worth noting that although the present paper focuses on linear unforced global modes, the GLE has seen a wide range of applications in flow instability -and pattern formation in general. The GLE has been shown to be able to (i) qualitatively model fully nonlinear global modes (Pier & Huerre 1996;Pier et al. 1998;Couairon & Chomaz 1999), (ii) capture secondary instabilities in wakes (Eckhaus and Benjamin-Feir instabilities) (Leweke, Provansal & Boyer 1993;Leweke & Provansal 1994), (iii) predict the nature of bifurcations (Chomaz 1992), (iv) predict the selection of solutions (van Saarloos & Hohenberg 1990;Chomaz & Couairon 2000), (v) capture the transient growth and forced response of convectively unstable flows (Deissler 1987a;Cossu & Chomaz 1997;Gupta & Wan 2019) and (vi) capture turbulent bursts in boundary layers (Deissler 1987b) and transition to turbulence in open shear flows (Moon et al. 1983). Extending the model such that it can reproduce these features in real flows remains an important open challenge. where L = −iU · ∇ + (i/Re) .
where C = −i(V∂ y + W∂ z ). It should be noted that the terms involving V, W, the derivative of the basic flow in x (∂ x U, ∂ x V and ∂ x W), and the 1/Re terms (D and the operator L 2 ) are considered to be of order ε or ε 2 . Consequently, these terms are not included in the zeroth-order solution of Monkewitz et al. (1993). Most other studies include the 1/Re terms in the zeroth-order solution but leave the other terms out. This, however, is a matter of choice. If these terms are accurately known, they can be included in the zeroth-order solution at no extra cost, as shown by Siconolfi et al. (2017).