Second-order adjoint-based sensitivity for hydrodynamic stability and control

Adjoint-based sensitivity analysis is routinely used today to assess efficiently the effect of open-loop control on the linear stability properties of unstable flows. Sensitivity maps identify regions where small-amplitude control is the most effective, i.e. yields the largest first-order (linear) eigenvalue variation. In this study an adjoint method is proposed for computing a second-order (quadratic) sensitivity operator, and applied to the flow past a circular cylinder, controlled with a steady body force or a passive device model. Maps of second-order eigenvalue variations are obtained, without computing controlled base flows and eigenmodes. For finite control amplitudes, the second-order analysis improves the accuracy of the first-order prediction, and informs about its range of validity, and whether it underestimates or overestimates the actual eigenvalue variation. Regions are identified where control has little or no first-order effect but a second-order effect. In the cylinder wake, the effect of a control cylinder tends to be underestimated by the first-order sensitivity, and including second-order effects yields larger regions of flow restabilisation. Second-order effects can be decomposed into two mechanisms: second-order base flow modification, and interaction between first-order modifications of the base flow and eigenmode. Both contribute equally in general in sensitive regions of the cylinder wake. Exploiting the second-order sensitivity operator, the optimal control maximising the total second-order stabilisation is computed via a quadratic eigenvalue problem. The approach is applicable to other types of control (e.g. wall blowing/suction and shape deformation) and other eigenvalue problems (e.g. amplification of time-harmonic perturbations, or resolvent gain, in stable flows).


Introduction
Over the past decades, adjoint-based sensitivity analysis has become a standard tool for estimating the effect of flow control. The key underlying idea is to compute the gradient of a quantity of interest with respect to control by solving so-called adjoint equations, only once. This approach contrasts with the brute-force method, where the gradient is obtained by solving the direct equations (e.g. Navier-Stokes equations) once for each control degree of freedom. When the control has many degrees of freedom, for instance when it depends on space or time, the adjoint method dramatically reduces the computational cost. This efficient calculation is crucial in iterative gradient-based methods for optimal control, where the gradient is repeatedly evaluated. This is true in general for systems governed by partial differential equations (Lions 1971), and in particular for a wide range of problems in fluid mechanics: shape optimisation for E. Boujo aerodynamics or mixing (Jameson et al. 1998;Mohammadi & Pironneau 2001); optimal wall actuation for turbulent drag reduction (Bewley et al. 2001) or mixing (Foures et al. 2014); optimal kinematics for thin-film coating (Boujo & Sellier 2019); and optimal perturbations (initial perturbations undergoing the largest possible transient growth), especially for time-varying base flows or nonlinear amplification (Schmid 2007), the latter being relevant to transition to turbulence (Pringle & Kerswell 2010;Monokrousos et al. 2011).
Adjoint equations also appear naturally in fluid mechanics when investigating how linear stability properties (growth rate and frequency, characterised by a linearised eigenvalue problem) are affected by flow control (Luchini & Bottaro 2014). Sensitivity maps are obtained that allow one to identify the most sensitive regions at a glance and thus to design effective controls easily. This approach is very efficient: unlike trialand-error techniques, it never actually solves for controlled flows, and only requires one adjoint calculation. The method has been applied extensively in the flow past a circular cylinder, a prototypical globally unstable open flow: the sensitivity of the leading eigenvalue has been computed with respect to passive control (namely, a model of a small secondary cylinder acting on both the base flow and the perturbations) (Hill 1992), to a localised feedback force proportional to the perturbation flow velocity (Giannetti & Luchini 2007), and to flow modification and steady forcing in the bulk (Marquet et al. 2008). To some extent, these studies correctly identified restabilising regions where vortex shedding is suppressed by a small secondary cylinder, first identified by the systematic experiment of Strykowski & Sreenivasan (1990). Other studies include sensitivity to base flow modification in the parallel Couette flow (Bottaro et al. 2003), a compressible axisymmetric body wake (Meliga et al. 2010) controlled with steady forcing in the bulk (with sources of mass, momentum or energy) and steady wall control (with blowing/suction or heating), the wake past a spheroidal bubble (Tchoufag et al. 2013), a three-dimensional T-junction (Fani et al. 2013), and a thermoacoustic system (Magri & Juniper 2013).
Because standard sensitivity analysis computes a gradient, it is linear by nature and expected to provide meaningful results in the limit of infinitesimal flow control only. For finite-amplitude control, nonlinear effects come into play, and the actual variation of the quantity of interest inevitably departs from the sensitivity prediction. This is illustrated in figure 1, which shows the effect of a localised body force on the leading growth rate λ r of the cylinder flow. At Re = 50 the uncontrolled flow is slightly unstable, λ r ( = 0) > 0. In all four control locations considered, the body force has a stabilising effect: the growth rate computed about the nonlinearly controlled base flow (symbols) initially decreases. Sensitivity analysis (dashed lines) perfectly captures the slope of the growth rate reduction at zero amplitude, dλ r /d | =0 . It does not, however, provide any information about finite amplitudes > 0: depending on the control location, sensitivity analysis is accurate up to smaller or larger amplitudes, and may or may not predict well the critical stabilising amplitude; it may also underestimate or overestimate the actual growth rate variation. This information cannot be obtained except with nonlinear calculations of the controlled flow.
Given this limitation, it is tempting to investigate whether adding one or more higherorder terms in the sensitivity analysis can improve the prediction accuracy for small but finite amplitudes. In some scientific fields, second-order sensitivity is sometimes calculated as a means to speed up the convergence of iterative gradient-based optimisation, where the modified state and the sensitivity need to be repeatedly recomputed. In hydrodynamic stability, iterative optimisation is seldom performed, and only first-order sensitivity is routinely calculated. One notable exception concerns the three-dimensional control of nominally two-dimensional (or axisymmetric) flows: when the control is periodic in the spanwise (or azimuthal) direction, the standard first-order sensitivity is exactly zero, and at leading order the effect of the control is quadratic (Hwang et al. 2013;Del Guercio et al. 2014a,b,c). In other words, expressing the eigenvalue variation with the control amplitude as λ = λ 0 + λ 1 + 2 λ 2 + . . ., the aforementioned periodic configuration is such that λ 1 = 0, and one needs to compute λ 2 . This has triggered a number of studies that either evaluated the second-order variation induced by a given control (Cossu 2014;Tammisola et al. 2014), or computed optimal spanwise-periodic flow modification or control (Tammisola 2017;Boujo et al. 2015. To the best of the author's knowledge, the second-order sensitivity of eigenvalues has never been computed in non-parallel flows subject to external control in the general case where λ 1 = 0, although the steps of the derivation are similar. Very recently, a related approach was proposed by Mensah et al. (2020) to compute second-and higher-order eigenvalue variations λ n induced by some scalar parameter modification. That approach, which explicitly computes eigenvector modifications, was applied to the parallel Poiseuille flow for variations of the Reynolds number, and to a two-dimensional time-delayed thermoacoustic system for variations of the time delay. The first aim of the present study is to propose a method for computing efficiently the second-order sensitivity of an eigenvalue with respect to control, in the context of hydrodynamic instability. Some emphasis is put on exploiting adjoint operators to derive a sensitivity that is valid for any control, instead of simply evaluating the second-order variation λ 2 for a specific control. Specifically, and postponing rigorous definitions to § 2, it might help to recall that the first-order coefficient of the eigenvalue variation can be expressed as λ 1 = ( S 1 | F), the inner product of a control F with a first-order sensitivity S 1 that depends only on the uncontrolled base flow; therefore, as S 1 is independent of the control, it can be computed once and for all, without computing controlled base flows and eigenmodes. Similarly, the present study will express the second-order variation as λ 2 = ( F | S 2 F), where the second-order sensitivity S 2 depends only on the uncontrolled base flow. The method will be illustrated with the global instability of the two-dimensional cylinder flow, controlled by a steady localised force or by a small control cylinder. The second aim of this study is to leverage second-order sensitivity to find the optimal control for stabilisation, i.e. the control yielding the largest growth rate reduction up to second order, λ 1r + 2 λ 2r .
The paper is organised as follows. Section 2 introduces the theoretical framework for the first-and second-order sensitivities of eigenvalues with respect to control ( § § 2.1-2.3). It also discusses the generalisation to higher orders ( § 2.4) and the computational cost of the method ( § 2.5). Section 3 presents the flow configuration and numerical methods. Results for the growth rate of the leading eigenmode of the cylinder flow at Re = 50 are given in § 4: sensitivity to a steady force ( § 4.1), sensitivity to a small control cylinder ( § 4.2), and an analysis of the stabilisation induced by the small control cylinder when located nearly optimally ( § 4.3). Finally, section 5 deals with optimal controls that maximise the growth rate reduction at first or second order separately, and first and second orders simultaneously. In addition, Appendix A briefly presents results for the sensitivity of the leading mode's frequency; Appendix B details the derivation of the sensitivity operators; and Appendix C outlines an extension of the method to the sensitivity of another quantity defined as an eigenvalue problem: the linear amplification of time-harmonic forcing (resolvent gain).

Base flow and stability analysis
Consider a steady fluid flow satisfying the stationary incompressible Navier-Stokes (NS) equations where P (x) is the pressure field and U(x) = (U, V ) T or (U, V, W ) T is the velocity vector in two or three dimensions. Equations are made dimensionless with a characteristic velocity U ∞ , a characteristic length scale D and the fluid kinematic viscosity ν, thus defining the Reynolds number Re = U ∞ D/ν. In the following, all velocity fields are incompressible and the continuity equation is omitted. The linear stability of the base flow is determined by the temporal evolution of small perturbations u . Considering, in particular, the normal mode ansatz u (x, t) = u(x)e λt + c.c., the (complex) eigenmodes u(x) are solutions of the linearised NS equations ( 2.3) The real and imaginary parts of an eigenvalue λ represent the linear growth rate λ r and linear frequency λ i of the associated eigenmode. The base flow is linearly unstable if at least one mode has a positive growth rate. In compact form, equations (2.1)-(2.2) for the steady base flow and (2.3) for the eigenmodes can be expressed as Here, N and A are the nonlinear and linearised Navier-Stokes operators, and I is the identity operator: (2.7)

Eigenvalue sensitivity to small-amplitude steady control
Assume now that a small-amplitude control is applied via a body force acting on the steady base flow: where ||F|| = 1 and 0 < 1. This control modifies the base flow, eigenmodes and eigenvalues, which can be expressed as power series expansions (Hinch 1991): (2.9) u = u 0 + u 1 + 2 u 2 + . . . , (2.10) Injecting the expansion (2.9) into the base flow equation (2.8) yields the following at orders 0 , 1 and 2 : Although the focus of this study is on first and second orders, note that the steady force F modifies the base flow at higher orders too, due to the nonlinear term of the NS operator. Similarly, injecting the expansions (2.9)-(2.11) into the eigenvalue problem (2.5) yields the following at orders 0 , 1 and 2 : where the operators A 1 and A 2 are linear in U 1 and U 2 , respectively, and do not depend on any other field, Before moving on to determining the first-and second-order eigenvalue variations λ 1 and λ 2 , note that the operator λ 0 I + A 0 is singular, since (2.15) holds. Therefore, according to the Fredholm alternative, commonly used in the context of weakly nonlinear expansions (see e.g. Sipp & Lebedev (2007)), (2.16)-(2.17) can be solved for u 1 and u 2 if and only if their right-hand sides have no component in the direction of the eigenmode u 0 , i.e. no projection on the adjoint mode u † 0 . Recall that the adjoint mode is a solution of where the overbar stands for complex conjugation, and A † 0 is the adjoint NS operator for the L 2 inner product ( a | b) = a T b dx for any a, b, In particular, projecting the left-hand E. Boujo side of (2.16)-(2.17) on u † 0 necessarily yields zero: Choosing the normalisation u † 0 u 0 = 1, the eigenvalue variations are obtained by projecting (2.16)-(2.17) on u † 0 (Hinch 1991;Chomaz 2005;Giannetti & Luchini 2007): (2.23) Any arbitrary component along u 0 can be added to u 1 and (2.16) will still hold, because of (2.15). This arbitrary component does not influence λ 2 , because (2.16) also implies that u † 0 (λ 1 I + A 1 )u 0 = 0. The two terms in (2.23) correspond to different mechanisms: term I is the effect of the second-order base flow modification U 2 (via the first-order flow modification and the nonlinear term of the NS operator); term II is the effect of the interaction between the first-order flow modification U 1 and the first-order eigenmode modification u 1 . As will be discussed in § 4, these two terms can either compete or collaborate.
Given a steady force F, one can compute the base flow modifications U 1 and U 2 from (2.13)-(2.14), build A 1 and A 2 , and use expressions (2.22)-(2.23) to estimate the eigenvalue variation up to first order, λ = λ 0 + λ 1 + O( 2 ), and up to second order, λ = λ 0 + λ 1 + 2 λ 2 + O( 3 ). This involves solving linear systems only, which avoids computing the nonlinear controlled flow U and solving the eigenvalue problem for the controlled mode u, thus reducing the computational cost. For instance, the dashed lines in figure 1 may be obtained by computing λ 1r this way. However, the procedure must be repeated every time a different force F is considered, which may become prohibitively expensive. More useful expressions for λ 1 and λ 2 can be obtained that do not depend explicitly on U 1 and U 2 , as explained in the next section.

Sensitivity operators
Because the operator A 1 is linear in U 1 , which itself depends linearly on F, the firstorder eigenvalue variation (2.22) can be recast as where the vector field S 1 is the usual sensitivity to a steady force (Marquet et al. 2008;Meliga et al. 2010), and depends only on the uncontrolled base flow U 0 and the uncontrolled direct and adjoint modes u 0 and u † 0 (see Appendix B). This formulation offers a significant advantage: S 1 can be calculated once and for all, and then used to predict the first-order effect of any steady force. Since no base flow modification U 1 is ever calculated, evaluating λ 1 for a large number of steady forces becomes straightforward. For instance, figure 5(a) shows the real part of the streamwise component of S 1 . The value displayed at each location x c is the first-order sensitivity of the growth rate to a steady force F = (δ(x−x c ), 0) T localised at that point and oriented along the streamwise direction.
In a similar way, because (2.23) is quadratic in U 1 and thus in F, the second-order eigenvalue variation can be recast as where S 2 is a linear operator that, again, depends only on the uncontrolled fields U 0 , u 0 and u † 0 . The derivation steps from (2.23) to (2.25) introduce suitable adjoint operators (see Appendix B), following the same steps as  for spanwise-periodic controls in nominally spanwise-invariant flows (where λ 1 = 0 and the expression for S 2 is slightly simpler). Again, this formulation suppresses the need to calculate the base flow modifications U 1 and U 2 . Once S 2 is available, λ 2 can be readily evaluated for any steady force. The dashed lines in figure 1 can now be obtained simply by probing S 1 (x c ) at each control location x c of interest.
Second-order variations are obtained just as easily, and results for the few control locations considered earlier are shown as solid lines in figure 2. The predicted growth rate variation is generally improved. In figure 2(b, d) for instance, the second-order prediction follows closely the actual growth rate variation up to much larger amplitudes than the first-order prediction. In other locations, however, like in figure 2(a, c), the improvement is less significant, owing to higher-order variations. Figure 2(e)-(h) highlights these higherorder variations, confirming their importance (figure 2e, g) or lack thereof (figure 2f, h).
As will be discussed in § 4, sensitivity maps can be produced that allow one to identify at a glance regions where a steady force alters the eigenvalue most effectively. Further, the relative signs and magnitudes of the first-and second-order eigenvalue variations will characterise the usefulness of the first-order sensitivity.
Before moving on to the next sections, it is worth mentioning that the method can be applied to the sensitivity of other quantities, as soon as they are defined as eigenvalue problems. To illustrate this point, Appendix C derives the first-and secondorder sensitivities of the linear amplification of harmonic forcing (resolvent gain).

Higher-order sensitivity
Higher-order terms U n for the base flow modification are governed by (2.28) Similarly, higher-order terms u n for the eigenmode modification are governed by which, upon projection onto u † 0 , yields the eigenvalue variations (Hinch 1991;Mensah et al. 2020): (2.34) Just like λ 1 and λ 2 are linear and quadratic in U 1 , respectively, each of the above expressions is exactly proportional to U n 1 , and thus to F n . In principle, one can therefore generalise expressions (2.24)-(2.25), which involve the vector S 1 (tensor of order one) and the matrix S 2 (tensor of order two), and introduce tensors S n of order n such that λ = λ 0 + λ 1 + 2 λ 2 + 3 λ 3 + . . . + n λ n + . . .
with Einstein notation for repeated indices. Conceptually, the method for obtaining the higher-order sensitivity operators S n is similar to that described in Appendix B, and involves a combination of the following steps: (i) redefine linear forms like A n u 0 , (λ m I + A m )u n−m , etc., so as to make explicit the dependence on the first-order flow modification U 1 , and eventually on the force F = A 0 U 1 ; (ii) introduce adjoint operators so as to isolate F, and identify the remaining control-independent operator as the sensitivity S n . It should be noted that adding more terms to the power series (2.11) does not necessarily improve its accuracy, and it certainly does not for amplitudes larger than the radius of convergence r of the expansion. In general, r depends on both the type and location of the control. In order to rigorously assess the validity of a second-order or higher-order sensitivity prediction, one must therefore compute the eigenvalue λ of the actual nonlinear controlled flow, similar to validation calculations for first-order sensitivity prediction.

Fully nonlinear
First-order Second-order controlled eigenvalue eigenvalue variation eigenvalue variation Table 1. Computational cost for the eigenvalue variation induced by a steady force, in a system discretised with N degrees of freedom and forced at M locations. The dominant contribution is derived assuming 1 M N . Recomputing the controlled base flow and the corresponding eigenvalue for each forcing location is substantially more expensive than evaluating the sensitivities.

Computational cost
The cost of computing the effect of a steady force on the eigenvalue is estimated in table 1. Different methods are compared: computing the fully nonlinear controlled flow and associated eigenvalue λ; and computing the first-and second-order eigenvalue variations λ 1 and λ 2 using sensitivity operators. In what follows, N is the total number of degrees of freedom after numerical discretisation, and M is the number of independent forcing locations. The uncontrolled base flow U 0 and leading eigenmode u 0 are computed prior to considering any control.
For the sake of simplicity, it is reasonable to assume that 1 M N when estimating the leading-order computational cost. That is, M must be rather large so as to obtain sufficiently fine-grained sensitivity maps, and N must be large enough to compute the eigenvalue and its variation accurately. To fix ideas, 10 different values for both x c and y c already yield M = 100 control locations to be evaluated. Further, with a finite-element method, a minimum of N = 10 3 to 10 4 degrees of freedom seem necessary. In this study, M 10 4 and N 6 × 10 5 .
The computational cost of the different methods is as follows.
• Recomputing the fully nonlinear controlled flow and the corresponding eigenvalue λ for each forcing location (second column of table 1) involves two steps: (i) computing M nonlinear base flows U, for instance with a Newton method requiring k linear system resolutions (typically five to ten iterations) of complexity O(N 3 ); (ii) solving M eigenvalue problems for u, for instance with an implicitly restarted Arnoldi method, of complexity proportional to O(N 2 ). Omitting constant factors for simplicity, the total cost scales like M × O(N 3 ).
• Estimating the first-and second-order eigenvalue variations with (2.24)-(2.25), i.e. with sensitivity operators (third and fourth columns of table 1), involves the following steps (see details in Appendix B): (i) computing once and for all the (uncontrolled) adjoint mode u † 0 , with a cost proportional to O(N 2 ); (ii) computing once and for all the lower-upper (LU) decompositions of complexity O(N 3 ) of A † 0 and (λ 0 I + A 0 ) for λ 1 and λ 2 , respectively; (iii) evaluating a few matrix-vector products, with a cost O(N 2 ) for each forcing location. The total cost therefore scales like O(N 3 ), for both λ 1 and λ 2 .
In conclusion, computing λ 2 involves an additional cost similar to that of computing λ 1 . It is much smaller than that of recomputing the nonlinear eigenvalue λ for each forcing location. The advantage of adjoint methods therefore applies to both first and second orders. Of course, this is true only when a large number M of control locations are considered, e.g. when constructing sensitivity maps. Conversely, when only a few control locations are of interest, calculating the actual eigenvalue λ is more accurate and not significantly more computationally expensive.
In the above analysis, memory requirements have not been considered. Storage is not an issue for two-dimensional configurations, and for spatial discretisation methods that yield sparse matrices (e.g. finite-element method), but may become prohibitive for three-dimensional configurations or methods that yield dense matrices (e.g. spectral methods). This is a practical limitation of the proposed approach. For standard eigenvalue calculations and first-order sensitivity analysis, one can use matrix-free time-stepping techniques as an alternative to matrix-based techniques (Tuckerman & Barkley (2000)). Whether such an approach is possible for second-order sensitivity analysis remains to be determined.

Flow configuration and numerical method
The two-dimensional, incompressible flow past a circular cylinder of diameter D with free-stream velocity (U ∞ , 0) T is considered. In the remainder of this study, the Reynolds number is set to Re = 50 unless otherwise stated.

Base flow
A two-dimensional triangulation of the domain is generated with the finite-element software FreeFem++ (Hecht 2012), resulting in approximately 136'000 elements. Velocity and pressure fields are discretised with P2 and P1 Taylor-Hood elements, respectively, yielding a total of N 615 000 degrees of freedom. All discrete operators are built from their continuous expressions (see details in Appendix B) in variational form. In particular, this means that the "differentiate then discretise" approach is used for adjoint operators, as opposed to the "discretise then differentiate" approach.
The uncontrolled steady base flow U = U 0 is obtained by solving (2.4) with a Newton method, iterated until residuals are smaller than 10 −12 . Boundary conditions are imposed as follows: uniform free-stream velocity at the inlet, no-slip boundary condition on the cylinder wall, outflow boundary condition −P n + Re −1 ∇U · n = 0 (with n the normal vector) at the outlet, and symmetry condition on the lateral sides of the domain. Figure 3 shows the vorticity ω = ω 0 = ∂ x V 0 − ∂ y U 0 of the base flow at Re = 50. Shear layers of opposite vorticity are created on both sides of the cylinder. The recirculation region (dashed line) extends over three diameters downstream.
Controlled base flows U are computed for validation purposes, solving (2.8) with the same method. For steady forces F that are localised in space, Dirac delta functions are smoothed out numerically into Gaussians of variance 0.0025.

Stability analysis
The eigenvalue problem (2.5) is solved with Matlab using an implicitly restarted Arnoldi method with shift-and-invert preconditioning. This study focuses on the leading eigenmode u = u 0 , which becomes unstable at Re 47 via a Hopf bifurcation, as a pair of complex conjugate eigenvalues cross the imaginary axis, as illustrated in the halfplane λ i > 0 in figure 4(a) (the other half is symmetric with respect to λ i = 0). The leading eigenmode at Re = 50, associated with the eigenvalue λ 0.0173 + 0.7797i is shown in figure 4(b). It is largest a few diameters downstream of the recirculation region, as perturbations are advected by the base flow. With its wave packet structure and its complex eigenvalue, this mode breaks both the spatial and temporal symmetries, leading to periodic vortex shedding and to the Bénard-von Kármán street in the cylinder wake.
The adjoint problem (2.19) is solved with the same method. The leading adjoint mode u † 0 shown in figure 4(c) is largest in the recirculation region, and adjoint perturbations travel upstream, a consequence of upstream advection in the adjoint NS operator.

Sensitivity
First-and second-order sensitivity maps are computed for localised control forces F. The control is moved within the subdomain x ∈ [−2, 6], y ∈ [0, 3], with a step size ∆x = ∆y = 0.05, leading to approximately M 10 000 control locations.
The second-order sensitivity operator S 2 defined by (2.25) contains inverse operators (see detailed expression in Appendix B) and is therefore not formed explicitly. Instead, the LU decomposition of each operator to be inverted is precomputed once and for all, such that each subsequent matrix inversion is replaced with two simple matrix-vector products. (Note that S 2 is a second-order tensor; by contrast, the first-order sensitivity S 1 defined by (2.24) is a vector that can be formed explicitly and plotted without further difficulty. In this study, sensitivity maps for λ 1 and λ 2 are evaluated location by location.)

Second-order sensitivity of the growth rate
This section investigates the effect of control on the first-and second-order variations of the leading growth rate λ r . (For the effect on the linear frequency λ i , see Appendix A.)

Sensitivity to a steady body force
Let us consider first a generic steady body force. Figure 5(a) shows the real part of the x component of the first-order sensitivity S 1 (x) to such a steady force. As shown by (2.24), the value at each location x = x c is also the value of the first-order variation λ 1r when choosing a localised force along the x direction, F = (δ(x − x c ), 0) T . The sensitivity is large and negative on the sides of the cylinder and inside the recirculation region, and positive on the sides of recirculation region, in agreement with Marquet et al. (2008) (figure 9(a) therein). Note that changing the sign of F x changes the sign of λ 1 , such that stabilising regions (λ 1r < 0, blue) become destabilising (λ 1r > 0, red) and vice versa.
The second-order sensitivity S 2 (x) is visualised in figure 5(b), which shows the secondorder growth rate variation λ 2r evaluated according to (2.25) for the same localised force F = (δ(x − x c ), 0) T . Overall, and in absolute value, sensitive regions are similar at first and second orders, namely the domain approximately delimited by 0 x 4, |y| 1, and containing the sides of the cylinder, the recirculation region and the shear layers. Note that, unlike S 1 , the sign of S 2 does not change with the sign of F x .
With these two maps available, it is now possible to explain the results of figure 2. The three control locations x c = (1, 0.7), x c = (1, 1) and x c = (1, 0.6) lie in a region of similar first-order sensitivity (figure 5a), and therefore induce similar first-order reductions λ 1r (figure 2a-c). The second-order variations, however, differ substantially between these three locations (figure 5b): small in x c = (1, 0.7), negative in x c = (1, 1) and positive in x c = (1, 0.6). As a result, the second-order prediction is not much different from the first-order one in figure 2(a), and yields a larger growth rate reduction in figure 2(b) and a smaller one in figure 2(c). The second-order prediction generally follows more closely the nonlinear results than the first-order one. In the last control location, x c = (3.5, 0.8), the first-order sensitivity is small (figure 5a), yielding a weak first-order variation in figure 2(d). The second-order variation, however, is clearly negative (figure 5b), and the actual growth rate reduction is well captured by the second-order prediction (figure 2d).
Considering the large differences observed between different control locations, and the potential impact on flow restabilisation, it would be useful to find a simple way to address the following questions: (i) What is the range of control amplitude where the first-order sensitivity yields an accurate prediction? (ii) Outside this range, does it underestimate or overestimate the actual variation? One step towards answering the first question is possible with the ratio of first-to second-order variations. Recalling the expansion it appears that the second-order correction 2 |λ 2 | is of the same order of magnitude as the first-order variation |λ − λ 0 | = |λ 1 | for the threshold amplitude For small enough amplitudes t , the first-order variation predicts the actual variation accurately, as the second-order correction is negligible. Conversely, for large amplitudes t , the second-order variation dominates the first-order one. In between, the secondorder variation becomes important and cannot be ignored when the control amplitude reaches some fraction of the threshold amplitude, say t /10.
Obviously, the analysis needs to be refined when λ 2 = 0. Taking into account 3 λ 3 or the next non-zero higher-order correction n λ n , the threshold amplitude becomes t = |λ 1 /λ n | 1/(n−1) . Note that the threshold amplitude decreases as the relative importance of 2 λ 2 grows; this latter term becomes the leading-order term in the limiting case λ 1 = 0 (e.g. for the spanwise-periodic control of spanwise-invariant flows), and the threshold amplitude then becomes t = |λ 2 /λ n | 1/(n−2) .
Figure 5(f ) shows the threshold amplitude (4.2), i.e. the ratio of the maps in panels (a) and (b), in logarithmic scale. Focusing on regions where λ 1r and λ 2r are not both small, it appears that the first-order prediction is especially accurate up to large amplitudes (log 10 ( t ) > 0, green) near the cylinder, downstream of the cylinder on the symmetry axis up to x = 2, and in a thin strip running along and outside the recirculation region. Conversely, the second-order prediction must be taken into account (log 10 ( t ) < −1, yellow and red) in other regions both inside and outside the recirculation region, particularly in a thin strip running along and inside it. The proximity of those two strips warns about locating a steady force near the separatrix, or in any region where t has a strong gradient: slight, unintentional shifts can dramatically increase the amplitude of the second-order variation and ruin the accuracy of the first-order prediction.
Figure 5(f ) confirms observations from figure 2: t is large and the first-order prediction is accurate over a wide range of control amplitudes in x c = (1, 1) and x c = (1, 0.7), while t is small and the second-order variation quickly becomes important in x c = (1, 0.6) and x c = (3.5, 0.8).
The second question above is answered by considering the signs of λ 1r and λ 2r . If both signs are identical, the second-order correction strengthens the effect of the first-order variation: when λ 1r , λ 2r < 0 the flow is stabilised even more than predicted by λ 1r alone (and destabilised even more when λ 1r , λ 2r > 0), such that a smaller control amplitude is actually sufficient to obtain the desired effect. Conversely, if the signs are opposite, the effect is weakened: for example, when λ 1r < 0 and λ 2r > 0, the flow is not stabilised as efficiently as predicted by λ 1r alone, such that a larger control amplitude is actually required to obtain the desired effect. As a way to distinguish between those two situations, figure 5(e) shows the sign of the product λ 1r λ 2r . Focusing again on regions where λ 1r and λ 2r are not both small, this map indicates that 'safe' regions where λ 1r λ 2r > 0 (green) are rather few and apart (mainly near the cylinder and along the separatrix), the rest being 'dangerous' regions where λ 1r λ 2r < 0 (red).
Consider again the four control configurations of figure 2, where F x < 0 (recall that the sign of λ 1r changes when the sign of F x is changed, which swaps the 'safe' and 'dangerous' regions). Figure 5(e) confirms that the first-order prediction underestimates the growth rate reduction (compared to first-and second-order predictions together) in x c = (1, 1), x c = (1, 0.7) and x c = (3.5, 0.8), and overestimates it in x c = (1, 0.6).
Let us come back to figure 5(c, d), which shows the two terms I and II in the secondorder sensitivity equation (2.23), i.e. the effects of U 2 and of the U 1 -u 1 interaction, respectively. The map in figure 5(b) is the sum of those two maps, and all three colour scales are identical. Overall, terms I and II are of the same order of magnitude. Both terms display regions of positive and negative sensitivity. They collaborate to yield positive sensitivity near the downstream end of the recirculation region, and negative sensitivity on the side of the recirculation region. Conversely, they compete on part of the symmetry axis inside the recirculation region, and on part of the separatrix, resulting in a weak total sensitivity. Although the map of term I bears an overall qualitative similarity to the map of total sensitivity, term II makes a significant contribution everywhere; in other words, the steady control force modifies the growth rate at second order by changing not only the base flow but also the eigenmode that develops on that base flow.

Sensitivity to a small control cylinder
The sensitivity analysis is now applied to a practical flow control strategy, namely inserting a small passive device in order to reduce the growth rate of the leading mode. Following Hill (1992), and later Marquet et al. (2008) and Meliga et al. (2010), the effect of a small circular cylinder of diameter d located in x c is modelled as a steady force acting on the base flow, equal and opposite to the drag force that would be felt by that cylinder in a uniform flow with the local velocity The drag coefficient C d of the control cylinder depends on the local Reynolds number  (Finn 1953;Tritton 1959) and in-house numerical simulations in the range of interest 1 Re d 15. In the following, results are illustrated with d = 0.1, i.e. a control cylinder 10 times smaller than the main cylinder.
The first-order growth rate variation induced by the control cylinder is displayed in figure 6(a). The map shows a destabilising region on the sides of the main cylinder, stabilising regions on the sides of the recirculation region, and more weakly stabilising regions on the symmetry axis both upstream and downstream of the main cylinder. This is in agreement with the map obtained by Marquet et al. (2008) (figure 11(a) therein), and is consistent with the map of figure 5(a), since the force (4.3) is oriented mainly along −x outside the recirculation region and mainly along x inside. Figure 6(b) shows the second-order growth rate variation. The main destabilising and stabilising regions appear rather similar to those of the first-order variation of panel (a). This means that, where the signs of those regions do correspond, the second-order variation tends to strengthen the effect of the first-order one. A closer look at figure 6(e) reveals that, where λ 1r and λ 2r are not both small, they generally have the same sign. Therefore, for a small control cylinder, and considering the variation of λ r up to second order, the situation is one of the following almost everywhere: (i) both λ 1r and λ 2r are small, so the control cylinder does not modify the growth rate substantially; (ii) only λ 2r is small, so the effect of the control cylinder is well predicted by λ 1r alone; (iii) λ 1r and λ 2r are not small and have the same sign, so the effect of the control cylinder is stronger (more destabilising or more stabilising) than predicted by λ 1r alone. One exception is the narrow region where λ 2r > 0 and λ 1r is small: although first-order sensitivity predicts no effect, the control cylinder is actually destabilising. The decomposition of λ 2r into terms I and II in figure 6(c)-(d) shows that the secondorder destabilising effect is primarily due to U 2 , while the second-order stabilising effect is due both to U 2 and to the U 1 -u 1 interaction. Figure 7(a)-(b) shows the contours where inserting a small control cylinder of diameter d = 0.1, as described above, is predicted to make the leading mode neutrally stable. Several Reynolds numbers Re 50 are considered. Inside the regions delimited by these contours, the mode is stable and vortex shedding is expected to be suppressed. In figure 7(a), only the first-order sensitivity prediction is considered, λ 0r + λ 1r = 0, while in figure 7(b), the second-order correction is included too, λ 0r + λ 1r + 2 λ 2r = 0. The results compare qualitatively well with the experimental observations of Strykowski & Sreenivasan (1990) (figure 20 therein): stabilisation is achieved on the side of the recirculation region, in an area that is rather wide at Re = 50 and that becomes smaller as the Reynolds number increases, until shrinking to a single point and vanishing when restabilisation is not possible any more. Compared to the first-order sensitivity, however, the second-order sensitivity seems to better capture the results of Strykowski & Sreenivasan (1990): in particular, it predicts a wider stabilising area at Re = 60, and a larger value of the maximum stabilisable Reynolds number 70 < Re < 80.
For completeness, figure 7(c)-(d) shows stabilising contours for a pair of control cylinders located symmetrically in (x c , y c ) and (x c , −y c ), still with d = 0.1. In the sensitivity framework, the two cylinders are assumed not to influence each other, which is not satisfied close to the symmetry axis y = 0. Unsurprisingly, the main stabilising region is wider but still located on the side of the recirculation region. Although conclusions should be drawn with care at larger Reynolds numbers, as the uncontrolled flow becomes linearly unstable to a second two-dimensional mode at Re 100 (Verma & Mittal 2011) and to a three-dimensional mode at Re 190 (Barkley & Henderson 1996), restabilisation can be achieved up to Re 100 and Re > 100 according to first-and second-order sensitivity predictions, respectively.

Analysis of the stabilisation induced by a small control cylinder located optimally
In this section, the effect of a small control cylinder is investigated in more detail for the specific control location x c = (1, 1), close to the location of largest first-and second-order stabilising effects identified in § 4.2. Figure 8 shows the eigenspectrum of the flow controlled with a secondary cylinder of increasing diameter d. The leading mode is restabilised for diameters d 0.004. Other modes remain stable for the whole range of diameters investigated. As seen in the close-up (figure 8b), the second-order sensitivity (thick solid line) follows closely the actual path of the leading eigenvalue in the complex plane (symbols), accurately capturing both the growth rate and the frequency, and improving on the first-order prediction (dashed line).
Let us now focus on the diameter d = 0.1. As apparent from figure 6, the predicted first-and second-order growth rate variations are comparable, (4.4) and the two second-order contributions are of similar order of magnitude: 2 λ 2r,I = −0.0258, 2 λ 2r,II = −0.0167. Figure 9 depicts the base flow modification. At first order, the control cylinder induces a strong velocity deficit U 1 < 0 in its wake, and a slight acceleration U 1 > 0 between the two cylinders (figure 9a). As a result, two layers of opposite vorticity emanate from the control cylinder (figure 9b) in a roughly symmetric way. At second order, velocity and vorticity are modified more weakly, with a more complicated spatial pattern ( figure 9c, d). The net effect of the control cylinder is best illustrated by the velocity and shear profiles in figure 9(e). In x = 0.8, just upstream of the control location x c = 1, the flow modification U 1 (red dashed line) smooths the velocity profile and reduces the maximum shear. Downstream (x = 1.5 and 2), the induced velocity deficit further reduces shear in the lower shear layer emanating from the control cylinder, where the positive vorticity ω 1 (figure 9b) counteracts the negative base flow vorticity ω 0 (figure 3). These observations are consistent with those of Marquet et al. (2008). In the upper shear layer, however, the negative ω 1 adds up to the negative ω 0 , and shear is strongly increased, well beyond the maximum uncontrolled shear. The second-order modification (green solid line) tends to yield an additional reduction in maximum shear, both upstream and downstream of x c , albeit much smaller. In light of these observations, shear alone does not seem to explain entirely (i) why U 1 is stabilising, and (ii) why U 2 brings an additional stabilisation as large as U 1 .
Some complementary insight can be gained by looking at regions that contribute to the growth rate variation. Recalling that λ 1 and λ 2 are defined by (2.22)-(2.23) as inner products, it is natural to look at the integrands of λ 1 , λ 2,I and λ 2,II , shown in figure 10(a, c, e). For a more quantitative picture, it is useful to consider the contribution from each streamwise location x: let us integrate those integrands vertically and define one-dimensional densities, By construction, the cumulative integral x −∞ l 1 (x )dx tends to λ 1r as x → ∞. Similarly, the limits of the cumulative integrals of l 2,I (x) and l 2,II (x) are λ 2r,I and λ 2r,II , respectively. These densities and cumulative integrals are shown in figure 10(b, d, f ) as dash-dotted lines and solid lines, respectively. All three densities are positive at the control cylinder location; farther downstream they become negative, in a longer region and with a similar intensity, finally resulting in λ 1r < 0 and λ 2r < 0. The two-dimensional integrands are mostly positive in the early wake of the control cylinder, and negative in a wider region running downstream along the separatrix. This region can therefore be identified as the main stabilising one when λ 1r and λ 2r are understood as inner products expressed in terms of modifications A 1 , A 2 of the linearised NS operator, and eigenmode modification u 1 . It is also possible, and perhaps more informative, to consider alternative expressions for λ 1r and λ 2r where the base flow modification U 1 appears explicitly. The interested reader is referred to (B 5) and (B 18) in Appendix B, where the sensitivity operators are  x −∞ l2,II(x ) dx → λ2r,II Figure 10. Integrands of the first-and second-order growth rate variations λ1r and λ2r expressed as (2.22)-(2.23), for a small control cylinder (d = 0.1, xc = (1, 1)) at Re = 50. (a) Integrand of (4.5); (b) density l1(x) (black dash-dotted line) and its cumulative integral (red solid line). (c) Integrand of (4.6); (d) density l2,I(x) (black dash-dotted line) and its cumulative integral (green solid line). (e) Integrand of (4.7); (f ) density l2,II(x) (black dash-dotted line) and its cumulative integral (green solid line).
derived. The corresponding integrands are shown in figure 11(a, c, e), and the densities (4.10) in figure 11(b, d, f ). The density l 1 (x) is qualitatively similar to l 1 (x): positive around x c and negative in a longer region downstream. The two-dimensional integrand, however, exhibits a more complicated structure with alternating positive regions (separatrix and control cylinder wake) and negative regions (especially the recirculation region). This reveals that the main first-order stabilising contribution in terms of flow modification  x −∞ l 2,II (x ) dx → λ2r,II Figure 11. Same as figure 10, for λ1r and λ2r expressed with sensitivities to first-order base flow modification U1. (a) Integrand of (4.8); (b) density l 1 (x) (black dash-dotted line) and its cumulative integral (red solid line). (c) Integrand of (4.9); (d) density l 2,I (x) (black dash-dotted line) and its cumulative integral (green solid line). (e) Integrand of (4.10); (f ) density l 2,II (x) (black dash-dotted line) and its cumulative integral (green solid line).
U 1 comes from the inside the recirculation region, not directly from the control cylinder wake. Again, this is consistent with the observations of Marquet et al. (2008). Turning now to second order, it appears that l 2,I and l 2,II are mostly negative or zero, and only marginally positive. The integrand of λ 2r,I is strongly negative immediately upstream of the control cylinder, while the integrand of λ 2r,II is negative in the control cylinder wake and along the separatrix. Therefore, the main second-order stabilising contribution from the flow modification (quadratic effect of U 1 ) comes directly from the control cylinder and its wake.

Optimal control
Previous sections have investigated the second-order sensitivity to given, localised controls. One may wonder about how to design an optimal distributed control F opt so as to maximise the growth rate reduction. This section first recalls how to compute optimal controls targeting separately the first-and second-order variations, λ 1r and λ 2r , and then presents a method for computing the optimal control targeting at once the total second-order variation, λ 1r + 2 λ 2r . This method, borrowed from general linear algebra and applied mathematics, seems rather new in the field of hydrodynamic stability.

Optimising first-and second-order variations separately
When only the first-order variation λ 1r is considered, the optimal unit control is proportional to the sensitivity itself (Bottaro et al. 2003;Boujo et al. 2015): This classical result can be obtained with a Lagrangian method, or as a direct consequence of the Cauchy-Schwarz inequality becoming an equality for two linearly dependent vectors.
When only the second-order variation λ 2r is considered (which is relevant when λ 1r = 0, for example for the spanwise-periodic control of spanwise-invariant flows), the largest growth rate reduction is i.e. the optimal unit control F opt 2 is the eigenvector associated with the smallest eigenvalue µ of the following symmetric eigenvalue problem (Boujo et al. 2015:

Optimising the total second-order variation
If now the total second-order variation is to be minimised, one can introduce the Lagrangian where β is an as yet unknown Lagrange multiplier enforcing the normalisation ||F|| = 1. From the stationarity condition ∂L/∂F = 0, one obtains the following equation for the optimal unit control F opt 1+2 : One can verify that: (i) in the limit of small control amplitudes, 1, the optimal control reduces to F opt 1 proportional to S 1r , as given by (5.1); (ii) in the limit of vanishing firstorder sensitivity, S 1r = 0, the optimal control reduces to the F opt 2 solution of an eigenvalue problem equivalent to (5.3). In both cases, F opt 1 and F opt 2 are independent of the control amplitude .
In general, however, equation (5.6) for F opt 1+2 is neither a linear system nor an eigenvalue problem, and F opt 1+2 depends on the amplitude considered. Together with the associated constrained minimisation problem (5.4), it appears in some least-squares problems, constrained eigenvalue problems and trust-region problems. It has been studied extensively in the literature, and several solution techniques are available. For instance, Gander et al. (1989) give an iterative method based on solving a so-called explicit secular equation for β, but it involves a full diagonalisation of the operator S 2r + S T 2r , which is not tractable in the present study. Another method consists in finding the smallest β solution of the implicit secular equation In either case, the optimal control F opt 1+2 is obtained by substituting the obtained value of β in (5.6).
Here, yet another approach from Gander et al. (1989) is used. First, F is expressed from (5.6) as F = − 2 S 2r + S T 2r − 2βI −1 S 1r , and the (unit) norm of F becomes because the operator in square brackets is symmetric. Second, defining the vector (5.10) finally obtaining the quadratic eigenvalue problem (5.11) to be solved for the smallest eigenvalue β. The associated eigenvector γ yields the optimal control F opt 1+2 via (5.9). In practice, the quadratic eigenvalue problem is transformed into an equivalent linear one, (5.12) which has twice the dimension but can be solved with standard methods. As mentioned earlier, the optimal control F opt 1+2 is a function of the amplitude considered because is a parameter of the optimisation problem (5.12), which one is free to choose. In the following, let us denote by * the optimisation amplitude. A given optimisation amplitude * yields an optimal unit control F opt 1+2 , and the associated values λ 1r and λ 2r . The control * F opt 1+2 is therefore optimal for this amplitude. When implementing this optimal unit control with another amplitude = * , the second-order effect of F opt 1+2 will be λ r = λ 0r + λ 1r + 2 λ 2r . By construction, this effect will be optimal only for = * . Figure 12(a, b) compares the linear variation of the leading growth obtained with the first-order optimal control F opt 1 (dashed line), and the quadratic variation obtained with the total second-order optimal control F opt 1+2 (solid lines) for several optimisation amplitudes * (symbols). In all cases, the second-order effect is stabilising (λ 2r < 0). For Re = 50 (figure 12a), changing * makes little difference. For Re = 80 (figure 12b), however, the impact of * is clearly visible: controls optimised for larger amplitudes * perform better at large , but worse at small (see inset). This highlights the flexibility of the method, which allows one to select a control amplitude and optimise for that specific amplitude. Figure 12(c)-(h) shows the unit optimal control for several values of * , at Re = 50 (left column) and Re = 80 (right column). Each panel compares the first-order optimal control F opt 1 (upper half), and the total second-order optimal control F opt 1+2 (lower half). Contours show the magnitude of the vector F, streamlines show its orientation. For small  Figure 12. Optimisation of the total second-order variation. (a, b) Quadratic variation of the leading growth rate λ0r + λ1r + 2 λ2r induced by the optimal control F opt 1+2 . Each solid line corresponds to a different optimisation amplitude * (symbols). Dashed line: linear variation for the first-order optimal F opt 1 . Inset: close-up of the small-amplitude region, also showing the linear variations (slopes in = 0). (c)-(h) Optimal unit control for first-order growth rate variation only ( * = 0, upper half) and for total first-and second-order growth rate variation ( * > 0, lower half). Colour, magnitude; streamlines, local orientation. Optimisation amplitude: (c, d) * = 0.02, (e, f ) * = 0.1, and (g, h) * = 0.5. Reynolds number: (a, c, e, g) Re = 50, and (b, d, f, h) Re = 80.

E. Boujo
amplitudes, F opt 1+2 is very similar to F opt 1 , as seen in figure 12(c, d) for * = 0.02. As the optimisation amplitude increases ( * = 0.1 in figure 12e,f and * = 0.5 in figure 12g,h), the optimal control becomes weaker in the recirculation region and immediately outside, and stronger in a new area outside the recirculation region, while its overall orientation is preserved. For finite control amplitudes, making small changes to a control may thus be important, as this can improve its second-order effect.

Conclusion
A second-order sensitivity operator has been derived and used to predict quadratic eigenvalue variations induced by flow control. Introducing suitable adjoint operators, this second-order sensitivity is made independent of the control. First-and second-order sensitivity maps have been obtained for the control of the cylinder wake with a steady body force and a model of a small control cylinder, at a much lower computational cost than by recomputing nonlinear controlled flows and eigenmodes. Considering finiteamplitude control, the range of validity of the first-order sensitivity is characterised with a map of 'threshold amplitude'. Regions where the first-order sensitivity underestimates or overestimates the eigenvalue variation up to second order are also conveniently visualised with another dedicated map. The effect of a small control cylinder tends to be underestimated, such that regions where the flow is fully restabilised become larger when including second-order effects at all the Reynolds numbers investigated. Decomposing the second-order variation into two contributions (second-order base flow modification, and interaction between first-order base flow and eigenmode modifications, respectively) reveals that both contributions are equally important in the most sensitive regions. Analysing the effect of a small control cylinder located nearly optimally shows that stabilising effects arise from flow modifications in different regions: inside the recirculation region for first-order stabilisation, immediately upstream of the control cylinder and in its wake for second-order stabilisation.
Finally, with the second-order sensitivity operator available, the optimal control (distributed body force) for stabilisation up to second order is computed. While the first-order optimal control is directly proportional to the first-order sensitivity (and independent of the control amplitude), the total second-order optimal control is obtained via a quadratic eigenvalue problem and depends on the amplitude. As the amplitude increases, this control becomes stronger on the sides of the cylinder and the recirculation region, and weaker inside the recirculation region. Therefore, given a desired amplitude, it is possible to fine-tune the control.
While first-order sensitivity perfectly captures the effect of infinitesimal control on linear stability properties, this study shows that adjoint-based second-order sensitivity provides a range of useful information for finite-amplitude control, at little extra computational cost. At some locations in the cylinder flow, e.g. in the shear layers, it seems that higher-order terms would improve the sensitivity prediction. This has not been investigated systematically in the present study, so several questions remain open, including the following ones: In which regions are higher-order effects λ n stronger? How does the radius of convergence r of the power expansion vary in space? Is it possible to relate spatial distributions of λ n and r to any physical mechanism, in this and other flows?
The present approach can easily be extended to other types of control such as wall blowing/suction and shape deformation. It is expected to be useful for the passive control of other globally unstable flows, and may be applied to stable flows too since the resolvent gain (amplification of time-harmonic perturbations) can be expressed as an eigenvalue problem and treated in a similar framework. It could also be used to speed up the convergence of gradient-based optimisation when iteratively designing practical controls aiming for flow stabilisation.

Declaration of Interests
The author reports no conflict of interest. (ii) Term I is dominant on the sides of the cylinder, whereas term II is dominant on the sides of both the cylinder and the recirculation region.
(iii) The second-order effect is approximately one order of magnitude smaller than the first-order effect. This contrasts with the growth rate (first-and second-order effects of the same order of magnitude; see figure 6).

B.1. First-order sensitivity operator
Recall the first-order eigenvalue variation (2.22) induced by a steady force F: Next, define the linear operator L, which depends only on u 0 , such that Substituting into (B 1) yields where the adjoint operator of L reads The first-order sensitivity to base flow modification is therefore and since U 1 is a solution of (2.13), the first-order sensitivity to a steady force F is One recognises the usual sensitivities to flow modification and to a steady force (Marquet et al. 2008;Meliga et al. 2010).

B.2. Second-order sensitivity operator
Recall the second-order eigenvalue variation (2.23) induced by a steady force F: As explained in § 2.2, u 1 is defined up to any component along u 0 (i.e. u 1 = u 1 + αu 0 such that u 1 has no component along u 0 ). Injecting and developing yields so the arbitrary component αu 0 does not modify λ 2 . The second term can be rewritten in terms of u 0 by recalling that u 1 is a solution of (2.16): Next, defining the linear operators T and M, which depend only on u 0 and u † 0 , respectively, such that (λ 1 I + A 1 )u 0 = λ 1 u 0 + U 1 · ∇u 0 + u 0 · ∇U 1 = TU 1 , (B 10) substituting into (B 7) and noting that A 2 u 0 = LU 2 yields where the adjoint operator of M reads Since U 2 is a solution of (2.14), one can rewrite the first term, where U † is a solution of Finally, introducing the linear operator E. Boujo allows one to rearrange the first term: The second-order sensitivity to base flow modification is therefore and since U 1 is a solution of (2.13), the second-order sensitivity to a steady force F is Like in (2.23), term I is the effect of U 2 and term II is the effect of the U 1 -u 1 interaction, the base flow, the linear response and the resolvent gain are modified and can be expressed as power series expansions, U = U 0 + U 1 + 2 U 2 + . . . , (C 8) u = u 0 + u 1 + 2 u 2 + . . . , (C 9) G 2 = G 2 0 + G 2 1 + 2 G 2 2 + . . . , (C 10) and one would like to predict the first ans second-order gain variations G 1 and G 2 . Two cases should be distinguished: either (i) the harmonic forcing f is prescribed, or (ii) the optimal gain is of interest and the optimal forcing f opt is itself modified by the control as f opt = f 0 + f 1 + 2 f 2 + . . . .

(C 11)
Let us consider for now the most general case (ii). Injecting the above expansions in (C 6)-(C 7) yields equations (2.12)-(2.14) for U 0 , U 1 and U 2 , and the following equations analogous to (2.15)-(2.17) for the response: In the derivation of the above equations, the expansion R = R 0 + R 1 + 2 R 2 + . . . has been injected into R † R, giving and the expansion A = A 0 + A 1 + 2 A 2 + . . . has been injected into R = (iωI + A) −1 , allowing one to identify (C 20) Projecting (C 13)-(C 14) on the adjoint forcing f † = f (note that R † R is self-adjoint) and choosing the normalisation ( f 0 | f 0 ) = 1 yields the expressions of the desired gain variations, similar to (2.22)-(2.23): For a given control F, one can easily compute the base flow modifications U 1 and U 2 , build the operators A 1 , A 2 , R 1 and R 2 , compute the forcing modification f 1 , and finally calculate the first-and second-order gain variations G 2 1 and G 2 2 . More interestingly, it is possible to recast these variations as G 2 1 = ( S 1 | F) , (C 23) where the sensitivity operators S 1 and S 2 depend only on the uncontrolled base flow U 0 and the forcing f 0 . The derivation involves introducing suitable adjoint operators, along the same lines as the derivation of the sensitivity operators for λ 1 and λ 2 . The final result reads for the first-order sensitivity, where one recognises the usual sensitivity to a steady force , and for the second-order sensitivity, where K, L, M and T are now defined by Comparing with the second-order eigenvalue sensitivity (B 19), it appears that term II (from the U 1 -u 1 interaction) is directly analogous, while term I (from U 2 ) contains an analogous part depending on K but also an additional part.
Coming back to case (i), where the harmonic forcing f is fixed, the second-order gain variation becomes G 2 2 = f | (R † R) 2 f , term II is null, and the second-order sensitivity operator reduces to