Low-frequency resolvent analysis of the laminar oblique shock wave/boundary layer interaction

Abstract Resolvent analysis is used to study the low-frequency behaviour of the laminar oblique shock wave/boundary layer interaction (SWBLI). It is shown that the computed optimal gain, which can be seen as a transfer function of the system, follows a first-order low-pass filter equation, recovering the results of Touber & Sandham (J. Fluid Mech., vol. 671, 2011, pp. 417–465). This behaviour is understood as proceeding from the excitation of a single stable, steady global mode whose damping rate sets the time scale of the filter. Different Mach and Reynolds numbers are studied, covering different recirculation lengths $L$. This damping rate is found to scale as $1/L$, leading to a constant Strouhal number $St_{L}$ as observed in the literature. It is associated with a breathing motion of the recirculation bubble. This analysis furthermore supports the idea that the low-frequency dynamics of the SWBLI is a forced dynamics, in which background perturbations continuously excite the flow. The investigation is then carried out for three-dimensional perturbations for which two regimes are identified. At low wavenumbers of the order of $L$, a modal mechanism similar to that of two-dimensional perturbations is found and exhibits larger values of the optimal gain. At larger wavenumbers, of the order of the boundary layer thickness, the growth of streaks, which results from a non-modal mechanism, is detected. No interaction with the recirculation region is observed. Based on these results, the potential prevalence of three-dimensional effects in the low-frequency dynamics of the SWBLI is discussed.


Introduction
An accurate understanding of the dynamics of the shock wave / boundary layer interaction (SWBLI) is needed in many aerospace and aeronautical applications to predict, for example, flows around transonic airfoils, supersonic air intakes, or deflected control surfaces of vehicles at transonic or supersonic speed.These interactions can lead to an increase of drag, to separation and to loss of performances.Moreover, shock wave / turbulent boundary layer interaction (SWTBLI) generally produces a low-frequency unsteadiness of the shock system.This unsteadiness can for example modify the thermal load or induce fatigue of the structure.It is well known that the flow separation may be at the origin of this unsteadiness in any regime, from the incompressible (Weiss et al. 2015;Mohammed-Taifour & Weiss 2016;Le Floc'h et al. 2020) to the hypersonic (Gaitonde 2015; Priebe & Martín 2021), taking a particular form for compressible flows when the separation is induced by a shock wave.A considerable amount of work has been carried out to investigate the steady and unsteady aspects (see Gaitonde (2015) for a review).Experimental research of the SWBLI started in the mid-1940s with the work of Ackeret et al. (1947) and has shown a continuing interest since then.At this time, most of the experiments only measure steady aerodynamic quantities such as pressure distribution, skin friction, and heat transfer rates.Detailed investigations of the phenomenon and its dependence on flow and boundary layer parameters have been already published Adamson & Messiter (1980); Délery & Marvin (1986); Stanewsky (1988).
It was only at the end of the 90s that the study of the dynamics of a SWTBLI was more systematically undertaken (Dolling 2001;Knight et al. 2003).It has been shown that the dynamics is rich, spanning over several frequency decades.Two main length scales are at play: the vorticity thickness of the shear layer coming from the separated region and linked to the disturbances developing along the shear layer and, on a more global scale, the separation length  that is related to lower frequencies and is associated with the dynamics of the separation bubble (Dupont 2012;Gaitonde 2015).Three main frequency ranges can be identified.First, high frequencies, characterised by Strouhal numbers based on  typically in the range of   1, are associated with the transitional or turbulent boundary layer dynamics (Bonne et al. 2019).Secondly, a vortex shedding resulting from the shear layer dynamics usually takes place at mid-frequency (0.1   1).This can be observed either in incompressible recirculation bubble (Weiss et al. 2015;Mohammed-Taifour & Weiss 2016) or in the SWBLI (Thomas et al. 1994;Agostini et al. 2012).Finally, a third dynamics is observed in the low-frequency range (  < 0.1) (Clemens & Narayanaswamy 2014).This is the frequency domain we will mainly focus on from now on.
The low-frequency range has been intensely studied since the early 2000s in order to unveil its physical origin (Dussauge et al. 2006).Dupont et al. (2005Dupont et al. ( , 2006) ) have shown that, for an interaction strong enough, low-frequency unsteadiness is observed around   0.03 and is characterised by oscillations of the separated shock (Dussauge et al. 2006;Humble et al. 2009).Scenarios attempting to model this low-frequency dynamics have been devised.The first of them is related to the interaction of large structures upstream of the boundary layer with the shock leading to low-frequency response of the shock (Ganapathisubramani et al. 2007).However, Dussauge & Piponniau (2008) have shown that the influence of downstream conditions, especially in the recirculation zone, is more significant than the upstream conditions with respect to low-frequencies when the interaction is strong.Piponniau et al. (2009) proposed a model based on a mass balance of the system "shear layerseparated zone" where coherent structures in the shear layer feed the recirculation zone which increases up to a critical size beyond which it empties, causing a breathing of the bubble and consequently the movement of separation of the shock (Agostini et al. 2015).Touber & Sandham (2011) furthermore characterised the system using a combination of theoretical and numerical model.They showed that the separated shock foot acts as a low-pass filter with respect to white noise.This result is in agreement with the linear interaction analysis (Ribner 1953;Robinet & Casalis 2001) of the shock response to a harmonic infinitesimal perturbation.
The three-dimensional dynamics resulting from the development of 3D disturbances, independently of aspects strictly related to the effects of turbulence, has been studied more recently.Mainly three types of phenomena have been identified.When the upstream boundary layer is turbulent and the SWBLI weak (incipient separated zone), streaks from the upstream boundary layer and associated with longitudinal vortices can interact with the shock system (Ganapathisubramani et al. 2007(Ganapathisubramani et al. , 2009;;Di Renzo et al. 2021).However, when the SWBLI is strong enough to generate a significant separated zone, the latter can induce a streamline curvature of the boundary layer and initiate a centrifugal Görtler-like instability (Pasquariello et al. 2017;Zhuang et al. 2018).When the SWBLI is sufficiently strong, the separated zone becomes three-dimensional and the emergence of a global non-oscillating instability can be observed.This characteristic is very general and can be observed from subsonic (Theofilis et al. 2000;Rodriguez & Theofilis 2010;Rodríguez et al. 2021) to supersonic regime (Robinet 2007;Hildebrand et al. 2018).The role that these 3D effects may play in the low-frequency unsteadiness remains however unclear.
From a numerical point of view, a substantial body of work has been done over the years, including direct numerical simulations and large eddy simulations.Touber & Sandham (2009a,b); Pirozzoli & Bernardini (2011)'s works made it possible for the first time to characterise the dynamics with reasonable statistical convergence, showing a very good agreement with IUSTI's experiments (Dupont et al. 2006(Dupont et al. , 2008)).In particular, these numerical studies also showed the broadband nature of the low-frequency dynamics.Priebe et al. (2009) for a reflected shock and Wu & Martin (2008) for a compression ramp detailed the low-frequency dynamics and recovered the scaling proposed in the literature.They also observed a large amplification of the turbulence through the SWTBLI.Priebe & Martin (2012) for the compression ramp have studied the physical mechanism that drives the shock motion.In their simulations, the flow undergoes low-frequency changes of topology in the interaction region, including the breaking-up of the recirculation bubble and the shedding of vortical structures.In addition, the growth of energetic turbulent structures in the shear layer was found to be modulated at low frequency.This could imply a modulation of the shear layer entrainment rate which is consistent with the scenario of Piponniau et al. (2009).These authors thus suggest that the low frequency dynamics is related to the dynamics of the separated zone.Aubard et al. (2013) and Priebe & Martin (2012) then showed that the dynamics of the separated zone occurs at medium frequency and is linked to Kelvin-Helmholtz instabilities in the shear layer and that the low-frequency dynamics is only a modulation of the latter and corresponds to the breathing of the separated zone.
Resolvent analysis, or resolvent-based modelling, has been intensively used in the fluid dynamics community in the recent years.This approach is based on the singular value decomposition (SVD) of the resolvent operator that provides an optimal orthonormal basis for an external forcing (input) of the system and its associated linear response (output).The optimality is here defined as the way for a given energy input, at a given frequency, to force the system in order to trigger the largest possible energy growth.This is therefore particularly suited to study convective instabilities in noise amplifier flows.In such flows, ranges of frequency of external perturbations are amplified in space and time through linear mechanisms, while others are damped (Huerre & Monkewitz 1990).Resolvent analysis was successfully applied, for example, to boundary layer (Åkervik et al. 2008;Sipp & Marquet 2013) and jet (Nichols & Lele 2011;Garnaud et al. 2013) flows.Because nonmodal effects are taken into account (Schmid 2007), algebraic instabilities like streaks or non-modal mechanisms such as the Orr mechanism can furthermore be modelled via this approach (Monokrousos et al. 2010).Resolvent analysis has also been intensively used in the modelling of turbulent flow since McKeon & Sharma (2010) noticed it could be used to identify coherent, large scale structures.While turbulence is a highly non-linear phenomenon, the (linear) resolvent analysis yet captures the growth of these structures.The interpretation, and consequently the relevance of the resolvent analysis, is that non-linear forcing terms, whose knowledge is not ultimately required here, provide background perturbations that grow into energetic coherent structures through linear mechanisms.Several authors later discussed the potential of the resolvent approach in reduced-order modelling of turbulent flows (Semeraro et al. 2016), noticing, in particular, the link between resolvent and spectral-POD modes (Towne et al. 2018).
Resolvent analysis has been recently used to study the dynamics of the SWBLI.(Sartor et al. 2015) used a turbulent mean flow of a transonic SWBLI over a bump and compared it to their experimental results.They obtained similar features at low and medium frequencies, supporting the idea of a forced dynamics as their mean flow was found globally stable.(Bonne et al. 2019) studied a transitional case, in which they identified three ranges of frequencies and proposed a scenario for the low-frequency dynamics.In our paper, we focus exclusively on the low-frequency dynamics of the SWBLI, that we analysed by means of the resolvent analysis.Unlike the two papers cited above, a laminar base flow will be considered, for which similar features of the low-frequency dynamics can be observed (Yao & Sandham 2002;Boin et al. 2006;Sansica et al. 2014;Sandham et al. 2014;Sansica et al. 2016).Our objective is first to understand to what extent a resolvent-based approach can model the behaviour of the SWBLI for this range of frequency, where the dynamics does not result from convective or global instabilities.This will also provide insights on the physics of the system that will be discussed.In particular, by conducting our analysis for different Mach and Reynolds numbers, we aim at characterising the scales involved in this dynamics.Besides, despite considering a laminar flow, this analysis could, to some extent, support the idea that the turbulent SWBLI has a forced dynamics at low frequency.Indeed, if the intrinsic forcing emanating from non-linear interactions in the turbulent case reasonably projects onto the optimal forcing found in resolvent analysis, it is then possible to assume that a similar response will be found.All of the above will be carried out considering 2D perturbations only, similarly to the two previously cited papers.The second objective is to explore the resolvent analysis of the SWBLI for 3D perturbations for which, to our best knowledge, no work has been published.By 3D perturbations, it is meant that Fourier modes will be considered in the homogeneous spanwise direction, characterised by their wave number .The aim is to investigate the role that streaks or any 3D effects may play in the low-frequency dynamics of the system.
The paper is structured as follows.The governing equations and the theoretical and numerical approaches are described in § 2. The methodology to generate the set of base flows understudied is then detailed in § 3.In section § 4, the resolvent analysis is carried out for 2D perturbations on different base flows.A link with global stability analysis will be proposed, leading to a model of the optimal and sub-optimal gain obtained from the SVD of the resolvent operator.The scaling of key quantities will be shown, before discussing these results with respect to the literature.In section § 5, 3D perturbations will be considered in the resolvent analysis.Only one base flow will be studied, for which the newly obtained features of the dynamics will be classified into two regimes depending on their wave numbers.This section ends with a discussion on the 3D effects.A general conclusion can be found in § 6.

Governing equations
The compressible non-linear Navier-Stokes equations are considered for the conservative variables q = (, u, )  , where , u and  are the density, the velocity vector and the specific total energy, respectively.The components of the velocity vector are noted ,  and  in the streamwise direction , the wall-normal direction  and the spanwise direction , respectively.The variables , , , and , standing for temperature, pressure, dynamic and kinematic viscosity and thermal conductivity, are also introduced.All variables are non-dimensional using their value at infinity, except the pressure that is normalised by  ∞  2 ∞ .The reference length scale and time scale are  * 0 and  * 0 / ∞ , where  * 0 is the compressible displacement thickness at a chosen streamwise location , defined as More details on the choice of  * 0 for our problem will be given in section 3.2.The governing equations read The fluid is assumed to be a perfect gas whose equation of states and total energy read Furthermore assuming the fluid to be Newtonian, the viscous stress tensor verifies Sutherland's law is used for the dynamic viscosity where   = 110.4K and  ∞ = 288.The thermal conductivity  is also assumed to follow Sutherland's law (Toro 2013).Three non-dimensional numbers have been introduced that control the flow: ,  and , the Reynolds, Mach and Prandtl numbers respectively, defined as where  ∞ is the speed of sound at infinity and   the heat capacity.While  is set to  = 0.72 in order to model the fluid as air, the influence of  and  will be studied.Note that a geometrical parameter  will also be introduced.It is defined as the angle between the streamwise direction and the shock wave impinging the boundary layer.Finally, equations (2.1) will be recast in the dynamical system form as where N is the non-linear differential operator of the Navier-Stokes equations.

Global stability
In this paper, the shock wave boundary layer interaction is studied via two types of linear analysis.The linearised Navier-Stokes equations, considered in their semi-discrete form, are written as q where J and q are the Jacobian matrix and the state-vector perturbations around the base flow q, respectively.The base flow q(, ) corresponds to a steady 2D laminar solution of the non-linear equations (2.1).A Fourier transform in both time and space (in the spanwise direction ) can then be performed, allowing the perturbations to be studied as the following global modes q (, , , ) = q(, )  (−) + .. (2.8) with  the angular frequency and  and the spanwise wave number.Considering  ∈ C and  ∈ R, the global stability (Theofilis 2011) of the system can be studied by solving the eigenvalue problem If an eigenvalue is found such that the imaginary part of  is positive, then the flow is globally unstable.In this paper, only globally stable flows will be considered: if not externally sustained, any linear perturbations will eventually decay.

Resolvent analysis
Even though globally stable, flows such as boundary layer or mixing layer can exhibit convective instabilities (Huerre & Monkewitz 1990).In such flows, termed as noise amplifier flows, studying the linear forced dynamics is relevant to understand and model their behaviour.This can be achieved using resolvent analysis (Sipp et al. 2010) which has been shown to provide a good model for convective instabilities over ranges of frequency where large energy growths take place (Sipp & Marquet 2013).Besides, this approach gives the location and structure of the optimal harmonic forcing of the flow, i.e. the forcing field that triggers the largest energy growth.More generally, resolvent analysis provides a transfer function of the system, characterising the response to a forcing for each frequency (Schmid 2007).
The linear forced dynamics can be studied by introducing a harmonic forcing term f in equation (2.9) as with  ∈ R which now represents the forcing frequency that one wants to study.Introducing the resolvent matrix R = (−I − J ) −1 , the relationship between any forcing and its linear response is then given by One can now look for the forcing field that triggers the most energetic response by computing the optimal gain  defined as Here, two norms must be introduced to measure the energy of the forcing and the response fields.The norm introduced by Chu (1965) is used for the response while the  2 -norm is used for the forcing.More details can be found in Bugeat et al. (2019).It can be shown that the optimal gain  corresponds to the largest singular value of the resolvent matrix modified by the norm matrices (Sipp & Marquet 2013).The associated optimal forcing and response are then the first left and right singular vectors, respectively.From equation (2.12),  can be seen as a transfer function as it measures the energy of the response with respect to that of an input forcing.In the following, the largest (or first) singular value will be called the optimal gain and will be referred to as  or  1 without any differences.The subsequent singular values  2 ,  3 , ... such that  1 >  2 >  3 > ... are called sub-optimal gains and are associated with the subsequent singular vectors which form an orthonormal basis.

Base flow computation
In order to get a 2D steady base flow, the non-linear Navier-Stokes equations (2.1) are solved using a finite volume CFD solver.Spatial discretisation of convective fluxes is performed using AUSM+ scheme (Liou 1996) associated with a fifth-order MUSCL extrapolation (Boin et al. 2006).Viscous fluxes at cell interfaces are obtained by a second-order centered finite difference scheme.The unsteady equations are marched in time until a steady state is reached.An implicit dual time stepping method with local time step is used (Jameson 1991).More details about this solver and its validation can be found in (Boin et al. 2006).A rectangular numerical domain is defined as in figure 1.The streamwise length   is set such that the distance between the impinging shock  imp and the outlet of the domain  out is 70 * 0 .The height   is set such that the reflecting shock leaves the domain at the outlet boundary.This value depends on the physical case considered.In any case, the baseflow eventually used in the stability and resolvent calculations is cropped to   = 50 * 0 .A Cartesian mesh is set with a geometrical progression from the wall.The number of points is 600 × 270 in the  and  directions, respectively.The independence of the results from the aforementioned numerical parameters is checked in appendix A. The following boundary conditions are used.At the inlet, oblique shock conditions are imposed: for  <  shock , a parallel flow is set using the Mach number at infinity while the corresponding downstream jump conditions are used for  >  shock .The latter condition is also set at top boundary of the domain.Note that  shock is calculated after choosing the angle of the shock  and the location of the impinging shock  imp .At the outlet, variables are extrapolated from the inside of the domain.At the bottom, an adiabatic flat plate is considered with no-slip conditions imposed on the velocity.

Linear solvers
The global stability and resolvent analysis are based on the linearised Navier-Stokes equations, which revolves around the computation of the Jacobian matrix J .In order to compute it, the discretised Navier-Stokes equation are linearised via a finite-difference

Validation
In order to validate the base flow calculation, the flow conditions are matched to those of the experimental investigation of Degrez et al. (1987).The freestream Mach number is  = 2.15 and the angle of the incident shock measured from the horizontal axis is  = 30.8• (corresponding to a flow deflection angle of 3.81 • ).The Reynolds number, based on the streamwise length from the leading edge at which the shock impinges the flat plate, is   = 0.96 × 10 5 .Comparisons with the experimental results reported by Degrez et al. (1987) are displayed in figure 2. This shows that our simulation is in good agreement with experimental data, both in terms of wall pressure measurements and location of the separation bubble.

Set of base flows studied
A set of base flows is computed in order to carry out stability and resolvent analysis for different separation lengths .Different Mach numbers are considered, ranging from 2.00 (incipient separation) to 2.35 (large separation length compared to the boundary layer length scale).First, the Reynolds number , based on the compressible boundary layer thickness  * 0 at which the shock impinges the plate, is kept constant at  = 1100 using the following protocol.For each Mach number, a boundary layer base flow is first computed without the impinging shock.At each -location,  * is then calculated by integrating the quantity 1 − ()/( ∞  ∞ ) in the normal direction.This allows us to detect the -location where  = 1100 and to use it as an input in a second base flow computation with a shock impinging at this point.Thus, the reference length scale  * 0 corresponds to the boundary layer thickness at the location of the impinging shock but prior to the shock interaction.In other words, it is not the effective thickness resulting from the shock interaction.The velocity field obtained at  = 2.20 and  = 1100, which will be used as the reference case, is shown in figure 3.By increasing the Mach number, the pressure gradient felt by the boundary layer is increased.The separation length / * 0 is readily computed by detecting the separation and reattachment points along the plate, which are the locations where / at the wall becomes negative and positive again, respectively.The separation length / * 0 is found very sensitive to  as its value increases over more than one order of magnitude for the Mach numbers considered in table 1.In order to assess the role of viscous effects in the low-frequency receptivity of the SWBLI, different Reynolds numbers are also considered, ranging from 600 to 2200.For one given Mach number, this corresponds to different -locations at which the shock impinges.Note that because  * ∼ √ , this location varies over a wide range of values, e.g. by a factor 10 at  = 2.20 for which the maximum value of  considered is 1900.As shown in table 1, increasing  leads to increasing the recirculation bubble as there is less streamwise momentum near the wall to resist the pressure gradient generated by the shock.Finally, note that the angle of the shock is kept constant at  = 30.8• for all the computations presented in this paper.

Results
The optimal gain is computed for different frequencies at  = 2.20 and  = 1100 (figure 4).At low frequency (  < 10 −1 ), the optimal gain monotonically decreases.The resolvent modes at   = 10 −4 are shown in figure 5.The optimal forcing is spread over a large portion of the domain, admitting maxima near the separation point and along the right part of the recirculation bubble.The optimal response is concentrated above the bubble, following its shape.Moreover, even though a drop of optimal gain is observed at low frequency, the optimal forcing and response vectors are found independent of the forcing frequency: the same velocity and pressure fields are observed at any low Strouhal numbers.Noticeable  differences start to appear for frequencies above   2 × 10 −2 .An interpretation will be proposed in section 4.1.3.At higher frequencies, a maximum of optimal gain is detected at   2 (figure 4).It is associated with a convective instability that starts developing in the recirculation region and continues to grow further downstream in the boundary layer (figure 6).This is reminiscent of the first mode instability of the supersonic boundary layer, which is the compressible counterpart of the Tollmien-Schlichting instability in incompressible boundary layers (Mack 1984).In order to efficiently trigger this instability, the optimal forcing is located upstream.Its tilted structure shows the concurrent action of the non-modal Orr mechanism, as usually observed in boundary layer flows (Ehrenstein & Gallaire 2008;Bugeat et al. 2019).

Global stability analysis
A global stability computation is performed for the set of base flows described in section 3.2.The eigenvalue spectrum at  = 2.20 and  = 1100 is presented in figure 7.Each  eigenvalue has a negative growth rate   meaning that the system is globally stable.This has also been verified for all the base flows considered in this study.The two least stable modes are found to be steady (S1 and S2,   = 0) while the third and fourth modes are unsteady (T1 and T2,   ≠ 0).Subsequent modes at higher frequency feature even larger damping rate.The velocity field of the mode S1 is mostly localised around the recirculation bubble, following its shape (figure 8-left).The mode T1 shares common characteristics but features a phase opposition between the upstream and downstream region of the bubble (figure 8-right).
The striking outcome of this global stability analysis is that the least stable global mode S1 is very similar to the optimal response found in the resolvent analysis in the previous section (figure 5-right), which was observed for any low frequencies.This means that, at low frequency, the optimal response results from the excitation of the least stable global mode.These results suggest that the receptivity at low frequency is a modal phenomenon, as opposed to the non-modal mechanisms usually observed in convective instabilities (Cossu & Chomaz 1997).Note that, because the system is stable, a continuous forcing is required to excite the flow around its base state.
It may appear somewhat unsettling that the receptivity, which pertains to the unsteady behaviour of the system, seems to be driven by a steady mode.This observation is analysed in the next section.

A low-pass filter model
To clarify how a steady mode could play a role in the low-frequency dynamics of the system, a model based on global stability analysis is proposed to describe the behaviour of the optimal gain.Two sources can cause the singular values of the resolvent operator to increase (Schmid 2007).The first one appears when the forcing frequency is close to the eigenvalues of the Jacobian matrix (as found by a stability analysis).The second source is related to the non-normality of the eigenvectors of the Jacobian matrix (i.e. the global modes).It is now  From the insight obtained from the previous stability analysis (section 4.1.2),we would like to test if the optimal gain proceeds from purely modal effects.Thus, non-modal phenomena are here discarded from the model.Furthermore, if one and only one global mode drives the optimal gain, then the singular value of the resolvent should vary as  ∼ 1/, where  is the distance in the complex plane between the forcing frequency   and the eigenvalue /2 = (  +   )/2 of this global mode.This distance is readily obtained as Because of the resemblance previously observed between the optimal response and the global mode S1, this mode is a natural candidate to test the purely modal receptivity model.Its eigenvalue is then plugged into equation (4.1).Since this mode is steady (  = 0), the low-frequency model of the optimal gain can eventually be recast as where  0 = lim   →0 (  ) is the value of the optimal gain when   goes zero.Here,  0 is not predicted a priori by the model, but is obtained by one computation of the optimal gain at zero frequency (or any small frequency where the optimal gain experiences a plateau).But since the absolute value of the optimal gain ultimately carries little significance (Sipp & Marquet 2013), this value is not essential to this model.The aim of this model is rather to detect the drop of optimal gain in order to analyse low-frequency receptivity of the system.It can be recognised that equation (4.2) is that of a first-order low-pass filter.The optimal gain being by definition analogous to a transfer function (see equation (2.12)), describing its behaviour as a filter is consistent.
In equation (4.2), the cut-off frequency is given by the damping rate  (1)

𝑖
of the global mode S1, which is the least stable global mode.Even though this mode is steady, a harmonic forcing of non-zero frequency is still able to excite it.When the forcing frequency   goes to zero, the gain does not depend on the frequency as the distance in the complex plane  tends to the constant value | (1)  /2|.However, increasing   until the order of  (1)  /2 and above affects the value of the optimal gain that then decreases as the forcing frequency is pushed away from the eigenvalue of the global mode it excites.Thus, the time scale appearing in this low-pass filter model is not set by the frequency of the global mode (which would be zero here), but by its damping rate.

Test of the model for different Mach and Reynolds numbers
The optimal gain model proposed in equation (4.2) is tested for different Mach and Reynolds numbers.Setting  = 1100, three Mach numbers are considered:  = 2.10,  = 2.20 and  = 2.30, with separation lengths ranging from 22 * 0 to 88 * 0 (table 1). Figure 9 shows that the low-pass filter model accurately detects the drop of optimal gain that appears between   = 10 −3 and 10 −2 .In other words, the cut-off frequency is correctly predicted by the damping rate of the global mode S1.This is observed for every Mach numbers.Disagreements appear above   2 × 10 −2 , where the low-pass filter underestimates the actual gain.At these frequencies, the structure of the optimal response starts to change.This means that the global mode S1 does not drive the dynamics of the flow alone any longer, but that other modes are getting involved.
Other Reynolds numbers, below and above the previous value  = 1100, are considered while keeping  = 2.2.Below  = 1100, the same agreement as previously described is observed (figure 10).As the Reynolds number increases ( = 1600), discrepancies start to appear for smaller frequencies than   = 10 −2 but the location of the drop of optimal gain is still correctly captured.At  = 1900, the model becomes irrelevant.The growing mismatch between the low-pass filter model and the optimal gain as  is increased can be caused by the action of other modes, previously damped by viscous effects.As a result, the model based on the excitation of the global mode S1 alone cannot predict the optimal gain any more, even though this mode does not directly depend on viscous effect (in the next section, it will be shown that the global mode S1 is independent of the viscous length scale  * 0 ).

Scaling
The low-pass filter model in equation (4.2) contains the damping rate  (1)

𝑖
of the mode S1.This parameter sets the unique time scale of the model.Its scaling relative to a reference length scale is now investigated.The non-dimensional damping rate  (1)  is presented in figure 11 as a function of the non-dimensional length of the recirculation bubble / * 0 .Note that each point is associated with a different pair of Reynolds and Mach numbers (see table 1).It is found that ∼  ∞ /.Because several Reynolds numbers (based on several  * 0 , as explained in section 3.2) have been tested, this , leading to values around   1 × 10 −2 .While of the same order of magnitude as the data of the literature, the cut-off frequency is slightly underestimated.This motivates to push the resolvent analysis a step further by looking at the sub-optimal modes.4.2.Resolvent analysis: sub-optimal forcing and response 4.2.1.Results So far, only the optimal gain  1 and its associated forcing and response modes have been considered.Subsequent singular values of the resolvent operator  2 and  3 , with  1 >  2 >  3 (i.e.sub-optimal gains), are now computed.Figure 12 compares these three gains as a function of the forcing frequency.At low frequency, a factor of about three is observed between  1 and  2 , and about two between  2 and  3 .A drop of sub-optimal gains occurs at higher frequencies that those noted for  1 .In particular, the decrease of  2 is observed at Strouhal numbers between 10 −2 and 10 −1 .These values lie in the range of the experimental data of the low-frequency unsteadiness more closely than those of  1 .Our attention will then be focused on the behaviour of  2 and its associated modes at low frequencies.The forcing mode follows the shape of the recirculation bubble but is mostly concentrated around the separation point as shown in figure 13-left.Unlike the optimal forcing in figure 5-left, the phase varies in the normal direction.This is also the case of the associated response (figure 13-right), whose phase varies in the streamwise direction.Finally, as previously observed for  1 , the sub-optimal modes associated with  2 do not depend on   in the low-frequency domain.This could suggest that the sub-optimal gain is also based on a modal excitation.In  the next section, this possibility is discussed by reconsidering the first-order low-pass filter model.

Sub-optimal gain: another low-pass filter
The response mode associated with  2 (figure 13-right) shares some features with the unsteady global mode T1 (figure 8-right) computed in section 4.1.2.Besides, the lowfrequency behaviour of  2 as a function of   is similar to that of  1 , the drop of gain being shifted to higher frequencies.It is then tempting to build a model of  2 based on the resonance of the global mode T1, derived from the distance in the complex plane between the real forcing frequency and its eigenvalue.However, because this mode is unsteady, a resonance model would lead to a peak of gain at the frequency of the mode T1 and a decrease of the gain as the frequency goes to zero.These two features are not observed in  2 (figure 12).Such a resonance model is then doomed to fail.Instead, equation (4.2), which is that of a first-order low-pass filter, satisfies the observation of  2 at low frequency.Therefore, we suggest to use this equation to model  2 by plugging the damping rate of the mode T1 instead of the mode S1 given its resemblance with the response mode associated with  2 .This equation becomes As previously mentioned, this approach is different from an approach based on the pure resonance of the global mode T1; instead, we only assume that the damping rate of this mode is the relevant time scale at play.The "purely modal dynamics" interpretation elaborated for  1 (see section 4.1.3)then no longer holds.Other modes, either by normal or non-normal effect, are understood to also contribute to the low-frequency dynamics.This is furthermore expected as the response mode of  2 is actually not strictly identical to the mode T1.This model is tested for the same set of base flow as previously investigated.At  = 1100, the location of the drop of the sub-optimal gain  2 , observed shortly after   = 10 −2 , is captured for the three Mach numbers considered as seen in figure 14.This means that the cut-off frequency of the low-pass filter, based on the damping rate of the global mode T1, is indeed relevant.The slope of the drop is also predicted by the first-order filter until   = 10 −1 , except for  = 2.10 for which the behaviour of  2 changes at lower frequencies.Now setting  = 2.20, the agreement of the cut-off frequency is also correctly captured for Reynolds numbers lower than or equal to  = 1100 (figure 15).At higher Reynolds number, the model over-estimates this frequency while still reasonably characterising the behaviour at  = 1600 .For  = 1900, the model completely fails (figure 15-bottom), as already noted for  1 .

Scaling
The damping rate  ( 1)  of the global mode T1 has been shown to capture reasonably well the low-frequency drop of  2 .Its scaling is considered by plotting its value against the non-dimensional separation length for the set of Reynolds and Mach numbers studied (figure 16).Similarly to section 4.1.5,the damping rate is found to scale with  ∞ / since  ( 1) This again means that the Strouhal number   associated with  cut-off frequency is constant for any Mach and Reynolds numbers.Its value is found around   4 × 10 −2 , in good agreement with the range of experimental data of the low-frequency unsteadiness.Before discussing the implication of this result, the bubble dynamics is studied in the next subsection.

Bubble dynamics
The first and second resolvent modes are now compared from the perspective of the bubble dynamics.In order to observe its motion, one response mode at one chosen frequency is linearly added to the base flow.The frequency considered here is   = 10 −3 but its value is actually not important since, as previously mentioned, the resolvent modes do not depend on   at low frequency.Snapshots at six different time steps of a periodic cycle, i.e. six different phases, are shown in figure 17.In each case, the distorted recirculation bubble is shown.It corresponds to the location in space where the total streamwise velocity is zero, calculated as ū +  û = 0.The scalar  is an arbitrary amplitude chosen as  = 0.3 to ease visualisation.The undistorted bubble of the base flow ( ū = 0) is also plotted as a reference.The optimal resolvent mode, associated with  1 , produces a breathing motion: the bubble periodically contracts and expands relatively to that of the basic state (figure 17-left).This dynamics is similar to the bubble motion found in the model of the low-frequency unsteadiness developed by Piponniau et al. (2009).The dynamics produced by the sub-optimal mode, shown in figure 17-right, is different.The bubble is now periodically displaced upstream and downstream from the reference state.No displacement is observed in the wall-normal direction.

Discussion
In their review paper, Clemens & Narayanaswamy (2014) ended their discussion on whether the SWBLI behaved as a forced dynamical system at low-frequency.Our resolvent analysis supports the idea that, indeed, the low-frequency behaviour results from a forced dynamics, as opposed to an intrinsic dynamics such as that of a vortex street behind cylinder.From the resemblance between their DMD and stable global stability modes, Nichols et al. (2017) suggested that the system could then be seen as a damped oscillator.Our study brings some new lights to this idea by showing that one stable global mode drives the resolvent analysis at low frequency.The fact that the optimal gain follows a first-order low-pass filter equation is furthermore consistent with the work of Touber & Sandham (2011) who were the first to model the low-frequency dynamics of the SWBLI as a filter of this kind.From this perspective, the dynamics requires a range of low frequencies that force the system and are filtered by the transfer function which, in the resolvent framework, is related to the singular values of the resolvent operator.Therefore, while carried out in the laminar regime, our analysis provides insights that the SWBLI in the turbulent regime could possess a forced dynamics at low frequency, in which fluctuations over a wide range of scales would force the flow (McKeon & Sharma 2010) and sustain its low-frequency behaviour.This would explain why low-frequency unsteadiness is observed in the turbulent regime rather that the laminar regime unless the later is externally forced as in Sansica et al. (2014).In the present study, because resolvent analysis assumes by nature the presence of a forcing field, the lowfrequency behaviour of the system can be qualitatively recovered even though the flow is laminar.
For both the optimal and sub-optimal gain, the damping rate of one stable global mode provides the time scale of the low-pass filter (the cut-off frequency).These frequencies scale as ( sep / ∞ ) −1 , which was expected as these modes are related to the bubble dynamic.This results in a constant cut-off Strouhal number   for any Mach and Reynolds numbers.This scaling is indeed observed both experimentally and numerically for different configuration of SWBLI, with frequencies reported from   3 × 10 −2 to   6 × 10 −2 (Dussauge et al. 2006).With values presently found around  (1)  1 × 10 −2 and  (2)  4 × 10 −2 for the optimal and sub-optimal gain respectively,  (2)   compares more favourably to the literature.However, the dynamics of the bubble of the optimal mode is, on one hand, related to a breathing motion that strongly resembles that at play in the low-frequency unsteadiness It is also strikingly similar to that observed experimentally in a turbulent incompressible recirculation bubble (Weiss et al. 2015), also occurring at Strouhal numbers of the order of 10 −2 .On the other hand, the sub-optimal mode that produces a very different type of motion, which is expected given the orthogonality of the resolvent modes.The resulting bubble motion is similar to that obtained by Nichols et al. (2017) in their global stability analysis; this is indeed the same unsteady stable mode.However, no conclusive element allows us to link it to the lowfrequency behaviour reported in the literature.This leads us to hypothesize that the optimal (rather than sub-optimal) gain and modes, associated with the bubble breathing, indeed describe the low-frequency dynamics of the SWBLI.The fact that  (1)  is underestimated in the present study may be due to the laminar nature of the SWBLI understudied while results in the literature usually deals with turbulent SWBLI.The shape of the bubbles being different between laminar and turbulent regimes this could account for the quantitative mismatch.Finally, while the range of values of separation length studied in the present work is large enough to characterise the scaling of the low-frequency time scale, the Mach number is restricted to a small range around M = 2.20.Similar analyses carried out at much different Mach numbers are then still needed to assess the universality of these findings.
The larger the gain separation  1 / 2 between the optimal and sub-optimal gain, the more likely the optimal response obtained through resolvent analysis is to model the system.More rigorously, this also depends on the projection of the effective forcing of the system (whether it is from external perturbations or intrinsic, turbulent non-linear interactions) onto the forcing resolvent modes (Beneddine et al. 2016).Here,  1 / 2 3 is found at low frequency.This value not being significantly greater than 1, if we assume that the optimal response drives the system (which we do based on the previous discussions on the bubble dynamics), then the actual forcing can be assumed to project reasonably well onto the optimal forcing mode.This gives us some insight on the actual forcing in the flow.The optimal forcing is indeed maximum around the separation point and along the separation line on the downstream part of the bubble (figure 5).Thus, some forcing is required in these regions in order to maintain the low-frequency bubble dynamics.The limit of the resolvent analysis is eventually reached when seeking for the origin of the forcing.Because of its linear nature, this framework does not allow us to discuss how the low-frequency structures can be produced by nonlinear interactions.Thus, both boundary layer structures convected to the interaction region (Ganapathisubramani et al. 2007) exciting the separation point and upstream travelling mechanisms (Sansica et al. 2016;Bonne et al. 2019) forcing the downstream part of the recirculation bubble could be valid candidates to sustain the low-frequency dynamics of the SWBLI.

Resolvent analysis
The low-frequency dynamics of three-dimensional perturbations, characterised by their nonzero spanwise wave number , is studied in this section.The analysis is here limited to the base flow calculated at  = 2.20 and  = 1100.The influence of  on the optimal gain at   = 10 −4 is shown in figure 18.A first peak is observed for low wave numbers around  = 0.25.The corresponding wavelength is   = 25 * 0 , which is of the order of the recirculation length (approximately half of it, see table 1).The associated optimal forcing and response can be visualised by plotting the isosurface of a selected amplitude in the 3D physical space (figure 19).The optimal response is located in the recirculation region and features a phase shift in the streamwise direction.A strong resemblance can be noted with the 3D global mode found by Robinet (2007).The possibility of a modal resonance to model the optimal gain, as developed previously for 2D perturbation, will be explored in the next section.The optimal forcing field is made of elongated structures that lay upstream from the recirculation bubble until the separation point.The energy density d () of the response can be defined at each streamwise location as d () = ∫ 0.5 ( * +  * +  * ) d.An analogous quantity can be defined for the forcing field.The profiles of these energy densities confirm that the response is localised in the separation region as it experiences a rapid decay downstream from it (figure 20).The forcing energy is shown to spread over a larger portion of the domain and reaches a peak in the vicinity of the separation point.
In addition of the local maximum at  = 0.25, the optimal gain features a peak of larger  magnitude around  = 2.The corresponding wavelength is   = 3 * 0 and is now of the order of the boundary layer thickness.Elongated streamwise vortices located upstream of the recirculation bubble are found in the optimal forcing field (figure 21).A further inspection of the wall normal organisation of its components (available in appendix B, figure 28) shows that they correspond to streamwise vortices.Downstream of the recirculation bubble, streaks of streamwise velocities are observed in the boundary layer, for which profiles reveal that  >>  and  >>  (figure 28).These features indicate the action of the lift-up mechanism that is ubiquitous in shear flows at low frequency and non-zero spanwise wave numbers (Landahl 1980).Besides, these results are very similar to the 3D optimal response in a supersonic boundary layer without impinging shock (Bugeat et al. 2019).The energy profiles reveal the convective nature of this instability as perturbations keep growing until the downstream end of the domain (figure 20).Most of the energy is contained downstream from the recirculation area whereas most of the forcing energy is localised upstream from it.These results hence show that the growth of streaks does not interact significantly with the recirculation region.

Link with global stability
Eigenvalue spectra of the global stability problem are given in figure 22 for the wave numbers  = 0.25 and  = 2 that correspond to maxima of the optimal gain.In both cases, the flow is found to be globally stable.Similarly to the analysis developed for 2D perturbations in section 4, the relation between the optimal gain behaviour at low frequency and the global stability features is explored for these two spanwise wave numbers.

Low 𝛽: excitation of a steady global mode
At  = 0.25, the least stable global mode is steady (mode S1) as shown in figure 22-left.Its structure is shown in figure 23.It strongly resembles that of the optimal response previously observed.For 2D perturbations, these observations led us to formulate the hypothesis of a purely modal dynamics at low frequency, driven by a single mode.The first-order low-pass filter model (equation (4.2)) that results from this assumption and that is based on the damping rate of the 3D global mode S1 is tested in figure 24.Once again, the model correctly predicts the cut-off frequency of the optimal gain.The slope of the gain drop is also reasonably captured until   10 −2 where additional modes participate to the forced dynamics.The low-frequency dynamics of 3D perturbations at low wave numbers thus behave the same way as that of 2D perturbations, which was expected as it is the limit when  → 0. Finally, another peak of optimal gain can be observed at higher frequency in figure 24.This maximum is associated with the convective instabilities that are known to exhibit larger growth rates for  ≠ 0 in supersonic flows (Mack 1984), resulting in larger values of the optimal gain (Bugeat et al. 2019).

Larger 𝛽: pseudo-resonance
The lift-up mechanism leading to the formation of streaks is a linear algebraic instability, as theorised by Ellingsen & Palm (1975).As such, linear stability theories based on a modal approach cannot capture it.Conversely, the resolvent analysis is particularly suited for their calculation as non-modal effects are taken into account (Schmid 2007).We now suggest to show the failure of the first-order low-pass filter model (based the excitation of a single mode) as a way to point out the non-modal character, or pseudo-resonance, of the optimal gain at wave numbers of the order of  = 2.The least damped mode of the spectrum is the unsteady modes T1, followed by the steady mode S1 (figure 22-right).Their structure is shown in figure 25.They share similarities with the optimal response already presented in figure 21-right as their growth takes place in the downstream region of the domain.Nevertheless, noticeable differences remain, which suggests that the optimal response does not proceed from the excitation of a unique global mode.Two low-pass filters model based on the modes S1 and T1 are tested in figure 26.Both models underestimate the cut-off frequency and underpredict the optimal gain, confirming the non-modal nature of the optimal response at these wave numbers.

Discussion
Depending on their spanwise wave number, the 3D optimal response at low frequency can result from two different mechanisms, underpinned by a modal or non-modal nature.At low wave numbers, whose length scales are of the order of  sep , optimal perturbations interact with the recirculation bubble.Their energy is localised in this region similarly to what was observed for 2D perturbations.Additionally, the low-pass filter model based on the excitation of a unique steady global mode reasonably describes the low-frequency behaviour of the optimal gain .Low wave numbers hence behave in continuity to the 2D results ( = 0) presented in section 4. Noticeably larger values of  are however reached for 3D perturbations.This raises the question of whether the dynamics of the system at low frequency and low wave number is dominated by 2D or 3D perturbations.On the one hand, larger values of the gain indeed indicates that more energetic 3D perturbations are expected in the flow, assuming the same forcing energy and the same projection of this forcing onto the resolvent modes for each wave number.These two assumptions are, at this point, impossible to verify.For example, if the forcing generated by non-linear interactions in a turbulent flow is preferentially 2D, then one may still expect the 2D response to be dominant even though it is associated with a lower gain value.Thus, no definitive conclusion can be drawn.On the other hand, the cut-off frequency based on the 3D global is   = 6 × 10 −3 .This underestimates the usual value of   = 3 × 10 −2 even more than the 2D perturbations, for which a value of   = 1 × 10 −2 was found.This is in favour of the prevalence of 2D perturbations at low frequency even though a perfect quantitative agreement between our laminar analysis and the experimental data on turbulent flow is not expected.This is also backed up by previous work which restricted their analysis to 2D perturbations (Touber & Sandham 2011;Sansica et al. 2014) The possibility that the low-frequency behaviour of the SWBLI may be a 3D phenomenon still remains a promising lead to explore further given some arguments found in our resolvent analysis.It could be interesting to investigate whether the optimal spanwise wave number would then scale with the recirculation bubble.A space-Fourier transform along the spanwise direction of the time-Fourier modes from experimental or numerical data would be required to explore this phenomenon.This, however, may not be manageable experimentally if the spanwise width of the facility is not wide enough compare to recirculation length scale.Recent studies have besides pointed out the influence of the side-walls in SWBLI experiments (Xiang & Babinsky 2019;Rabey et al. 2019).From a computational point of view, the same challenge of ensuring sufficiently wide domain is expected, additionally coupled with the costly computational time necessary to statistically converge the low-frequency time-Fourier transform.Such a study, if achievable, could help to characterise the 3D organisation of the forcing field and eventually understand its role in the low-frequency dynamics.
While low wave numbers exhibit a modal and localised (around the bubble) dynamics, optimal responses at larger wave numbers have been shown to ignore the recirculation region as they are driven by a convective instability growing downstream in the boundary layer, as a result of non-modal effects.These findings are consistent with those of Dwivedi et al.
(2020) who recently carried out a transient growth analysis in a SWBLI at  = 5.92.Such analysis also embeds non-modal effect and can in fact be seen as the temporal counterpart of the resolvent analysis (Sipp et al. 2010).The authors detected streaky structures evolving in time and convected in space, similar to those obtained in the present study.Our resolvent analysis showed that streaks are associated with larger optimal gain values than the lowwave number perturbations.This is not surprising as it results from a convective instability which efficiently extract energy from the base flow.The role played by the streaks in the low-frequency dynamics (  3 × 10 −2 ) is not obvious.Contrary to low wave number perturbations, streaks are not directly involved in this dynamics since they do not interact with the recirculation region and are expected to scale with  * 0 rather than .At first, it then seems that they could only be of interest regarding the transition to turbulence of the system.However, the large energy growth they trigger could quickly lead to non-linear interactions that may play a role in the forcing of the bubble through upstream feedback (Bonne et al. 2019).Studying the development of streaks in a turbulent SWBLI, seen as coherent structures, would then be of interest.Following the work of Sartor et al. (2015), this could be achieved by performing this 3D resolvent analysis on a turbulent mean flow.

Conclusion
A resolvent analysis was carried out on a laminar SWBLI in order to study the low-frequency behaviour of the system for both 2D and 3D perturbations.The optimal gain, that can be seen as a transfer function of the system, was computed over several decades of Strouhal numbers.A model of optimal gain based on the excitation of a single global mode was developed, leading to the equation of a first-order low pass-filter in agreement with the work of Touber & Sandham (2011).Our analysis was shown to be valid for 2D perturbations over different base flows calculated from different pairs of Reynolds and Mach numbers, but was not relevant at the largest Reynolds numbers tested.The time scale of the model, i.e. the cut-off frequency of the filter, was found to be determined by the damping rate of the least damped global mode.This damping rate scales as  ∞ /.Thus, the low-frequency dynamics is associated with a constant Strouhal number based on the recirculation length, as observed experimentally and numerically.The value obtained here is   1 × 10 −2 , which slightly underestimates the classical values of the literature.This mismatch could be attributed to the laminar nature of the flow studied in this paper.Ultimately, this work carried out in the laminar regime suggests that the low-frequency dynamics of the SWBLI in the turbulent regime could be a forced dynamics, sustained by the background turbulent fluctuations.The resulting bubble dynamics, which would be sustained by a background forcing produced by non-linear interactions in a turbulent flow, is then found to exhibit a breathing motion.
Assuming the spanwise direction as homogeneous, a 3D resolvent analysis showed that two regimes of spanwise wave numbers  could be identified.At low , the same modal mechanism as that of 2D perturbations, leading to a first-order low-pass filter, was observed.Larger values of optimal gain were obtained for wave numbers of the order of the recirculation length.The possibility of a prevailing 3D dynamics was then discussed.At larger wave numbers, streaks generated by the non-modal lift-up effect were detected.Their growth does not interact with the recirculation region, suggesting that streaks are not directly involved in the low-frequency dynamics of the SWBLI.Finally, from the perspective of transition to turbulence, 3D resolvent analysis remains an efficient tool to characterise the energy growth of streaks and provides their optimal wave number and forcing location.
Future work investigating the forcing terms (both 2D and 3D) produced by non-linear interactions in a turbulent SWBLI, either through DNS or experimental observations, would be a decisive follow-up to this work.While the low-frequency behaviour of the SWBLI was suggested to be understood as a forced dynamics, the physical origin of the forcing is missing.Ultimately, this could also provide some insights to assess the prevalence of 3D effects in the dynamics of the system.line of table 2. Independency with respect to the number of grid points   and   is verified in the cases A and B. The influence of the height of the domain   / * 0 and the length of the domain downstream from the impinging shock    / * 0 are assessed in the cases C and D. Overall, the unsteady mode is the most easily converged mode: relative differences of eigenvalues between the reference case and the other ones are equal to or below than 0.1% for the range of tested parameters.As for the steady mode, relative difference of eigenvalues lower than 1% are found, which is deemed acceptable.

Figure 1 :
Figure 1: Sketch of the numerical domain and typical structure of the OSWBLI

Figure 2 :
Figure 2: Validation against the experimental data from Degrez et al. (1987).Left: wall pressure normalised by the minimum pressure  0 upstream from the shock.Right: Skin friction coefficient   .The blue squares correspond to the experimental separation and reattachment points.

Figure 7 :
Figure 7: Global stability spectrum at  = 2.20,  = 1100.Steady and unsteady modes are referred to with the letter S and T, respectively.

Figure 9 :Figure 10 :
Figure 9: Comparison between the optimal gain (red squares) and the low-pass filter model from equation (4.2) (blue line) for different Mach numbers ( = 1100).

Figure 11 :
Figure 11: Damping rate of the global mode S1 for different Mach and Reynolds numbers, as a function of the separation length.Here,  * 0 is the boundary layer thickness at  = 1100 rather than that at each Reynolds number.

Figure 14 :
Figure 14: Comparison between the sub-optimal gain  2 (red squares) and the low-pass filter model in equation (4.3) (blue line) for different Mach numbers ( = 1100).

Figure 15 :
Figure 15: Comparison between the sub-optimal gain  2 (red squares) and the low-pass filter model in equation (4.3) (blue line) for different Reynolds numbers ( = 2.20).

Figure 16 :
Figure 16: Damping rate of the global mode T1 for different Mach and Reynolds numbers, as a function of the separation length.Here,  * 0 is the boundary layer thickness at  = 1100 rather than that at each Reynolds number.

Figure 17 :
Figure17: Bubble dynamics at  = 2.20 and  = 1100 resulting from the addition of the optimal (left) and sub-optimal (right) low-frequency resolvent modes to the base flow.The blue dashed lines is the distorted recirculation bubble corresponding to the location where ū +  û = 0 (using  = 0.3).Its evolution at 6 different phases of a periodic cycle is shown from top to bottom.The red line is the recirculation bubble of the base flow ( ū = 0).

Figure 19 :Figure 20 :
Figure 19: Optimal resolvent mode at  = 0.25 and   = 10 −4 , real part of .Iso-surfaces at 40% and −40% the maximum absolute value are plotted in red and blue.The dark grey surface is the recirculation bubble of the base flow (homogeneous in ).

Figure 21 :
Figure 21: Optimal resolvent mode at  = 2 and   = 10 −4 , real part of  of the forcing (left) and real part of  of the response (right).Isosurfaces at 20% and −20% of the maximum absolute value are plotted.

Figure 23 :Figure 24 :
Figure 23: Global mode S1 at  = 0.25,  = 2.20 and  = 1100.Isosurfaces of 40% and −40% of the maximum absolute value of  are plotted in red and blue.

Figure 25 :Figure 26 :
Figure 25: Global modes S1 (left) and T1 (right) at  = 2 and   = 10 −4 , real part of .Isosurfaces of 20% and −20% of the maximum absolute value amplitude are plotted in red and blue.

Table 1 :
Set of base flows studied with the length of the separation region  obtained in each case.

Table 2 :
Name        / * 0   / * 0  (1) × 10 4 |Δ (1) ( 1) × 10 3 |Δ ( 1) Mesh convergence of the most unstable mode at  = 2.30,  = 1100.The reference set of numerical parameters, used in all the calculation of the paper, is the first one of the table. Δ ref corresponds ts to the difference between the eigenvalue of the reference case and that obtained in each other case.