# Modal analysis of the triadic interactions in the dynamics of a transitional shock wave boundary layer interaction

Ismaël Ben Hassan Saïdi, Stéphane Wang, Guillaume Fournier, Christian

Tenaud, Jean-Christophe Robinet

► To cite this version:

Ismaël Ben Hassan Saïdi, Stéphane Wang, Guillaume Fournier, Christian Tenaud, Jean-Christophe Robinet. Modal analysis of the triadic interactions in the dynamics of a transitional shock wave boundary layer interaction. Journal of Fluid Mechanics, 2025, 1009, pp.A43. 10.1017/jfm.2025.236. hal-05063074

HAL Id: hal-05063074

https://hal.science/hal-05063074v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Modal analysis of the triadic interactions in the dynamics of a transitional shock wave boundary layer interaction

Ismaël Ben Hassan Saïdi¹ ID, Stéphane Wang², Guillaume Fournier³ ID,
Christian Tenaud²,⁴ ID and Jean-Christophe Robinet¹ ID

¹École Nationale Supérieure des Arts et Métiers, Laboratoire DynFluid, Paris 75013, France

²Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire EM2C, Gif-sur-Yvette 91190, France

³Université Paris-Saclay, Université d'Évry, LMEE, Évry 91020, France

⁴CNRS, Fédération de Mathématiques de Centrale Supélec, Université Paris-Saclay, France

**Corresponding author:** Jean-Christophe Robinet, Jean-christophe.ROBINET@ensam.eu

(Received 15 March 2024; revised 22 December 2024; accepted 17 February 2025)

This work is a numerical study of a transitional shock wave boundary layer interaction (SWBLI). The main goal is to improve our understanding of the well-known low-frequency SWBLI unsteadiness and especially the suspected role of triadic interactions in the underlying physical mechanism. To this end, a direct numerical simulation is performed using a high-order finite-volume scheme equipped with a suitable shock capturing procedure. The resulting database is then extensively post-processed in order to extract the main dynamical features of the interaction zone dynamics (involved characteristic frequencies, characteristics of the vortical structures, etc.). The dynamical organisation and space-time evolution of the flow at dominant frequencies are then further characterised by mean of an spectral proper orthogonal decomposition analysis. In order to study the role of triadic interactions occurring in the interaction region, a bispectral mode decomposition analysis is applied to the database. It allows us to extract the significant triadic interactions, their location and the resulting physical spatial modes. Strong triadic interactions are detected in the downstream part of the separation bubble whose role on the low-frequency unsteadiness is characterised. All the results of the various analyses are then discussed and integrated to formulate a possible mechanism fuelling low-frequency SWBLI unsteadiness.

**Key words:** high-speed flow, supersonic flow, shock waves
# 1. Introduction

Situations in which a shock wave interacts with a boundary layer are numerous in the aeronautical and space industries. These interactions can exist on external surfaces (transonic profiles, junctions of surfaces, etc.) or in internal aerodynamics devices (supersonic air intakes, cascade of turbine blades, nozzles, etc.). Under certain conditions (high Mach number, large shock wave angle, etc.), these interactions can generate a transient separation bubble causing increased drag force, heat fluxes and pressure fluctuations. It is also known since decades that the separation bubble and reflected shock wave are subject to low-frequency motion, known as ‘shock wave boundary layer interaction (SWBLI) unsteadiness’. For turbulent interactions, with a turbulent incoming boundary layer, the characteristic frequency of these oscillations is two orders of magnitude lower than the characteristic frequencies of the incoming boundary layer (Dussauge, Dupont & Debiève 2006). The SWBLI unsteadiness can be detrimental to engineering system performances and can expose structures to oscillating loads, potentially damaging the solid structure’s integrity (Dolling 2001; Délery & Dussauge 2009; Babinsky & Harvey 2011; Clemens & Narayanaswamy 2014; Gaitonde 2015). During the last decades, attention has been focused on numerical and experimental studies of unswept SWBLIs (where the shock impingement line is orthogonal to the incident boundary layer), employing various analysis techniques like Fourier analysis of signals (e.g. pressure or velocity probe signals), modal decomposition (proper orthogonal decomposition (POD), Sirovich 1987; Shinde et al. 2019), spectral POD (SPOD) (Towne, Schmidt & Colonius 2018), dynamic mode decomposition (DMD) (Schmid 2010; Priebe et al. 2016)) and stability analysis (Theofilis 2003, 2011; Robinet 2007; Sartor et al. 2015; Guiho, Alizard & Robinet 2016; Song & Hao 2023).

The research community identifies two main mechanisms for SWBLI unsteadiness: (A) upstream, caused by the advection of large-scale structures of the incoming boundary layer (Ganapathisubramani et al. 2007, 2009) and (B) downstream, in which the dynamics of the separation bubble generates disturbances that drive the oscillations of the reflected shock. As stated in Clemens & Narayanaswamy (2014), following Souverein et al. (2010), a consensus view emerged in the community that strong interactions, which exhibit large separation bubbles, are primarily driven by a downstream instability (mechanisms of type B), whereas weakly separated interactions (also called incipient interactions) can be strongly influenced by fluctuations in the upstream boundary layer (mechanisms of type A).

Regarding strong interactions, a large number of studies have been undertaken in the last decades to document the low-frequency dynamics of the separation bubble and to identify the exact physical mechanism of type B underlying the SWBLI unsteadiness. It is now well established that the dynamics of the flow around this mean configuration is characterised by several unsteady phenomena whose characteristic scales spread over a large broadband spectrum range. As illustrated in figure 1, the main dynamical features of the flow can be classified according to three frequency ranges spreading over two decades: high, medium and low frequencies. For turbulent interactions, the incoming boundary layer is turbulent with the most energetic fluctuations at high frequencies characterised by a Strouhal number $St_δ = fδ/U_e \sim 1$ (where f, δ and U_e are respectively the main frequency of the fluctuations, the boundary layer thickness before the interaction and the free-stream velocity). Various studies carried out to characterise the dynamics of subsonic turbulent separated and reattached flows have highlighted a now well-documented medium-frequency dynamics of the separation bubble (Weiss, Mohammed-Taifour & Schwaab 2015; Wu, Meneveau & Mittal 2020). In this range of frequencies the shear layer, bounding the upper part of the bubble, is subjected to a global flapping motion of the separated
<figure><img src="image_3.png" /><figcaption>Figure 1. Illustration of the documented dynamical features of SWBLIs. The figures are from Ben Hassan Saidi (2019), Kirk & Yarusevych (2017) and Buffin-Bélanger et al. (2013).</figcaption></figure>

shear layer at a Strouhal number $St_L = fL/U_e \sim 0.1–0.2$ (where L is the mean separation length) and a higher medium-frequency instability linked to the shedding of large-scale vortical structures at $St_L \sim 0.3–0.8$ near the reattachment. The flapping mode consists of successive enlargement and shrinkage of the recirculation bubble, the shrinkage being associated with the shedding of a large vortex downstream of the recirculation bubble. These two medium-frequency unsteadinesses correspond to the modes that have long been documented for subsonic fixed-separation bubbles (Cherry, Hillier & Latour 1984; Kiya & Sasaki 1985). For these fixed-separation bubbles, the flapping frequency was reported as the lowest frequency present in the flow. Nevertheless, in turbulent SWBLIs, a low-frequency flapping mode was also observed at a Strouhal number of $St_L = fL/U_e = 0.02–0.05$. In the following, we refer to this low-frequency mode as the ‘breathing’ of the separated zone. This mode is accompanied by the oscillations of the reflected shock wave, in the same frequency range, which constitute the so-called SWBLI unsteadiness. Thus, the low-frequency separation bubble oscillations documented for SWBLIs were one order of magnitude lower than those documented for fixed-separation subsonic separation bubbles. This led to a series of research projects aimed at explaining the low-frequency oscillations of SWBLIs as being a specificity of these flows, namely dependent on the compressible nature of the flow or the presence of the shock wave (Pirozzoli & Grasso 2006; Piponniau et al. 2009; Touber & Sandham 2011; Aubard, Gloerfelt & Robinet 2013; Priebe et al. 2016; Adler & Gaitonde 2018). This view of the problem has been challenged in the recent work of Weiss et al. (2015) studying the dynamics of pressure induced subsonic turbulent separation bubbles on flat test surfaces in which a breathing mode at a frequency of the same order of magnitude than the breathing of SWBLI induced separation bubbles have been documented. This result tends to demonstrate that the low-frequency breathing mode is an intrinsic characteristic of the dynamics of separation bubbles, where the position of the separation point is not imposed by the geometry, independently of the compressible nature of the flow.

Another furrow dug for decades to identify the physical mechanism behind the low-frequency SWBLI unsteadiness is the simulation and analysis of transitional SWBLIs with forced or non-forced laminar incoming boundary layers in order to remove the presumed influence of large-scale turbulent structures in the incident boundary layer. This approach allows us to carefully analyse the dynamics of the separation bubble that is suspected to be the source of the low-frequency SWBLI unsteadiness (Robinet 2007;
I. Ben Hassan Saïdi, S. Wang, G. Fournier, C. Tenaud and J.-C. Robinet

Sansica, Sandham & Hu 2014; Guiho et al. 2016; Diop, Piponniau & Dupont 2019; Bugeat et al. 2022; Mauriello, Larcheveque & Dupont 2022; Mahalingesh, Piponniau & Dupont 2023; Song & Hao 2023). Part of these works involved characterising the linear dynamics of the flow by means of stability analyses. The global (self-sustained dynamics) linear stability analysis performed by Robinet (2007) on a transitional incident oblique shock wave SWBLI on a flat plate at $M = 2.15$ demonstrated that the fixed point of this flow is stable for two-dimensional (2-D) perturbations. This result has then been confirmed in Guiho et al. (2016). For three-dimensional (3-D) perturbations, the flow exhibited a global instability when the shock angle exceeded a certain critical value. This result has then been confirmed by several studies (Song & Hao 2023). This mode is however stationary and, therefore, not directly related to the self-sustained low-frequency dynamics. As the low-frequency unsteadiness could not be directly linked to any global mode around the fixed point, the receptivity of the flow has been studied by means of resolvent analysis. In particular, the results obtained by Bugeat et al. (2022) show that the global resolvent analysis no frequency selectivity at low frequency and, therefore, behaves as a low-pass filter with respect to external disturbances and no clear evidence of quasi-resonance has been found. In light of the results of both global stability and resolvent analyses cited above, potential purely linear mechanisms explaining the low-frequency unsteadiness rely on the existence of low-frequency external perturbations of the flow. This brings us to the limits of a purely linear analysis of the low-frequency dynamics of SWBLIs. Indeed, the low-frequency SWBLI unsteadiness has been documented in flow configurations with a laminar incoming boundary layer, without any imposed low-frequency forcing (Sansica, Sandham & Hu 2016; Ben Hassan Saidi 2019; Mauriello et al. 2022; Mahalingesh et al. 2023). The search for nonlinear mechanisms playing a role in the generation of low frequencies is therefore a promising avenue. This approach is encouraged in particular by the work of Sansica et al. (2016), who studied a case of transitional interaction forced upstream by unstable high-frequency modes. In this SWBLI, analysis of the Fourier spectra of the parietal pressure showed that the low frequencies appear to be generated in the transition zone towards turbulence (suggesting the nonlinear nature of the mechanism), close to downstream of the reattachment and be travelling upstream in the separation bubble. In a more recent study of another upstream forced (with high frequencies) transitional SWBLI (Mauriello et al. 2022), nonlinear spectral analyses (bispectrum and bicoherence) have been performed that pointed out quadratic coupling between frequencies as a possible mechanism involved in the low-frequency dynamics of the separation bubble. This work highlights the appearance of energy at low frequency near the reattachment point as a result of quadratic coupling with linear unstable modes developing within the shear layer. This energy at low frequency is then shown to be convected upstream through a feedback mechanism. These recent works point out the possible importance of triadic interactions in the mechanism driving the breathing of the separation bubble. However, these previous studies were based on local analyses of the signals, which did not allow us to determine precisely the spatial location of the nonlinear interactions in the interaction zone. In addition, the methods used did not allow the physical structure of the modes resulting from these nonlinear interactions to be determined.

The objective of the present work is to go further in the analysis of the triadic interactions in the separation bubble of SWBLIs and their involvement in the mechanisms underlying the low-frequency unsteadiness of SWBLIs. To this end, direct numerical simulation (DNS) of a non-forced transitional oblique shock wave SWBLI is performed. We study nonlinear coupling between modes at different frequencies using bispectral mode decomposition (BSMD) (Schmidt 2020), which is a modal decomposition method specially designed to study triadic interactions in flows. The use of this method is of
major interest compared with previous studies because it is designed to analyse triadic
interactions in 3-D signals. As a result, this algorithm allows us to highlight quadratic
coupling between modes at different frequencies. This approach allows us to associate a
spatial mode to each frequency resulting from a triadic interaction. Moreover, the location
of the interactions between two frequencies in the flow field can be represented in a so-
called interaction map. These two outputs of this method represent a major advantage
with respect to local methods. The main objective of this analysis is to document the
suspected nonlinear link between the dynamical activity at medium and low frequencies,
given that the presence of medium-frequency modes of the shear layer are already strongly
documented. The underlying objective is to demonstrate that the nonlinear evolution of the
medium-frequency shear layer dynamics create the frequencies characterising the SWBLI
unsteadiness. Spectral POD (Towne et al. 2018) has also been used as a support for the
physical interpretation of the involved modes. Indeed, for each frequency, SPOD performs
a decomposition of the flow in spatio-temporal modes ordered by their energy content,
highlighting the structure of the flow for a given frequency.

In our study we have deliberately chosen a transitional configuration in which the incident boundary layer is laminar and not forced. In this configuration, the transition to turbulence is naturally obtained by the imposed shock pressure gradient. In this respect, the present case study represents a limit case in which the mechanisms at the origin of the low-frequency dynamics appear to be present at low intensity. Moreover, the flow dynamics is not disturbed by the natural frequencies of the incident turbulence. The characteristic frequencies and mechanisms associated with the low and medium frequencies can therefore be studied independently of the mechanisms and frequencies specific to the turbulence in the upstream boundary layer. We therefore expect, so to speak, to find in the flow the main features and mechanisms of the dynamics we are interested in, but no more. We believe that this characteristic makes the present flow configuration a good candidate for identifying the mechanisms underlying the low frequency of the separation bubble by modal analysis; in particular, modal analysis associated with the search for triadic interactions in the interaction zone. Indeed, a more complex dynamics induced by the presence of the turbulence or the forcing of the incident dynamics would imply a spurious triadic interaction, not relevant for present study.

In the following, we first present the numerical simulation and database. The studied flow configuration is described as well as the governing equations and numerical strategy used for the simulation. The salient features of the mean flow and its dynamics is documented and the database used for both SPOD and BSMD analysis are precisely described. We then present the SPOD analysis. We particularly describe the dominant modes of the flow at medium and low frequencies. The BSMD analysis follows, in which the triadic interactions in the interaction region are documented. The physical relevance of the results regarding the mechanism at the origin of the SWBLI unsteadiness is the discussed. The paper ends with concluding remarks and perspectives.

## 2. Numerical simulation and database

2.1. *Flow configuration*

The flow consists of the interaction between an incident shock wave and a laminar
boundary layer developing on a flat plate. A 2-D sketch of the flow (2-D slice) is shown
in figure 2. The presence of the impinging shock wave (also called the incident shock
wave, noted as (S<sub>I</sub>)) imposes a sharp adverse pressure gradient to the boundary layer
of thickness δ (99 % thickness) evaluated just upstream of the interaction region. As
I. Ben Hassan Saïdi, S. Wang, G. Fournier, C. Tenaud and J.-C. Robinet

<figure><img src="image_4.png" /><figcaption>Figure 2. Sketch of a 2-D slice of the flow, where ($S_I$) is the incident shock wave of angle $\alpha$, $\delta$ is the boundary layer thickness just upstream of the interaction zone, ($S_R$) is the reflected shock wave, $\bar{x}_S$ is the mean separation point, $\bar{x}_R$ is the mean reattachment point, $x_{sh}$ is the location of the impingement shock wave (if the flow were non-viscous) and $L_{int} = x_{sh} - \bar{x}_S$ is the interaction length.</figcaption></figure>

accurately described by Délery & Dussauge (2009), in the configuration presently studied, the adverse pressure gradient is strong enough to initiate a separation of the boundary layer that reattaches further downstream forming a closed separation bubble. The mean separation and reattachment points are respectively noted as $\bar{x}_S$ and $\bar{x}_R$. The flow being subsonic in the wall vicinity, the pressure rise due to the incident shock wave is sensed upstream to the location, noted as $x_{sh}$, where the incident shock wave would impact the wall (if the flow were non-viscous) because of slow acoustic waves, explaining the location of the separation point, $\bar{x}_S$, upstream of the location of the incident shock wave impingement $x_{sh}$. At the separation point, the deviation of the supersonic flow due to the separation leads to the formation of the so-called reflected shock wave, noted as ($S_R$). In fact, the incident shock is reflected at the apex of the separation bubble as expansion waves. At the reattachment point, the deviation of the supersonic flow due to the presence of the wall leads to compression waves that can also coalesce to form the so-called reattachment shock wave. For given Reynolds and Mach numbers, the extent of the separation bubble is driven by the intensity of the incident shock wave, let us say the pressure ratio from each side of the shock, prescribed by the shock angle $\alpha$. In the following, we characterise this extent by the interaction length $L_{int} = x_{sh} - \bar{x}_S$.

The free-stream conditions are similar to the experimental and numerical test case documented in Degrez, Boccadoro & Wendt (1987), with a free-stream Mach number $M = 2.15$ and a free-stream Reynolds number $Re = 1600$, based on the boundary layer thickness just upstream of the interaction zone $\delta = 1.6 \times 10^{-2}$ m. Compared with the configuration of Degrez et al. (1987), where the whole flow is steady (fully laminar interaction), the incident shock wave angle $\alpha$ is increased up to 33.8°. The interaction is thus strengthened, the separation bubble is enlarged and the dynamics of the separation bubble becomes unsteady with transition to turbulence triggered in the bubble (Ben Hassan Saïdi 2019). Table 1 summarises the physical parameters of the simulation.

## 2.2. Governing equations

To perform the DNS, we consider the dimensionless compressible Navier-Stokes equations expressed in a Cartesian coordinate system. The length and time scales for non-dimensionalisation are respectively $x_{sh}$ and $U_\infty/x_{sh}$, with $U_\infty$ being the free-stream
<table><caption>Table 1. Flow parameters of the SWBLI.</caption><tr><td>Mach number</td><td><i>M</i> = 2.15</td></tr><tr><td>Reynolds number</td><td><i>Re</i> = 10<sup>5</sup></td></tr><tr><td>Prandtl number</td><td><i>Pr</i> = 0.71</td></tr><tr><td>Shock wave angle</td><td><i>α</i> = 33.8</td></tr><tr><td>Ratio of specific heats</td><td><i>γ</i> = 1.4</td></tr></table>

velocity. All the primitive variables are non-dimensionalised by the corresponding free-stream quantities except the thermodynamic pressure $P^*$ and the total energy per unit volume $\rho^*E^*$ that are non-dimensionalised by the free-stream dynamic pressure $\rho_\infty U_\infty^2$. The dimensionless Navier-Stokes system is therefore written as

$$ \frac{\partial U}{\partial t} + \nabla \cdot F(U) - \nabla \cdot F_v(U, \nabla U) = 0, \quad (2.1) $$

where $U$ is the dimensionless vector of conservative variables, $F(U)$ the dimensionless convective fluxes and $F_v(U, \nabla(U))$ the dimensionless diffusive fluxes that write respectively as

$$ U = \begin{bmatrix} \rho \\ \rho \mathbf{u} \\ \rho E \end{bmatrix}, \quad F = \begin{bmatrix} \rho \mathbf{u} \\ \rho \mathbf{u} \otimes \mathbf{u} + P \mathbb{I} \\ (\rho E + P) \mathbf{u} \end{bmatrix} \quad \text{and} \quad F_v = \begin{bmatrix} 0 \\ \frac{1}{Re}\sigma \\ \frac{\mu}{(\gamma-1)RePrM^2}\nabla T \end{bmatrix}. \quad (2.2) $$

Here $\mathbb{I}$ stands for the identity matrix, $\rho$ is the dimensionless density, $\mathbf{u} = [u, v, w]^T$ is the dimensionless velocity vector, $E$ is the dimensionless specific total energy, $P$ is the dimensionless thermodynamic pressure and $\lambda$ is the dimensionless thermal conductivity; $P$ relates to the dimensionless conservative variables through the relationship

$$ P = (\gamma - 1) \left( \rho E - \frac{1}{2} \frac{(\rho \mathbf{u}) \cdot (\rho \mathbf{u})}{\rho} \right), \quad (2.3) $$

and to the dimensionless static temperature $T$ through the following ideal gas dimensionless equation of state:

$$ T = \gamma M^2 \frac{P}{\rho}. \quad (2.4) $$

The dimensionless viscous stress tensor is expressed as

$$ \sigma = \mu (\nabla \mathbf{u} + \nabla^T \mathbf{u}) - \frac{2}{3} \mu (\nabla \cdot \mathbf{u}) \mathbf{I}. \quad (2.5) $$

We assume that the dynamic viscosity only depends on the temperature through Sutherland's law. Given the Prandtl number, the dimensionless thermal conductivity is then deduced from the dimensionless dynamic viscosity $\lambda = \mu C_p / Pr$.

## 2.3. Numerical methods

The DNS of this flow is performed using an in-house finite-volume based DNS and/or large-eddy simulation solver parallelized using the MPI (message passing interface) library. The numerical methods employed are presented in detail in Ben Hassan Saïdi et al. (2020), where the ability of this code to compute high Reynolds compressible (turbulent
and shocked) flows has already been demonstrated (including the boundary conditions presented in § 2.4). An one step (OS) procedure is employed that splits the resolution into the convective part and the viscous part of the problem. The convective fluxes are discretised by a high-order one step (OS) finite-volume scheme based on the Lax–Wendroff approach that is a seventh-order accurate coupled time and space approximation. It is coupled to an original monotonicity-preserving (MP) shock-capturing procedure to form the OSMP7 scheme (first introduced in Daru & Tenaud 2004). A second-order centred finite difference scheme is used for the spatial discretisation of the diffusive fluxes. The temporal discretisation of the diffusive fluxes is obtained by means of a second-order Runge–Kutta time integration. In Ben Hassan Saïdi et al. (2020) we showed that the use of a higher-order scheme for computing diffusive fluxes has a negligible effect on the accuracy of the SWBLI calculation and is therefore unnecessary.

The simulation is performed with a time step of $\Delta t \simeq 1.19 \times 10^{-4} \times L_{int}/U_{\infty}$ corresponding to a maximum Courant–Friedrichs–Lewy number of 0.5 in the computational domain.

## 2.4. Computational domain and boundary conditions

The geometry of the 3-D domain is $D = [0; 250\delta] \times [0; 125\delta] \times [0; 62.5\delta]$. The domain is discretised using a Cartesian mesh with non-uniform spacing in the direction normal to the wall (z). The mesh employs $800 \times 400 \times 202$ grid points in (x × y × z). In the normal-to-the-wall direction, the mesh is tightened close to the wall using a hyperbolic tangent law to obtain a minimum grid spacing over the plate of $\Delta z_{min} = 0.0125\delta$. A grid convergence study (see Appendix B) has been performed in order to assess that the employed mesh is fine enough to compute precisely the flow in the interaction region.

A uniform flow is prescribed as the inlet of the domain as prescribed by the guidelines of the 4th HiOCFD workshop (Wang et al. 2013; HiOCFD 2016) to compute the SWBLI in the configuration of Degrez et al. (1987), which has been validated in Ben Hassan Saïdi et al. (2020) against reference experimental and numerical results. The incident shock wave is created by imposing the Rankine–Hugoniot relationships in the inlet plane at a height z chosen so that the shock wave impinges the flat plate at the desired abscissa $x_{sh} = 62.5\delta$. No-slip and adiabatic wall conditions are prescribed at the flat plate location (z = 0). Outlet time-dependent non-reflecting boundary conditions (Thompson 1987) are imposed at the top surface and at the downstream outlet boundary of the computational domain. Periodic conditions are used in the spanwise direction (y).

The numerical strategy has already been validated for this type of simulation. Indeed, the same numerical schemes and boundary conditions have already been successfully employed to compute SWBLIs in the flat plate configuration. Please refer to Ben Hassan Saïdi et al. (2020) for more details.

## 2.5. Database

After the initial transient, corresponding to the initial propagation of the incident shock wave in the domain and the creation of the separation bubble and subsequent shock wave system, $3 \times 10^6$ iterations have been simulated with a time step of $\Delta t \simeq 1.19 \times 10^{-4} \times L_{int}/U_{\infty}$. A snapshot has been recorded every 1000 time steps so that a total of 3000 snapshots are available for analysis. Each snapshot consists in the state vector of primitive variables $q(\mathbf{x}_i, t_l) = [\rho, u, v, w, T]^T$ evaluated at each mesh point $\mathbf{x}_i$ (with $1 \le i \le M, M$ being the number of mesh points considered) in the 3-D domain containing the interaction region and at each time of snapshot recording $t_l$ (with $1 \le l \le N_t, N_t$ being the number of
<table><caption>Table 2. Sampling characteristics of time signals.</caption><tr><td>Signal length</td><td>$T_{tot} = 562 \times L_{int}/U_{\infty}$</td></tr><tr><td>Time step</td><td>$\Delta t \approx 0.19 \times L_{int}/U_{\infty}$</td></tr><tr><td>Number of samples</td><td>3000</td></tr></table>

snapshots recorded). All data are gathered in the snapshot matrix

$$ \mathbf{Q} = [\mathbf{q}_1 \dots \mathbf{q}_{N_t}] = \begin{pmatrix} \mathbf{q}(\mathbf{x}_1, t_1) \dots \mathbf{q}(\mathbf{x}_1, t_{N_t}) \\ \vdots & \vdots \\ \mathbf{q}(\mathbf{x}_M, t_1) \dots \mathbf{q}(\mathbf{x}_M, t_{N_t}) \end{pmatrix}. \quad (2.6) $$

All the results presented in the rest of this paper are based on the exploitation of the data contained in this snapshot matrix.

All the power spectral densities (PSDs) of signals presented in § 2.6 have been computed using the periodogram method implemented in the scipy.signal package of python (Virtanen et al. 2020) using Hann windowing. For all time signals studied in this section, the sampling characteristics are given in table 2. The total length of the signal considered is around 11 events at $St = 0.02$ and 22 events at $St = 0.04$.

For both SPOD and BSMD analysis algorithms introduced in §§ 3 and 4, the Welch's method (Welch 1967) is used to compute the discrete Fourier transform (DFT) of the state vector. The matrix of snapshots is split into $N_b = 10$ with an overlapping of 67 % between the blocks. With this slicing, the total length of the signal considered for the analysis is around 30 events at $St = 0.02$ and 60 events at $St = 0.04$.

In the following, for both SPOD and BSMD analyses, we consider the centred snapshot matrix

$$ \mathbf{Q}^c = [\mathbf{q}_1^c \dots \mathbf{q}_{N_t}^c] = \begin{pmatrix} \mathbf{q}^c(\mathbf{x}_1, t_1) \dots \mathbf{q}^c(\mathbf{x}_1, t_{N_t}) \\ \vdots & \vdots \\ \mathbf{q}^c(\mathbf{x}_M, t_1) \dots \mathbf{q}^c(\mathbf{x}_M, t_{N_t}) \end{pmatrix}, \quad (2.7) $$

with $\mathbf{q}^c(\mathbf{x}_i, t_l) = \mathbf{q}(\mathbf{x}_i, t_l) - \bar{\mathbf{q}}(\mathbf{x}_i)$, where $\bar{\mathbf{q}}(\mathbf{x}_i)$ is the vector of primitive variable averages in time and spanwise direction.

Note that, for the signal length, for the following SPOD and BSMD analyses, a convergence study has been performed. Indeed, these analyses have been performed with different numbers of snapshots (i.e. a different length of the signal and frequency resolution) leading to the same conclusions to those presented below.

## 2.6. Mean flow organisation and frequency content

The time mean longitudinal velocity averaged in the spanwise direction $\bar{u}$ is plotted in figure 3 in the interaction zone. The velocity is scaled by the free-stream velocity. The red line shows the isocontour $\bar{u} = 0$. It allows us to clearly identify the separation bubble in which the mean longitudinal velocity is negative. The separation bubble also corresponds to the region in which the mean skin friction coefficient plotted in figure 4 is negative. The green line in figure 3 represents the $z/\delta = 2.5$ plane in which we analyse the characteristic sizes of the flow structures in the transverse direction in figure 6(a). The green dots, referenced by numbers 1–4, mark the position of the probes placed in the mixing layer to analyse its characteristic frequencies by calculating the PSD of the velocity signal shown in figure 7.
<figure><img src="image_5.png" /><figcaption>Figure 3. Time mean longitudinal velocity field averaged in the spanwise direction $\bar{u}$. The red line shows the isocontour $\bar{u} = 0$. The green line is parallel to the wall at a height of $z/\delta = 2.5$. Coordinates of probes in the $(x, z)$ basis, from 1–4: (13.44$\delta$, 0.3$\delta$), (43.75, 3.5$\delta$), (68.75$\delta$, 2.5$\delta$), (92.50$\delta$, 0.3$\delta$).</figcaption></figure>

As the incident boundary layer is perfectly laminar, it is particularly prone to separation when subjected to the adverse pressure gradient imposed by the incident shock wave, compared with turbulent boundary layers. The resulting mean separation length is $L = \bar{x}_R - \bar{x}_S \approx 74.68\delta$, with $\bar{x}_S \approx 13.44\delta$ and $\bar{x}_R \approx 88.12\delta$ respectively the time mean separation and reattachment points averaged in the spanwise direction. The interaction length is $L_{int} = x_{sh} - \bar{x}_S = 49.06\delta$. The time mean height of the bubble, averaged in the spanwise direction, is $\bar{h} = 3.5\delta$. Moreover, the mean separation bubble is quasi-symmetric with respect to the vertical axis passing through the apex of the bubble. This length and aspect of the mean bubble is consistent with comparable simulations (Robinet 2007; Song & Hao 2023). For similar Mach and Reynolds numbers, the length of the separation bubbles for SWBLIs with laminar incident boundary layers is around one order of magnitude bigger than their counterpart for SWBLIs involving turbulent incident boundary layers; see, for example, Adler & Gaitonde (2018) where $\bar{L} = 4.0\delta$ or Ben Hassan Saïdi et al. (2020) where $\bar{L} = 4.92\delta$. For turbulent SWBLIs, the symmetry of the separation bubble is also broken, the second part of the bubble being less long than the first part, which is due to the increased mixing rate in the shear layer. The same symmetry breaking (less marked) is observed for SWBLIs with a sufficiently strongly forced laminar incident boundary layer (Sansica 2015; Diop et al. 2019; Mauriello et al. 2022).

The difference of topology of the mean separation bubble of the present flow configuration with respect to turbulent SWBLIs calls for a comment about the characteristic length used to build non-dimensionalised frequencies. Indeed, the main frequencies involved in the interactions have been introduced in the literature for turbulent interactions as Strouhal numbers based on the separation length $St_L = fL/U_e$. For the present flow configuration, as the bubble is symmetric, the interaction length $L_{int}$ is more suitable to be consistent with the characteristic frequencies highlighted in turbulent interactions. In the following, the frequencies will therefore be expressed as Strouhal numbers based on the interaction length $St_{L_{int}} = fL_{int}/U_e$.

A snapshot of the flow is shown in figure 5. The discriminant criterion, introduced in Chong et al. (1998), already used by Pirozzoli & Grasso (2006) in the context of SWBLI simulations, is used to identify the vortex structures present in the flow. Shock waves are highlighted by isosurfaces of $||\nabla P||$ coloured in black. The upstream boundary layer not being artificially forced by modes resulting from a stability problem (Sansica et al. 2014, 2016) or a free-stream turbulence method, the only fluctuations supplied to the flow come from shock waves (leading edge, separated or the incident shock waves) that disrupt the flow mainly two-dimensionally, although the natural most unstable normal modes are 3-D.
<figure><img src="image_6.png" /><figcaption>Figure 4. Distribution of the spanwise averaged skin friction along the flat plate (____). Blasius laminar boundary layer solution (____). Turbulent correlation from Cousteix (1989) (---).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5. Vortical structures highlighted by isosurfaces of the discriminant criterion coloured by the magnitude of the longitudinal component of the velocity. Shock waves are highlighted by isosurfaces of $\|\nabla P\|$. The upper left insert shows a zoomed view of the vortical structures around the separation bubble. The upper right insert shows an overhead view of the vortical structures in the separation bubble and in the downstream boundary layer. On the upper right insert, the vertical black lines indicate, from left to right, the mean separation line (—), the mean line of incident shock impingement on the apex of the separation bubble (---), the mean line of shock impingement on the flat plate (......) and the mean reattachment line (---).</figcaption></figure>

Figure 5 shows clearly this characteristic. In what we seek to show, the existence of a nonlinear interaction between medium-frequency modes generating low-frequency modes, this has no consequence whether the convective modes exited are 2-D or 3-D. Large spanwise vortices, corresponding to Kelvin-Helmholtz rolls that progressively develop, are clearly visible in the shear layer edging the separation bubble between the reflected and the incident shocks. After the incident shock impingement, the shear layer is populated by elongated structures in the streamwise direction. Inside the separation bubble, 3-D structures are visible in the downstream part of the separation zone. For the particular value of the discriminant criterion shown in figure 5, no coherent structure is visible in the early part of the interaction. The dynamical activity inside the separation bubble therefore
<figure><img src="image_8.png" /><figcaption>Figure 6. Spanwise length scales in the interaction zone. (a) Time averaged spanwise wavenumber of the longitudinal component of velocity *u* in a plane parallel to the flat plate at a height of *z*/*δ* = 2.5. The vertical dashed lines indicate the *β* = 0.25 and *β* = 2 wavenumbers (_______). The horizontal dashed lines indicate, from bottom to top, the separation point (___), the raising shear layer (_______), the descending shear layer ( ) and the reattachment point (___). (b) Isocontours of the longitudinal component of velocity *u* = −0.041 coloured by *u*. Slice (y plane) of the density *ρ*.</figcaption></figure>

seems to be mainly concentrated in its second part. The dynamic activity of the separation bubble causes the rapid transition of the boundary layer after the reattachment point. The turbulent nature of the boundary layer downstream of the interaction zone is confirmed by the downstream distribution of the mean skin friction coefficient, as visible in figure 4.

In order to physically characterise the longitudinal structures identified in figure 5, the spectral content of the velocity signal in the transverse direction is studied. To this end, the time averaged one-dimensional spatial Fourier transform in the spanwise direction of the longitudinal component of the velocity (u) was calculated in different planes parallel to the wall, defined by z = cste (z/δ = 0.3, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5). Figure 6(a) shows the evolution of these spectra as a function of the abscissa along the flat plate in the z/δ = 2.5 plane marked by the green line in figure 3. The wavenumber is scaled using the same scaling as in Bugeat et al. (2022) for comparison. The abscissa is normalised by using the mean interaction length Lint and the mean separation point x̄S. In this figure, the cyan and green horizontal dashed lines indicate the intersection between the plane under consideration and the shear layer bounding the separation bubble. The line at (x - x̄S)/Lint = 0.445 indicates the shear layer upstream of the impact of the incident shock, which we call the ‘rising shear layer’. The line at (x - x̄S)/Lint = 1.2 indicates the shear layer downstream of the impact of the incident shock, which we call the ‘descending shear layer’. The vertical dashed lines indicate the β = 0.25 and β = 2 wavenumbers that were highlighted as the two maxima of the G(β) curve, where G is the optimal gain of a low-frequency resolvent analysis, and that were respectively associated with Görtler-type structures and streaks in Bugeat et al. (2022). Just upstream of the descending shear layer (i.e. inside the bubble) we can clearly see a peak at β = 0.25. In the descending shear layer, we see a plateau containing the value β = 2. These same peaks are found in the spectra of the other planes (not shown here), indicating the presence of Görtler-type structures inside the bubble and streaks in the descending shear layer. This location of Görtler-type structures and streaks is in agreement with the results of Bugeat et al. (2022). The optimal modes of low-frequency resolvent analysis are therefore found in the DNS data of the present unforced transitional SWBLI. Analysis of the instantaneous 3-D fields highlights
these longitudinal structures and illustrates their location. For example, figure 6(b) shows
the isocontour of the longitudinal velocity field $u = -0.041$ next to a 2-D slice coloured
by $\rho$ to locate the separation bubble. The wavelength of the Görtler-type structures and
the streaks clearly manifest themselves in the isocontour of velocity in the zones identified
by the spectral analysis.

We now want to document the frequency content of the dynamics of the interaction
zone. The dynamic structures mainly visible in figure 5 seem to be concentrated in the
shear layer. We therefore begin by examining the frequency content of the mixing layer.
To do this, we examine the PSD of the longitudinal velocity signal measured at different
points in the shear layer. These points are shown in a 2-D plane in figure 3, where their
coordinates in the (x, z) plane are specified. In practice, for these coordinates in the (x, z)
plane, a probe was placed at each mesh point in the y direction to measure the velocity
signal and compute the PSD. The PSDs were then averaged over span.

We can see that in the zone of the separation point (probe 1), the energy level of
the oscillations in the shear layer is very low and localised in the low- and medium-
frequency ranges characterising the breathing and flapping phenomena. There are clear
peaks at frequencies $St_{L_{int}} = 0.04, 0.0667$ and $0.1$. This frequency content is amplified
in the rising mixing layer. Indeed, the signal measured by probe 2 has approximately the
same frequency content as probe 1, but at a higher energy level.

We saw above that the passage through the incident shock is associated with a distortion
of the vortex structures, which are mainly oriented in the spanwise direction upstream of
the shock and preferentially in the streamwise direction downstream of the shock. These
structures are rapidly amplified downstream of the incident shock. Just downstream of
the incident shock (probe 3), the energy of the fluctuations is amplified compared with
the level observed at probe 2. The peaks measured at low and medium frequencies are
still present, but the spectrum is filling out and these frequencies are less predominant. In
addition, higher frequencies appear. This signal is massively amplified in the descending
shear layer. In the attachment zone (probe 4), the energy of the signal is almost quadrupled
compared with probe 3. This amplification is also selective. The frequencies already
highlighted, $St_{L_{int}} = 0.04, 0.0667$ and $0.1$, are particularly peaked in the reattachment
zone. There is also a further filling in of the spectrum at high frequencies concomitant with
the generation of fine structures associated with the transition to turbulence in this region.

These measurements in the shear layer seem to indicate that the mixing layer is forced
from the separation point by the characteristic SWBLI frequencies at low and medium
frequencies. The rising shear layer amplifies these frequencies by a factor of around
10. However, the energy levels in the rising shear layer are low relative to those in the
descending shear layer. The latter shows a massive amplification of fluctuations in the
low- and medium-frequency ranges, and in particular, of the frequencies characteristic of
the SWBLI unsteadiness $St_{L_{int}} = 0.04, 0.0667$ and $0.1$. These observations are consistent
with various previous works, in particular, those of Sansica et al. (2014) and Mauriello
et al. (2022) who have shown that the low-frequency content of the interaction zone is
created by the dynamics of the shear layer and seems to result from the nonlinear evolution
of the latter in the downstream part of the recirculation bubble. These authors also showed
that disturbances at these frequencies then travel upstream inside the recirculation bubble.

In order to study the frequency content of this information feedback in the recirculation
bubble, we plot in figure 8(a) the evolution of the spanwise averaged premultiplied PSD
of the longitudinal velocity signal ($u$) measured in a plane parallel to the flat plate, in
very close proximity to the wall ($z/\delta = 0.3$). For each abscissa along the flat plate, the
PSD at this abscissa is normalised and premultiplied. This makes it possible to determine
the predominant frequencies for each abscissa, but the relative energy differences between
<figure><img src="image_9.png" /><figcaption>Figure 7. Power spectral densities of the spanwise averaged longitudinal velocity $u$ signal recorded on probes shown in figure 3. (a) Probe 1. (b) Probe 2. (c) Probe 3. (d) Probe 4.</figcaption></figure>

abscissas are erased by this representation. (The relative difference in energy between abscissa is however shown in figure 8(b), which is commented on below.) In this figure, the horizontal dashed lines indicate, from bottom to top, the positions of the separation point, the rising shear layer, the impact of the incident shock on the flat plate and the reattachment. We can see that the spectral content in the attachment zone is present throughout the bubble in the near wall, but with a progressive narrowing of the significant frequency range towards the low- and medium-frequency range as we move upstream in the bubble. Thus, the disturbances reaching the foot of the rising shear layer have an exclusively low- and medium-frequency spectral content with a cutoff frequency of approximately $St_{Lint} \approx 0.1$. We can therefore see that the rise of the disturbance inside the bubble selects the low and medium frequencies. This result had already been shown in Mauriello et al. (2022). In order to determine whether this feedback into the bubble is damped, we integrate the PSD between the $St_{Lint} = 0.02$ and 0.1 for each abscissa along the flat plate to obtain the streamwise evolution of the expected power in this frequency range $E(x) = \int_{St=0.02}^{St=0.1} \text{PSD}(u(x))dSt$. This power distribution is shown in figure 8(b). This figure clearly shows a large peak of power in the reattachment region, in agreement with the analysis of the spectra shown in figure 7. We can clearly see an exponential decrease in power inside the bubble. The upwelling of disturbances is therefore damped inside the bubble. The signal reaches the detachment zone with low but finite energy levels. The shear layer is therefore forced into the separation zone by a signal of finite energy and quasi-continuous spectrum in the frequency range $St_{Lint} \in [0.02, 0.1]$. However, figure 8(a) shows a high selectivity of the mixing layer. Indeed, the quasi-continuous spectrum
<figure><img src="image_10.png" /></figure>

Figure 8. Power spectral density and power of the longitudinal velocity signal close to the wall inside the separation bubble. (a) Distribution along the flat plate of the spanwise averaged premultiplied and normalised PSD of the longitudinal component of velocity *u* measured in a plane parallel to the flat plate at a height of $z/\delta = 0.3$. The vertical dashed lines indicate significant frequency peaks in the low- and medium-frequency range (______)_. The horizontal dashed lines indicate, from bottom to top, the reflected shock foot (_______), the separation point (____), the crossing of the rising shear layer (_______), the incident shock impingement location (_______) and the reattachment point (____). (b) Distribution along the flat plate of the expected power in the range $St_{Lint} \in [0.01, 0.1]$ of the longitudinal component of velocity *u* measured in a plane parallel to the flat plate at a height of $z/\delta = 0.3$. The vertical dashed lines indicate, from left to right, the separation point (____), the crossing of the rising shear layer (_______), the incident shock impingement location (_______), the crossing of the descending shear layer (____) and the reattachment point (____).

becomes strongly peaked in the shear layer at the frequencies already highlighted in figure 7(a), i.e. $St_{Lint} = 0.04, 0.0667$ and $0.1$.

We are also interested in the pressure forces exerted on the wall in the interaction zone. To this end, we analyse the evolution in the streamwise direction of the premultiplied PSD normalised to the wall pressure shown in figure 9(a). A similar evolution of the spectrum can be observed between the points of separation and reattachment in the sense that we move from the low/medium-frequency range to a wider range, extended to higher frequencies. However, we can clearly see that in the first part of the separation bubble, the wall pressure spectrum predominantly highlights the spectral signature of the shear layer, which is peaked at frequencies $St_{Lint} = 0.04, 0.0667$ and $0.1$. In both respects, the wall pressure spectrum is very similar to the spectra presented in Mauriello et al. (2022). In this simulation, we therefore record pressure oscillations at the foot of the reflected shock dominated by the low- and medium-frequency range $St_{Lint} \in [0.02, 0.1]$ that characterises the SWBLI unsteadiness. The evolution in the streamwise direction of the power of the parietal pressure signal in this frequency range is shown in figure 9(b). As with the near-wall velocity signal, there is a quasi-exponential attenuation of the wall pressure perturbations as we move upwards from the reattachment point. We then see an increase in the first quarter of the bubble as we get closer to the separation point and the shear layer moves closer to the wall. The signal power in the separation zone is an order of magnitude lower than the power in the reattachment zone.

The complex dynamics described above are also reflected in the breathing and flapping cycles of enlargement and shrinkage of the separation bubble, which will be described in § 3 documenting the SPOD analysis we carried out. This low- and medium-frequency dynamics of the bubble can also be seen by analysing the temporal evolution of the position of the separation and reattachment points, respectively, $x_S$ and $x_R$, shown in figure 10. It can be seen that the reattachment point is subject to oscillations whose maximum
<figure><img src="image_11.png" /><figcaption>Figure 9. Spanwise averaged PSD and power of the wall pressure along the flat plate. (a) Distribution along the flat plate of the spanwise averaged premultiplied and normalised PSD of the wall pressure. The vertical dashed lines indicate significant frequency peaks in the low- and medium-frequency range. The horizontal dashed lines indicate, from bottom to top, the reflected shock foot (____), the separation point (---), the incident shock wave impingement location (____) and the reattachment point (---). (b) Distribution along the flat plate of the expected power in the range $St_{Lint} \in [0.02, 0.1]$ of the wall pressure. The vertical dashed lines indicate, from left to right, the separation point (----), the incident shock impingement location (---) and the reattachment point (---).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10. History of the spanwise averaged abscissa along the flat plate of (a) the separation point $x_S$ and (b) the reattachment point $x_R$.</figcaption></figure>

extrusions correspond to about $7\delta$ while the separation point is subject to very weak oscillations with maximum amplitudes of the order of $0.4\delta$. This order-of-magnitude difference in the amplitudes of the oscillations of $x_S$ and $x_R$ is consistent with the order-of-magnitude difference observed in the power of the oscillations of the parietal pressure and velocity signals in these two regions. Similarly, in agreement with the wall pressure and velocity signals in these regions of the flow, the spectral content of $x_S(t)$ oscillations is mainly located in the low- and medium-frequency range with significant peaks at frequencies $St_{Lint} = 0.04, 0.0667$ and $0.1$, while the oscillations of $x_R(t)$ present a wider band spectrum with peaks also localised at these frequencies, as shown in figure 11.

In the present flow configuration, the low energy level of perturbations coming from the reattachment point reaching the separation point explains the low amplitude of the separation point and associated reflected shock foot (not shown here) oscillations.
<figure><img src="image_13.png" /><figcaption>Figure 11. Spanwise averaged PSDs of (a) the separation $x_s$ and (b) the reattachment point locations $x_R$.</figcaption></figure>

This point is a salient feature of the flow studied here, which makes it different from configurations with a forced laminar or fully turbulent boundary layer. As already stated in the introduction of this paper, in this respect, the present case study represents a limit case in which the mechanisms at the origin of the low-frequency dynamics appear to be present at low intensity and are not disturbed by the natural frequencies of the incident turbulence that, we believe, makes it a good candidate for identifying the intrinsic underlying mechanisms of this dynamics. In the cases of forced or turbulent incoming boundary layers, the more intense dynamics of the mixing layer (due to upstream forcing) has two consequences: (i) the energy of the perturbations at reattachment, and therefore, the energy of the signal emitted from reattachment is much higher; (ii) the separation length, i.e. the distance between the zone of signal emission and the separation point, is smaller. The combination of these two characteristics means that the disturbances received at the separation point are more energetic than in the case studied in this paper. In the forced and turbulent cases, disturbances are then more likely (for a given strength of interaction) sufficiently energetic to induce higher amplitude oscillations at the separation point.

### 3. Spectral proper orthogonal decomposition analysis

In order to get more insight in the low/medium-frequency dynamics of the flow, we perform a modal decomposition of the database formed by the centred snapshot matrix introduced in the previous section. In particular, we use SPOD, introduced in Towne et al. (2018) and Schmidt & Colonius (2020) (to which we refer to for more details about the methodology), to highlight the physical structure of the flow at each frequency. Indeed, SPOD can be schematically described as a POD applied in the frequency domain frequency by frequency. The first step of the algorithm is the computation of the DFT of the state vector. To this end, the matrix of snapshots (2.7) is split into $N_b$ segments of $N_f$ snapshots with an overlapping of $N_o$ snapshots; thus, $N_b$ snapshot matrices $Q^{(n)} = [\hat{q}_1^{(n)} \dots \hat{q}_{N_f}^{(n)}]$ are obtained, with $1 \le n \le N_b$. The DFT is applied to each snapshot matrix to build $N_b$ snapshot matrices in the frequency domain $\hat{Q}^{(n)} = [\hat{q}_1^{(n)} \dots \hat{q}_{N_w}^{(n)}]$, with $N_w = N_f/2 + 1$ the number of resolved frequencies. As a second step of the algorithm, the snapshot matrices in the spectral domain $\hat{Q}_k^{(n)}$ are sorted by frequency to extract the matrices of realisations by frequency $\hat{Q}_k^c = [\hat{q}_k^{c(1)} \dots \hat{q}_k^{c(N_b)}]/\sqrt{N_b}$, with $1 \le k \le N_w$. As already stated in § 2.5, in our case, $N_b = 10$ and the overlapping between
blocks is 67 %. The third step consists in performing a POD to each matrix of realisations $\hat{\mathbf{Q}}_k$. In our implementation, we use the snapshot POD. In the frequency domain, it consists in solving the eigenvalue problem

$$ \mathbf{C}_k \Theta_k^H = \Theta_k^H \Lambda_k, \quad (3.1) $$

where $\mathbf{C}_k = \hat{\mathbf{Q}}_k^{cH} \mathbf{W} \mathbf{Q}_k^c$ is the matrix of spectral densities, $\Theta_k \in \mathbb{C}^{N_b \times N_b}$ is the projection matrix in which each column contains the projection coefficients of one realisation over the spatial modes, $(.)^H$ is the conjugate transpose and $\Lambda_k \in \mathbb{R}^{N_b \times N_b}$ is the diagonal matrix containing the eigenvalue associated with each spatial mode with $\lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_{N_b} \ge 0$. In our implementation, the inner product is weighted using the diagonal matrix of spatial quadrature weights $\mathbf{W}$ introduced by Chu (1965). The matrix of spatial modes is then obtained by projection

$$ \psi_k = \hat{\mathbf{Q}}_k^c \Theta_k^H. \quad (3.2) $$

Each column of $\psi_k$ contains the $m$th spatial mode $\psi_m(f_k)$ associated with $\lambda_m$.

As a result, each spectral realisation $\hat{\mathbf{q}}_k^{(n)}$ at frequency $f_k$ is decomposed into a linear combination of $N_b$ spatial modes $\psi_m$ as

$$ \hat{\mathbf{q}}^{(n)}(\mathbf{x}_i, f_k) = \sum_{m=1}^{N_b} a_m^{(n)}(f_k) \psi_m(\mathbf{x}_i, f_k), \quad (3.3) $$

where the $a_m^{(n)}$ coefficients are the components of the projection coefficients at each frequency.

The spatial modes correlate space and time and they are ordered by their probability of existence in the flow quantified by their associated eigenvalue.

The knowledge of the $m$th spatial mode at frequency $f_k$, $\psi_m(\mathbf{x}_i, f_k)$, allows the construction of the associated $m$th spatio-temporal mode by inverse Fourier transform

$$ \Phi_m(\mathbf{x}_i, t) = \psi_m(\mathbf{x}_i, f_k) e^{2i\pi f_k t}. \quad (3.4) $$

The spatio-temporal modes represent structures that evolve coherently in space and time and are ordered according to the most likely dynamics for a given frequency (mode 1 being the most probable one, then mode 2, etc.).

## 3.1. Results

The SPOD spectrum of eigenvalues $\lambda_m$ is shown in figure 12. We can see that the SPOD spectrum describes a marked dynamic in the low- and medium-frequency range, with a clear drop-off beyond that. We can clearly see that the energy decreases with the number of the mode, so that the first mode contains about twice as much energy as the third. In the low- and medium-frequency range, almost 60 % of the flow energy is contained in the first two modes. Interestingly, the three main frequencies for mode 1 are the frequencies found in local DNS measurements (see figure 7): $St_{L_{int}} = 0.04, 0.0667$ and $0.1$. In particular, we have already shown that these three frequencies, which are present throughout the shear layer, are largely dominant in the rising shear layer. We can also see that mode 2 signs the frequency $St_{L_{int}} = 0.04$ and its first subharmonic $St_{L_{int}} = 0.02$ as well as the sum of the two ($St_{L_{int}} = 0.06$) and another frequency at $St_{L_{int}} = 0.073$. The three frequencies highlighted by the second mode (i.e. $St_{L_{int}} \simeq 0.02, 0.06$ and $0.073$) are well present in the spectra of the DNS signals measured locally in the shear layer shown in figure 7. However, unlike the frequencies highlighted by the first mode, these frequencies reach energy levels
<figure><img src="image_14.png" /><figcaption>Figure 12. The SPOD spectrum of the state vector $\hat{\boldsymbol{q}}$ for each mode ($m_j, j = 1, \dots, 10$) expressed as a relative energy given in percent.</figcaption></figure>

comparable to the former in the descending shear layer, downstream of the impact of the incident shock, and are comparatively less energetic in the upstream part of the separation bubble. The first mode therefore seems to bring out the dominant frequencies in the whole of the interaction zone, while the second mode brings out frequencies that are mainly present in the downstream part of the interaction zone.

We are now interested in the shape of the dominant modes in the entire bubble. As an illustration, figure 13 shows the modulus of the first mode at frequency $St = 0.04$ (characteristic frequency of the separation bubble breathing) for the three velocity components $u, v$ and $w$. For the longitudinal component $u$, we can see that the activity of the mode is mainly localised in the shear layer. In particular, large amplitudes of the mode are observed in the rising shear layer. In the descending shear layer, the amplitude increases the closer we get to the abscissa of the reattachment point. In the attachment zone, strong amplitudes are also observed close to the wall and amplify downstream of the interaction. We also observe activity in the longitudinal velocity mode localised on the reflected shock as well as the reattachment shock, materialising the oscillations of these shocks at this low frequency, characteristic of SWBLI unsteadiness. These oscillations of the reflected and reattachment shocks are even more visible in the shape of the vertical velocity mode $w$, whose activity is localised in these zones as well as in the shear layer, with less intensification than for the longitudinal velocity in the reattachment zone. The activity of the mode for the spanwise velocity component $v$ is completely localised inside the separation bubble.

The spatial modes at frequencies $St_{Lint} = 0.0667$ and $0.0934$ are qualitatively almost identical to those at frequency $St_{Lint} = 0.04$. For the sake of brevity, they are shown in Appendix A, in figures 29 and 32 respectively for $St_{Lint} = 0.0667$ and $St_{Lint} = 0.0934$.
<figure><img src="image_15.png" /><figcaption>Figure 13. First spatial SPOD mode ($m_1$) averaged in the spanwise direction at $St = 0.04$. The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$. (b) Spanwise velocity $v$. (c) Vertical velocity $w$.</figcaption></figure>

In figure 14 we plot the real part of the first spatio-temporal mode $\Phi_1(x_i, t)$ of the three velocity components for the Strouhal number $St = 0.04$, characteristic of the breathing of the recirculating bubble. We plot regularly spaced snapshots in a period $T = 1/St$. For each speed component, a supplementary movie of the animation of the spatio-temporal evolution of the associated mode is available at https://doi.org/.

In the longitudinal velocity mode $u$, we can clearly see that the oscillations of the reflected shock in the longitudinal direction are in phase with the oscillations of the shear layer in this direction, and therefore, in phase with the oscillations of the separation point. We can also see that the breathing cycles follow the following sequence: the value of the mode changes in the downstream part of the bubble before being propagated upstream in the first part of the bubble, eventually contaminating the entire shear layer. The upstream
<figure><img src="image_16.png" /><figcaption>Figure 14. Visualisation of $\Phi_1(x, t)$ velocity $u$ at $St = 0.04$, averaged in the spanwise direction. (a) Velocity component $u$. (b) Velocity component $v$. (c) Velocity component $w$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15. Representation of $\Phi_1(x, z=\text{cste}, t)$ for $z/\delta = 1.48$. (a) Velocity component $u$. (b) Velocity component $v$. The vertical black lines indicate, from left to right, the separation point, the crossing of the rising shear layer, the crossing of the incident shock wave, the crossing of the descending shear layer and the reattachment point.</figcaption></figure>

propagation of information in the bubble is shown in figure 15, where we plot versus time (for a single breathing period) the value of $\Phi_1^u$ along the flat plate measured at a constant height above it $z/\delta = 1.48$. In the first part of the bubble, we clearly distinguish regions of the same sign oriented with a negative slope in the $(x, t)$ frame, denoting an upstream
propagation of this sign. This behaviour is qualitatively observed at all altitudes in the separation bubble. The upstream propagation from the latter part of the separation bubble is also visible in the spanwise velocity mode as shown in figures 14(b) and 15(b). This result seems to confirm that the downstream part of the separation plays a decisive role in the breathing dynamics of the separation bubble. More precisely, it confirms the analysis of § 2.6: the low-frequency oscillations of the separation bubble are driven by oscillations at this frequency in its latter part (materialised by the change in sign of the mode in the reattachment zone) propagated upstream in the separation bubble.

A similar qualitative behaviour is observed in the results obtained for the frequencies $St_{Lint} = 0.0667$ and $St_{Lint} = 0.0934$, the latter being characteristic of the flapping of the separation bubble. The figures spatial evolution of the first mode at these frequencies can be found in Appendix A, in figures 30, 31 and 33, 34 respectively for $St_{Lint} = 0.0667$ and $St_{Lint} = 0.0934$. The similarity of the spatio-temporal modes for breathing and flapping may suggest a link to be determined between these two oscillation modes.

## 4. Bispectral mode decomposition analysis

In order to study possible nonlinear links between modes at different frequencies, the BSMD algorithm has been applied to highlight possible triadic interactions taking place in the interaction zone. Triadic interactions result from the quadratic nonlinearities of the Navier-Stokes equations. They are the fundamental mechanism of energy transfer in fluid flows, and manifest themselves in the frequency domain as triplets of three frequencies ($f_1, f_2, f_3$), the combination of which is zero, i.e. $f_1 \pm f_2 \pm f_3 = 0$. Such interactions can be identified from triple correlations of frequency components, i.e. from the bispectrum, which is a measure of the quadratic nonlinearities at the bifrequency ($f_1, f_2$). For a scalar signal $q$, it is defined by

$$S_{qqq}(f_1, f_2) = \lim_{T \to \infty} \frac{1}{T} E[\hat{q}(f_1)^* \hat{q}(f_2)^* \hat{q}(f_1 + f_2)], \quad (4.1)$$

where $(.)^*$ denotes the complex conjugate and $E$ the expectation operator.

### 4.1. The BSMD methodology

The BSMD algorithm is a modal decomposition that reveals the presence of triadic nonlinear interactions from multidimensional data. To this end, an integral measure of the pointwise bispectral density between frequency components of the vector of primitive variables is introduced as

$$b(f_k, f_l) = E \left[ \int_{\Omega} \hat{q}_k^* \circ \hat{q}_l^* \circ \hat{q}_{k+l} \, dx \right] = E \left[ \hat{q}_{k \circ l}^H W \hat{q}_{k+l} \right] = E \left[ (\hat{q}_{(k \circ l)}, \hat{q}_{(k+l)}) \right], \quad (4.2)$$

where $\hat{q}_{(k \circ l)} = \hat{q}(x, f_k) \circ \hat{q}(x, f_l)$ (with $\circ$ is the Hadamard product), $\hat{q}_{(k+l)} = \hat{q}(x, f_k)$ and $(.)^H$ is the conjugate transpose; $W$ is the diagonal matrix of spatial quadrature weights and $\Omega$ is the spatial domain over which the flow is defined.

Using the Welch's slicing of the snapshot matrix, two modal projections are constructed: the cross-frequency fields

$$\Phi_{k \circ l}^{[i]}(\mathbf{x}, f_k, f_l) = \sum_{j=1}^{N_b} a_i^{[j]}(f_{k+l}) \hat{q}_{(k \circ l)}^{[j]} \quad (4.3)$$
and the bispectral modes

$$ \mathbf{\Phi}_{k+l}^{[i]}(\mathbf{x}, f_{k+l}) = \sum_{j=1}^{N_b} a_i^{[j]}(f_{k+l}) \hat{\mathbf{q}}_{(k+l)}^{[j]}, \quad (4.4) $$

which share a common set of expansion coefficients $a_i^{[j]}$ with $1 \le j \le N_b$.

The goal of BSMD is to compute modes that optimally represent the data in terms of the integral bispectral density defined in (4.2). That is, we seek the set of expansion coefficients $a_1$ that maximises the absolute value of $b(f_k, f_l)$ as defined in (4.2) with $\hat{\mathbf{q}}_{(k+l)}$ and $\hat{\mathbf{q}}_{(k+l)}$ replaced by their modal expansion. We therefore target the set of expansion coefficients that verifies

$$ \mathbf{a}_1 = \arg \max_{\|\mathbf{a}_1\|=1} |E\left[\mathbf{\Phi}_{k\circ l}^{[1]} \mathbf{W}\mathbf{\Phi}_{k+l}^{[i]}\right]|, \quad (4.5) $$

where the coefficient vector is required to be a unit vector in order to guarantee boundedness of the expansion.

More details about the BSMD methodology can be found in Schmidt (2020), especially about the construction and solving of this optimisation problem.

As a result of the BSMD analysis, we obtain optimal cross-frequency fields and bispectral modes, respectively $\Phi_{k\circ l}^{[1]}$ and $\Phi_{k+l}^{[1]}$, as well as the mode bispectrum

$$ \lambda_1(f_k, f_l) = E \left[ \Phi_{k\circ l}^{[1]H} \mathbf{W} \Phi_{k+l}^{[1]} \right] \in \mathbb{C}. \quad (4.6) $$

The bispectrum mode is the fundamental result of the BSMD analysis. Indeed, significant triadic interactions are indicated by local maxima of the modulus of this quantity, also called bispectrum mode amplitude. Bispectral modes $\Phi_{k+l}^{[1]}$ are linear combinations of Fourier modes and can be interpreted as observable physical structures resulting from quadratic interactions between frequencies $f_k$ and $f_l$. The multiplicative cross-frequency fields $\Phi_{k\circ l}^{[1]}$, on the contrary, are maps of phase alignment between two frequency components that may not directly be observed. Once significant triadic interactions are identified in the map of bispectrum mode amplitude, interaction maps can be computed as

$$ \Psi_{k,l}(\mathbf{x}, f_k, f_l) = |\Phi_{k\circ l}^{[1]} \circ \Phi_{k+l}^{[1]}|. \quad (4.7) $$

This field indicates the location of the triadic interaction involving the triad $(f_k, f_l, f_{k+l})$ in the flow field, as it quantifies the average local bicorrelation between the frequencies in the domain.

To sum up, the BSMD analysis process involves calculating the map of mode bispectrum amplitude $\|\lambda_1(f_k, f_l)\|$ to identify significant quadratic interactions. These are the interactions that contribute most to the flow dynamics. Once these interactions have been identified, the associated interaction map $\Psi_{k,l}(\mathbf{x}, f_k, f_l)$ is studied. This allows us to identify the regions of the flow where this triadic interaction occurs. Finally, the associated bispectral mode $\Phi_{k+l}^{[1]}$, which is the physical mode resulting from this interaction, is visualised and interpreted.

In the following we present the BSMD analysis of the database introduced in § 2.5. The interpretation of the results will be performed following the process that has just been outlined above, focusing on the low/medium-frequency dynamics of the interaction. In particular, we will be interested in identifying and documenting nonlinear links between medium frequencies and low frequencies.
I. Ben Hassan Saïdi, S. Wang, G. Fournier, C. Tenaud and J.-C. Robinet

## 4.2. Mode bispectrum amplitude

The modulus of the mode bispectrum $||λ₁(St_k, St_l)||$ obtained when applying BSMD to the snapshot matrix (2.7) is shown in figure 16. Only the region of non-redundant information, as demonstrated by Schmidt (2020), is shown. Moreover, the interactions between under-resolved frequencies for which two or less periods are resolved through the Welch Fourier transform estimation have been greyed out on the graph. The region above the horizontal axis of the graph represents sum interactions $St_1 + St_2 = St_3$, the region below the horizontal axis represents difference interactions $St_1 - St_2 = St_3$. The interactions located on the horizontal axis correspond to interactions for which $St_2 = 0$ resulting in $St_3 = St_1$ and, therefore, correspond to the linear evolution of the flow around its mean. In our analysis, we do not consider local maxima located on the line $St_1 = -St_2$ as they result in modes of frequency $St_3 = 0$ corresponding to the mean flow distortion.

We can clearly identify strong quadratic couplings (highlighted by the high value of $||λ₁(St_k, St_l)||$) for $(St_1, St_2) \in [0, 0.20]^2$, namely in the low- and medium-frequency ranges. The intensity of interactions decreases sharply at higher frequencies. In the low- and medium-frequency range, the following three types of interesting couplings (highlighted by the high value of $||λ₁(St_k, St_l)||$) can be identified, which will be analysed in more details in the following paragraphs.

(i) The interactions in the difference interaction region correspond to a cascade of interactions generating low frequencies from medium frequencies. This cascade is indicated by the solid black arrow.

(ii) The interactions in the sum interaction region correspond to a cascade of sum interactions generating medium frequencies from low frequencies. This cascade is indicated by the dashed arrow.

(iii) The interactions circled correspond to interactions associated with the linear evolution of the flow around the mean flow that contribute significantly to the dynamics of the flow.

In the following, in order to further interpret the physical role of the triadic interactions that we have highlighted, for each type of interaction, we analyse in more detail the involved frequencies, as well as their characteristic interaction maps and the shape of the associated bispectral modes.

## 4.3. Cascade of interactions generating low frequencies from medium frequencies (type I)

In order to get more insight in the low- and medium-frequency range, the map of mode bispectrum amplitude is reproduced in figure 17 by zooming in these ranges.

In the region of the map below the horizontal axis (difference interaction region), there is a large continuous zone of interaction of significant amplitude defined by $(St_1, St_2) \in ]0, 0.18] \times ]0, -0.18]$. These interactions are difference interactions, each interaction therefore results in the creation of a lower frequency $St_3$ from two interacting frequencies $St_1$ and $St_2$. As this zone is continuous, a multitude of cascades of successive interactions exist, creating low frequencies from interactions in the mid-frequency range characteristic of the well-documented flapping of the shear layer $St \in [0.1, 0.18]$.

As an illustration, an example of such a cascade is shown in figure 18. Each interaction is given a number (between brackets) that is reported on the mode bispectrum amplitude map (figure 17). We would stress that the cascade considered here is just one example
<figure><img src="image_18.png" /><figcaption>Figure 16. Modulus of the complex mode bispectrum of the centred state vector $q^c$ in the $(St_1, St_2)$ plane. Relevant types of interactions are circled or indicated by arrows: interactions between medium frequencies creating low frequencies (arrow); interactions between low frequencies creating medium frequencies (dashed arrow); interactions reflecting the linear evolution of the flow around its mean. Interactions between unresolved frequencies are greyed.</figcaption></figure>

of a multitude of other possible paths. The first column (interaction of numbers (1), (2), (3) and (4)) contains interactions between medium frequencies present in the flow (typical frequencies of the well-documented flapping mode) that are written in red. The presence of these medium frequencies in the flow is sufficient to initiate a cascade of interactions that eventually creates low frequencies. Indeed the interactions of the first column generate lower frequencies. In column two, we report interactions of these generated lower frequencies with medium frequencies (interaction of numbers (5), (6) and (7)). These interactions generate low frequencies typical of the breathing mode, written in orange. We show in the third column that interactions between frequencies generated from interactions of the first and second columns also contribute to the generation of low frequencies characteristic of the breathing (interaction of numbers (8), (9), (10) and (11)).

This example of cascade shows how the presence in the flow of the well-documented medium frequencies characteristic of the flapping generates low frequencies through a multitude of cascades of difference interactions.

The interaction map of an emblematic triadic interaction of type I is shown in figure 19 for the three components of velocity. Other interactions of type I are not shown here for the sake of brevity, but they exhibit the same qualitative behaviour. We clearly identify that
<figure><img src="image_19.png" /><figcaption>Figure 17. Modulus of the complex mode bispectrum of the centred state vector $q^c$ in the ($St_1$, $St_2$) plane: zoom on medium/low frequencies. Sequences of interactions of each type are highlighted and numbered. Interactions between unresolved frequencies are greyed.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 18. Example of triadic interactions sequence between medium frequencies eventually creating low frequencies. The frequencies are expressed as $St_{int}$. Medium frequencies at the origin of the sequence are written in red. The resulting low frequencies are written in orange. Each interaction has a number between parentheses reported in figure 17.</figcaption></figure>

the nonlinear interactions constituting the cascade creating energy at low frequencies from medium frequencies are mainly localised in the second part of the separation bubble and especially in the region of reattachment. This point is consistent with the analyses and interpretations exposed in previous sections, namely the direct analysis of the DNS database (§ 2.5) and the SPOD analysis (§ 3), in which the low-frequency breathing dynamics of the separation bubble has been shown to be driven by dynamical activity taking place in the downstream part of the separation bubble. Moreover, this result confirms and refines the results of Sansica *et al.* (2014) and Mauriello *et al.* (2022). Indeed, both works suspected a role of nonlinear interactions in generating low frequencies in the downstream part of the
<figure><img src="image_21.png" /><figcaption>Figure 19. Spanwise averaged interaction maps for a triadic interaction (interaction (7) in figures 17 and 18) illustrating the sequence of interactions creating low frequencies from medium frequencies (interactions of type I. Frequencies involved: $St_1 = 0.12, St_2 = -0.08, St_3 = 0.04$. The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$. (b) Spanwise velocity $v$. (c) Vertical velocity $w$.</figcaption></figure>

interaction zone, although the precise determination of the location of significant triadic interactions was difficult due to the local nature of the methods employed.

The low-frequency bispectral mode resulting from the triadic interactions of type I evoked in figure 19 is shown in figure 20 for the three components of velocity. For all components of velocity, we clearly see that the bispectral mode is very similar to the first SPOD mode at the same frequency ($St = 0.04$, characteristic of the breathing) shown in figure 13. This result strongly suggests that the low-frequency breathing dynamics is fuelled by interactions between medium-frequency modes characteristic of flapping.
<figure><img src="image_22.png" /><figcaption>Figure 20. Spanwise averaged bispectral mode for a triadic interaction (interaction (7)) in figures 17 and 18 illustrating the sequence of interactions creating low frequencies from medium frequencies (interactions of type I. Frequencies involved: $St_1 = 0.12$, $St_2 = -0.08$, $St_3 = 0.04$. The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$. (b) Spanwise velocity $v$. (c) Vertical velocity $w$.</figcaption></figure>

4.4. Cascade of interactions creating medium-frequency modes from low frequencies
(type II)

By analysing figure 17, we identify a region of high amplitudes of the mode bispectrum above the horizontal axis (sum interactions), defined by ($St_1$, $St_2$) ∈ [0, 0.15]×[0, 0.1]. This region is symmetrical to the region of type I interactions analysed in the previous paragraph. It has an analog effect on flow dynamics. However, as it is populated by sum interactions, it consists of a multitude of cascades of interactions generating medium frequencies from interactions between low frequencies instead of the opposite.

In the same way as for interactions of type I, we show in figure 21 a particular path for a cascade of such interactions. Each interaction is given a number (between brackets)
<figure><img src="image_23.png" /><figcaption>Figure 21. Example of a cascade of triadic interactions creating medium frequencies from low frequencies. The frequencies are expressed as $St_{int}$. Low frequencies at the origin of the cascade are written in red. The resulting medium frequencies are written in orange. Each interaction has a number between parentheses reported in figure 17.</figcaption></figure>

that is reported on the mode bispectrum amplitude map (figure 17). To form this cascade, we only assume the existence of two low frequencies marked in red, i.e. here $St = 0.02$ and $St = 0.0267$. The first stage of the cascade consists of the self-interactions of these frequencies and the interactions of these frequencies with each other, creating higher frequencies. We then unfold a cascade of interactions between the frequencies created and show how medium frequencies (marked in orange) can be created in one or two stages. Again, as the region of high amplitude sum interactions is continuous, a multitude of cascades exist that has the same effect on the flow: creating medium frequencies, characteristic of the flapping, from low frequencies, characteristic of the breathing of the separation bubble.

The interaction map of an emblematic triadic interaction of type II is shown in figure 22 for the three components of velocity. Other interactions of type II are not shown here for the sake of brevity, but they exhibit the same qualitative behaviour. Again, these interaction maps clearly shows that the cascade of interactions creating medium-frequency modes from low frequencies mainly take place in the second part of the separation bubble and especially in the region of reattachment. This result is consistent with the SPOD analysis of the dominant dynamics at medium frequency (mode 1 at $St = 0.0934$) that shows a driving role of the reattachment region's dynamics in the flapping's dynamics. The direct analysis of the DNS database also revealed a strong intensification of the medium-frequency dynamics in the second part of the separation bubble.

The medium-frequency bispectral mode resulting from the triadic interactions of type II evoked in figure 22 is shown in figure 23 for the three components of velocity. The medium-frequency bispectral modes resulting from the cascade of interactions of type II are very similar to the medium-frequency flapping modes highlighted by the SPOD analysis and shown in Appendix A. This result suggests a feedback loop between the medium and low frequencies. Indeed, the analysis of the triadic interactions of type I (undertaken in § 4.3) suggest that low-frequency modes (breathing) are fuelled by triadic interactions between medium frequencies characteristic of flapping. In turn, the current analysis of interactions of type II suggest that the medium-frequency flapping dynamics is fuelled by a cascade of triadic interactions between low-frequency modes characteristic of the breathing of the separation bubble.

## 4.5. Modes arising from the linear evolution of the flow (type III)

We clearly see a range of high amplitudes along the horizontal axis of the map of the mode bispectrum amplitude. These points, circled in figure 16, correspond to interactions for
<figure><img src="image_24.png" /><figcaption>Figure 22. Spanwise averaged interaction map for a triadic interaction (interaction (18) in figures 17 and 21) illustrating the cascade of interactions creating medium-frequency modes from low frequencies (interactions of type II. Frequencies involved: $St_1 = 0.04667$, $St_2 = 0.04667$, $St_3 = 0.0933$. The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$. (b) Spanwise velocity $v$. (c) Vertical velocity $w$.</figcaption></figure>

which $St_1 \in [0.02, 0.18]$ and $St_2 = 0$. Hence, $St_3 = St_1$ and the corresponding bispectral modes can be interpreted as the linear evolution of the flow around its mean, according to Schmidt (2020). Interestingly, the range of frequencies of these modes corresponds to the range of low and medium frequencies characteristic of the SWBLI unsteadiness, including breathing dynamics and flapping dynamics. The flow therefore shows a strong linear tendency of the flow to oscillate around its mean field at medium and low frequencies characterising respectively the flapping and the breathing of the bubble. In particular, we observe high intensities at the frequencies identified by the direct analysis of the DNS data and the SPOD analysis (peak frequencies of the first mode), namely $St_1 = 0.04$,
<figure><img src="image_25.png" /><figcaption>Figure 23. Spanwise averaged bispectral mode for a triadic interaction (interaction (18) in figures 17 and 21) illustrating the cascade of interactions creating medium-frequency modes from low frequencies (interactions of type II. Frequencies involved: $St_1 = 0.04667$, $St_2 = 0.04667$, $St_3 = 0.0933$. The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$. (b) Spanwise velocity $v$. (c) Vertical velocity $w$.</figcaption></figure>

$St_1 = 0.0667$ and $St_1 = 0.0934$. These three modes are reported in figure 17 as numbers (27), (28) and (29).

We show below the bispectral modes for the modes oscillating at $St_1 = 0.04$ and $St_1 = 0.0934$ in figures 24 and 26, respectively. These modes are qualitatively very similar to the low-frequency breathing and medium-frequency flapping modes obtained by SPOD analysis. It shows that these linear modes are contributing to the breathing and flapping modes of the separation bubble.

The interaction maps corresponding to these modes are shown in figures 25 and 27. The amplitude of these maps are mainly significant in the second part of the separation bubble and especially in the reattachment region. It clearly shows that these medium- and
<figure><img src="image_26.png" /><figcaption>Figure 24. Spanwise averaged bispectral mode for a triadic interaction (interaction (28) in figure 17) illustrating the interactions expressing the linear evolution of the flow around its mean field (interactions of type III. Frequencies involved: $St_1 = 0.04$, $St_2 = 0$, $St_3 = 0.04$. The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$. (b) Spanwise velocity $v$. (c) Vertical velocity $w$.</figcaption></figure>

low-frequency modes highlighted are driven by the dynamical activity in the downstream part of the interaction, similarly to the modes resulting from interactions of type I and II.

## 5. Discussion

At this point, it is important to remember that the medium-frequency dynamics of separation bubbles is a well-documented phenomenon on which there is a consensus. The motivation of the preceding analyses was the documentation of the lower-frequency breathing-type dynamics that is associated with reflected shock oscillations in the case of
<figure><img src="image_27.png" /><figcaption>Figure 25. Spanwise averaged interaction map for a triadic interaction (interaction (28) in figure 17) illustrating the interactions expressing the linear evolution of the flow around its mean field (interactions of type III. Frequencies involved: $St_1 = 0.04$, $St_2 = 0$, $St_3 = 0.04$. The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$. (b) Longitudinal velocity $v$. (c) Longitudinal velocity $w$.</figcaption></figure>

SWBLIs, and whose mechanisms are still under debate. In this respect we are particularly interested in highlighting the nonlinear link between these two ranges of frequencies.

The analyses conducted in the preceding sections lead to the following factual conclusions.

(a) The dominant dynamics of this flow are in the low- and medium-frequency range, with breathing and flapping modes associated with reflected shock oscillations at these frequencies.

(b) The dynamics of these modes are qualitatively similar. The associated frequencies develop in the mixing layer to reach maximum amplitudes in the downstream zone of the separation bubble close to the reattachment.
<figure><img src="image_28.png" /><figcaption>Figure 26. Spanwise averaged bispectral mode for a triadic interaction (interaction (28) in figure 17) illustrating the interactions expressing the linear evolution of the flow around its mean field (interactions of type III. Frequencies involved: $St_1 = 0.0934$, $St_2 = 0$, $St_3 = 0.0934$. The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$. (b) Spanwise velocity $v$. (c) Vertical velocity $w$.</figcaption></figure>

(c) In addition, for both modes, there is an upstream propagation of information inside the bubble that forces the separation zone.

(d) The downstream zone of the separation bubble is the site of significant triadic interaction cascades providing energy to the breathing dynamics from interactions between existing frequencies of the shear layer flapping. In return, these low-frequency modes interact nonlinearly, also in the downstream part of the bubble, and form cascades of triadic interactions feeding energy into the flapping modes of the shear layer.

In light of these elements, we can propose a mechanism underlying the SWBLI unsteadiness. This mechanism is illustrated in figure 28.
<figure><img src="image_29.png" /><figcaption>Figure 27. Spanwise averaged interaction map for a triadic interaction (interaction (28) in figure 17) illustrating the interactions expressing the linear evolution of the flow around its mean field (interactions of type III. Frequencies involved: St₁ = 0.0934, St₂ = 0, St₃ = 0.0934. The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity *u*. (b) Longitudinal velocity *v*. (c) Longitudinal velocity *w*.</figcaption></figure>

For the description, the flow is divided into three zones: zone 1, before the interaction; zone 2, the interaction zone itself and zone 3, downstream of the interaction. The first zone consists of the boundary layer upstream of the interaction zone that amplifies free-stream perturbations. The flow then enters zone 2 by passing through the reflected shock wave. For given Mach and Reynolds numbers of the free-stream flow, the reflected shock strength does not depend on the incident shock wave strength as predicted by the free interaction theory. Through the shock, the disturbances are slightly amplified. At this stage, for sufficiently small free-stream upstream disturbances, it is reasonable to assume that the disturbances are still in the linear regime. These perturbations act as a forcing of the shear layer bounding the separation bubble in the separation region. These perturbations
<figure><img src="image_30.png" /><figcaption>Figure 28. Diagram explaining the suspected mechanism underlying the low-frequency unsteadiness of the transitional SWBLI.</figcaption></figure>

are then amplified in the shear layer in which convective instabilities also develop. If the interaction is strong enough, the separation bubble is large enough and the shear layer is long enough so that the linear shear layer modes saturate and exhibit a nonlinear dynamics with cascades of triadic interactions generating the low-frequency activity. In the present case, after the reattachment, the flow enters into zone 3 (downstream of reattachment), where the vortical structures created in the interaction zone are propagated and finally lead the boundary layer to transition to turbulence, which is not the subject of the present study. Back in zone 2, there is an upstream transport of a part of the dynamical activity at medium and low frequency from the reattachment zone until the separation zone. This upstream propagation has been shown to be damped and selective, as it damps high frequencies and selects low and medium frequencies. These upstream transported perturbations act as finite amplitude perturbations in the receptivity process of the shear layer, along with the upstream boundary layer perturbations. The shear layer is therefore forced by perturbations at low and medium frequencies characteristic of the breathing and flapping.

The mechanism thus described is a feedback loop whose key elements are the generation of low frequencies by the nonlinear dynamics of the shear layer and the upstream propagation of these disturbances within the mixing layer.

# 6. Conclusions and perspectives

## 6.1. *Conclusions*

The aim of this paper was to improve our physical understanding of the low-frequency dynamics of SWBLIs, and more specifically, the phenomenon of SWBLI unsteadiness. Pursuing this objective, we turned our attention to the case of strong interactions, creating large separation bubbles, for which the consensus is that the low-frequency dynamics of the interaction is an intrinsic feature of the recirculation bubble. The specific interest of our study was to document the possible role of triadic interactions in the mechanism behind low-frequency instationnarities.
A DNS of a strong impinging shock wave SWBLI at $M = 2.15$ on a flat plate have therefore been performed and analysed. In order to avoid the possible influence of the upstream boundary layer dynamics on the interaction dynamics, we have chosen to simulate an interaction with a laminar incident boundary layer, without forcing.

First, the DNS database was analysed to determine the average flow topology, its characteristic frequencies and the nature of the vortex structures within the interaction. This analysis showed that the flow exhibited a large, quasi-symmetric mean separation bubble with respect to the vertical axis passing through the apex of the bubble. Instantaneous field analysis showed that the vortical activity is mainly concentrated in the shear layer and in the downstream part of the separation bubble. Görtler-types structures have been identified inside the separation bubble downstream of its apex, whereas the downstream part of the shear layer (after the incident impingement) have been shown to be populated by longitudinal streaks. This dynamical activity of the late part of the interaction leads to transition after the reattachment point. The analysis of the frequency content across the interaction zone showed that the dominant dynamical activity of this flow is in the low- and medium-frequency ranges characterising respectively the breathing of the separation bubble and the flapping of the shear layer. The associated frequencies develop in the mixing layer to reach maximum amplitudes in the downstream zone of the separation bubble close to the reattachment. Part of this signal is emitted from the downstream part of the separation bubble and travels upstream until eventually forcing the shear layer in the separation region. This upstream propagation has also been shown to be damped and selective, as it selects low and medium frequencies.

To get more insight in the physical structure of the flow at these frequencies, a SPOD analysis has been performed. Interestingly, the peak frequencies of the dominant mode have been found to match the frequencies highlighted by the direct analysis of the DNS database. The corresponding spatial modes consist in separation bubble oscillations typical of the breathing and flapping modes accompanied by oscillations of the reflected shock. Moreover, the SPOD modes at these frequencies were found to be piloted by the downstream part of the interaction and the upstream propagation of perturbations originating from this zone inside the separation has been highlighted in the space-time evolution of these dominant modes eventually reaching the separation point.

A BSMD analysis has then been performed in order to characterise the triadic interactions at play in the interaction, especially in the downstream part of the shear layer, and to which such existent interactions could play a role in the onset of the breathing dynamics, given that the flapping dynamics of the shear layer is now well documented and is a consensus. This analysis has shown that the strong dynamical activity in the downstream part of the interaction is largely fuelled by strong triadic interactions in this zone. Indeed, this region is the seat of strong triadic interactions between frequencies characterising the flapping, contributing to the breathing activity. In turn, in the same region of the flow, triadic interactions between frequencies in the range characterising the breathing contribute to the flapping activity.

All these results strongly suggest a mechanism underlying the breathing motion of the separation bubble involving the nonlinear saturation of the convective instabilities of the shear layer, creating energy at low and medium frequencies in the downstream part of the separation bubble. These perturbations are then transported upstream inside the bubble and force the shear layer in the upstream part of the separation bubble. Such a feedback loop could participate in the onset of the low-frequency breathing of the bubble and, therefore, the low-frequency unsteadiness of the SWBLI.
I. Ben Hassan Saïdi, S. Wang, G. Fournier, C. Tenaud and J.-C. Robinet

## 6.2. Perspectives

Another important feature of the flow, highlighted by the BSMD analysis, is the existence of linear modes oscillating around the mean flow in the low- and medium-frequency range characterising breathing and flapping (type III interactions). This suggests other possible contributions to the low-frequency unsteadiness that were not further investigated as it was out of the scope of the present paper that was focused on the role of triadic interactions. Indeed, these results indicate the possibility of the existence of unstable or slightly stable global modes around the mean flow field in the low- and medium-frequency range. In the case of slightly stable modes, the perturbations arising from the nonlinear dynamics of the downstream part of the bubble could then act as a finite amplitude forcing of these modes at the appropriate frequency, thus sustaining their oscillations. Such a scenario is not inconsistent with the fact that globally unstable modes are not documented in the literature for stability analyses around the fixed point. Indeed, these analyses around the fixed point neglect to take into account at least part (the distortion) of the nonlinear effects linked to the convective dynamics responsible for the transition. A stability analysis around the mean field, on the other hand, takes partial account of the effect of nonlinearities, which play a major role in the described scenario. Such a stability study could be carried out in future work to confirm the presumptions based on the results of the present BSMD analysis. It should be noted, however, that such a neutral mode has been highlighted by a linear stability analysis around the mean flow of an axisymmetrical compression ramp in Lugrin et al. (2022). Moreover, in a recent work (Cura et al. 2024) the authors performed a global stability analysis around the mean flow of an incompressible separation bubble (with a non-fixed-separation point), whose results suggest that the breathing motion is driven by such a forced modal mechanism. This result lends further credence to our hypothesis that calls for a future work in which the global stability analysis of the flow around its mean field should be undertaken.

Finally, in this work we considered an idealised flow in which the upstream boundary layer is a laminar unforced boundary layer. This case is a limit case of more realistic flows in which the incoming boundary layer is forced by environmental disturbances and has been chosen to serve as a reference. As perspectives of this work, we will perform analog analyses for the interaction with a forced incident laminar boundary layer in order to characterise the effect of this forcing on the low-frequency dynamics of transitional SWBLIs. Future studies should also be undertaken involving turbulent incident boundary layers to study the fate of the mechanisms highlighted in transitional interactions when submitted to an intense turbulent forcing.

**Supplementary movies.** Supplementary movies are available at https://doi.org/10.1017/jfm.2025.236.

**Acknowledgements.** The authors thank the “Compressible Flows” team at IUSTI Marseille and more particularly L. Larchevèque and P. Dupont for their discussions around the dynamics of compressible separation.

This work was granted access to the HPC resources of IDRIS under the allocation 2024- A0152A07195 made by GENCI.

**Declaration of interests.** The authors report no conflict of interest.

## Appendix A. SPOD spatial modes

We reproduce here the first SPOD modes for each velocity components at $St = 0.0667$ (figure 29) and $St = 0.0934$ (figure 32) in order to support the assertion in § 3 that these modes are qualitatively identical to the first mode at $St = 0.04$. The similarity of the predominant dynamics at frequencies $St = 0.04$, $St = 0.0667$ and $St = 0.093$ can also be
<figure><img src="image_31.png" /><figcaption>Figure 29. First spatial SPOD mode ($m_1$) averaged in the spanwise direction at $S_t = 0.0667$. The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$. (b) Spanwise velocity $v$. (c) Vertical velocity $w$.</figcaption></figure>

seen in the time evolution of these modes (shown below for frequencies $S_t = 0.0667$ and $S_t = 0.0934$ in figures 30, 31, 33 and 34). In particular, the upstream propagation of information from the reattachment zone inside the bubble is highlighted in figures 31 and 34.

## Appendix B. Mesh convergence

Three meshes have been considered to assess the quality of the mesh employed in this paper:
<figure><img src="image_32.png" /><figcaption>Figure 30. Visualisation of $\Phi_1(x, t)$ velocity $u$ at $St = 0.0667$, averaged in the spanwise direction. (a) Velocity component $u$. (b) Velocity component $v$. (c) Velocity component $w$.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Figure 31. Representation of $\Phi_1(x, z = \text{cste}, t)$ at $St = 0.0667$ for $z/\delta = 1.48$. The vertical black lines indicate, from left to right, the separation point, the crossing of the rising shear layer, the crossing of the incident shock wave, the crossing of the descending shear layer and the reattachment point. (a) Velocity component $u$. (b) Velocity component $v$.</figcaption></figure>

a. mesh used in the present paper: 800 x 400 x 202 grid points in (x × y × z) and the minimum grid spacing over the plate is $\Delta z_{min} = 0.0125\delta$ in the normal-to-the-wall direction;

b. coarse mesh: 400 x 200 x 102 grid points in (x × y × z) with the same minimum grid spacing over the plate in the normal-to-the-wall direction;
<figure><img src="image_34.png" /><figcaption>Figure 32. First spatial SPOD mode ($m_1$) averaged in the spanwise direction at $S_r = 0.0934$. The mean flow is indicated by isolines of the mean density field. (a) Longitudinal velocity $u$. (b) Spanwise velocity $v$. (c) Vertical velocity $w$.</figcaption></figure>

c. finer mesh: 1600 × 800 × 202 grid points in ($x \times y \times z$) with the same minimum grid spacing over the plate in the normal-to-the-wall direction.

The grid spacing evolution along the wall is shown in figure 35, in wall units, for each coordinate direction for the mesh used in this study and the finer mesh. The two meshes share the same grid spacing at the wall in the normal-to-the-wall direction. The distribution of $\Delta z_+$ obtained with the two meshes is therefore similar throughout the whole flat plate and remains less than 1 everywhere (and even less than 0.5 in the interaction zone), indicating a sufficient resolution of the mesh in the normal-to-the-wall direction. For the mesh used in this paper, in the longitudinal and transverse directions, $\Delta x_+$ and $\Delta y_+$ remain at low values throughout the entire interaction region (less than 10 before the reattachment point $\bar{x}_R \simeq 88.12\delta$). After the reattachment, where the boundary layer becomes turbulent, the value of $\Delta x_+$ and $\Delta y_+$ increase until a larger
<figure><img src="image_35.png" /><figcaption>Figure 33. Visualisation of $\Phi_1(x, t)$ velocity $u$ at $St = 0.0934$, averaged in the spanwise direction. (a) Velocity component $u$. (b) Velocity component $v$. (c) Velocity component $w$.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Figure 34. Representation of $\Phi_1(x, z=\text{iste}, t)$ at $St = 0.0934$ for $z/\delta = 1.48$. The vertical black lines indicate, from left to right, the separation point, the crossing of the rising shear layer, the crossing of the incident shock wave, the crossing of the descending shear layer and the reattachment point. (a) Velocity component $u$. (b) Velocity component $v$.</figcaption></figure>

value, which remains however under 20. Notably, the value of 15 is reached around $x = 100\delta$, i.e. more than 10 boundary layer thicknesses downstream of the reattachment point. The finer mesh has been designed to lower the distribution of $\Delta x_+$ and $\Delta y_+$ and see if it influences significantly the accuracy of the result. To this end, the finer mesh
<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /><figcaption>Figure 35. Distribution along the flat plate of the grid spacing at the wall in each direction. Results are shown for (a) $\Delta y_+$, (b) $\Delta x_+$ and (c) $\Delta z_+$. Refined mesh (____) and mesh of the present paper (____).</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Figure 36. Time mean longitudinal velocity field averaged in the spanwise direction $\bar{u}$. The coloured flood corresponds to the mesh used in the paper whereas the isocontour green lines correspond to the finer mesh.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Figure 37. Distribution of the spanwise averaged wall pressure along the flat plate: refined mesh (____), mesh of the present paper (____), coarse mesh (____).</figcaption></figure>

has been obtained by doubling the number of points in the streamwise and spanwise directions.

The influence of the mesh refinement on the solution can be evaluated in figures 37 and 38 that show the distribution of the spanwise averaged wall pressure and the distribution of the spanwise averaged skin friction along the flat plate for the three considered meshes, respectively.

We clearly see in figure 37 that the pressure distribution at the wall is not affected by the mesh refinement.

In figure 38 we see that the coarser mesh in not fine enough to catch precisely the distribution of $C_f$ in the downstream part of the separation bubble. Moreover, we clearly see that the finer mesh solution and the paper's mesh solution agree completely until
<figure><img src="image_42.png" /><figcaption>Figure 38. Distribution of the spanwise averaged skin friction along the flat plate: refined mesh (___), mesh of the present paper (—), coarse mesh (---). Blasius laminar boundary layer solution (___). Turbulent solution (----).</figcaption></figure>

$x \simeq 110\delta$ (more than 20 boundary layer thicknesses downstream of the reattachment point). The two solutions show some small discrepancies downstream this abscissa but both follow the turbulent solution. The turbulent skin friction correlation plotted in this figure is taken from Cousteix (1989):

$$Cf = \frac{0.0172}{Re_{\theta}^{1/5}} f^{6/5}. \qquad (B1)$$

Here the compressibility function is

$$f = \left( \frac{\mu^*}{\mu_e} \right)^{1/6} \left( \frac{T_e}{T^*} \right)^{5/6}, \qquad (B2)$$

where subscript * refers to the free-stream value (outside the boundary layer) and the superscript * refers to the reference state given by

$$T^* = T_e + 0.7(T_w - T_e) \qquad (B3)$$

for an adiabatic boundary layer.

In order to have a comparison of the fields obtained in the whole domain with the mesh used in this paper and the finer mesh, we show in figure 35 the comparison of the time mean longitudinal velocity field averaged in the spanwise direction obtained with the mesh of the present paper and the finer mesh. It shows good agreement between the two obtained mean flow fields throughout the entire domain.

These comparisons between the results obtained with the mesh employed in this paper and the finer mesh show that a mesh refinement, multiplying the computation cost by 4 (we multiply the number of points by 2 in the streamwise and spanwise direction), has only an influence on the skin friction distribution starting from 20 boundary layer thicknesses downstream of the reattachment point, which is a part of the flow field not investigated in our study. Moreover, it shows that the mesh used in this paper already shows a fairly good agreement with the turbulent boundary layer solution after the transition. For these reasons, the mesh employed in this paper has been selected as a good compromise between accuracy and simulation cost.

REFERENCES

ADLER, M.C. & GAITONDE, D.V. 2018 Dynamic linear response of a shock/turbulent-boundary-layer interaction using constrained perturbations. *J. Fluid Mech.* **840**, 291–341.

AUBARD, G., GLOERFELT, X. & ROBINET, J.-C. 2013 Large-eddy simulation of broadband unsteadiness in a shock/boundary-layer interaction. *AIAA J.* **51** (10), 2395–2409.

BABINSKY, H. & HARVEY, J.K. 2011 *Shock Wave-Boundary-Layer Interactions*. Cambridge University Press.
BEN HASSAN SAÏDI, I., FOURNIER, G. & TENAUD, C. 2020 On the behaviour of high-order one-step monotonicity-preserving scheme for direct numerical simulation of shocked turbulent flows. *Intl J. Comput. Fluid Dyn*, **34**(9), 671–704.

BEN HASSAN SAÏDI, I. 2019 *Numerical Simulations 1077 of the Shock Wave-Boundary Layer Interactions*. *Theses*, Université Paris Saclay (COmUE), France.

BUFFIN-BÉLANGER, T., ROY, A.G. & DEMERS, S. 2013 9.6 turbulence in river flows. In *Treatise On Geomorphology* (ed. SHRODER, J.F.), pp. 69–86, Academic Press.

BUGEAT, B., ROBINET, J.-CH, CHASSAING, J.-C. & SAGAUT, P. 2022 Low-frequency resolvent analysis of the laminar oblique shock wave/boundary layer interaction. *J. Fluid Mech.* **942**, A43.

CHERRY, N.J., HILLIER, R. & LATOUR, M.E.M.P. 1984 Unsteady measurements in a separated and reattaching flow. *J. Fluid Mech.* **144**, 13–46.

CHONG, M.S., SORIA, J., PERRY, A.E., CHACIN, J., CANTWELL, B.J. & NA, Y. 1998 Turbulence structures of wall-bounded shear flows found using DNS data. *J. Fluid Mech.* **357**, 225–247.

CHU, B.T. 1965 On the energy transfer to small disturbances in fluid flow (Part I). *Acta Mechanica* **1** (3), 215–234.

CLEMENS, N.T. & NARAYANASWAMY, V. 2014 Low-frequency unsteadiness of shock wave/turbulent boundary layer interactions. *Annu. Rev. Fluid Mech.* **46** (1), 469–492.

COUSTEIX, J. 1989 *Aérodynamique Turbulence et Couche Limite*. Cépaduès éditions.

CURA, C., HANIFI, A., CAVALIERI, A.V.G. & WEISS, J. 2024 On the low-frequency dynamics of turbulent separation bubbles. *J. Fluid Mech.* **991**.

DARU, V. & TENAUD, C. 2004 High order one-step monotonicity-preserving schemes for unsteady compressible flow calculations. *J. Comput. Phys.* **193** (2), 563–594.

DEGREZ, G., BOCCADORO, C.H. & WENDT, J.F. 1987 The interaction of an oblique shock wave with a laminar boundary layer revisited. An experimental and numerical study. *J. Fluid Mech.* **177**, 247–263.

DIOP, M., PIPONNIAU, S. & DUPONT, P. 2019 High resolution Ida measurements in transitional oblique shock wave boundary layer interaction. *Exp. Fluids* **60** (4), 57.

DÉLERY, J. & DUSSAUGE, J.-P. 2009 Some physical aspects of shock wave/boundary layer interactions. *Shock Waves* **19** (6), 453–468.

DOLLING, D.S. 2001 Fifty years of shock-wave/boundary-layer interaction research: what next? *AIAA J.* **39** (8), 1517–1531.

DUSSAUGE, J.-P., DUPONT, P. & DEBIÈVE, J.-F. 2006 Unsteadiness in shock wave boundary layer interactions with separation. *Aerosp. Sci. Technol.* **10** (2), 85–91.

GAITONDE, D.V. 2015 Progress in shock wave/boundary layer interactions. *Prog. Aerosp. Sci.* **72**, 80–99.

GANAPATHISUBRAMANI, B., CLEMENS, N.T. & DOLLING, D.S. 2007 Effects of upstream boundary layer on the unsteadiness of shock-induced separation. *J. Fluid Mech.* **585**, 369–394.

GANAPATHISUBRAMANI, B., CLEMENS, N.T. & DOLLING, D.S. 2009 Low-frequency dynamics of shock-induced separation in a compression ramp interaction. *J. Fluid Mech.* **636**, 397–425.

GUIHO, F., ALIZARD, F. & ROBINET, J.-C. 2016 Instabilities in oblique shock wave/laminar boundary-layer interactions. *J. Fluid Mech.* **789**, 1–35.

HiOCFD Committee 2016 4th international workshop on higher order CFD methods.

KIRK, T.M. & YARUSEVYCH, S. 2017 Vortex shedding within laminar separation bubbles forming over an airfoil. *Exp. Fluids* **58** (5), 43.

KIYA, M. & SASAKI, K. 1985 Structure of large-scale vortices and unsteady reverse flow in the reattaching zone of a turbulent separation bubble. *J. Fluid Mech.* **154**, 463–491.

LUGRIN, M., BENEDDINE, S., GARNIER, E. & BUR, R. 2022 Multi-scale study of the transitional shock-wave boundary layer interaction in hypersonic flow. *Theor. Comput. Fluid Dyn.* **36** (2), 277–302.

MAHALINGESH, N., PIPONNIAU, S. & DUPONT, P. 2023 Effects of shock-induced separation on boundary layer transitional mechanisms. In *AIAA SCITECH 2023 Forum*. American Institute of Aeronautics and Astronautics.

MAURIELLO, M., LARCHEVEQUE, L. & DUPONT, P. 2022 Non-linearities in the low-frequency dynamics of transitional shock wave / boundary layer interactions. In *56th 3AF International Conference on Applied Aerodynamics*.

PIPONNIAU, S., DUSSAUGE, J.P., DEBIÈVE, J.F. & DUPONT, P. 2009 A simple model for low-frequency unsteadiness in shock-induced separation. *J. Fluid Mech.* **629**, 87–108.

PIROZZOLI, S. & GRASSO, F. 2006 Direct numerical simulation of impinging shock wave/turbulent boundary layer interaction at M = 2.25. *Phys. Fluids* **18** (6), 065113

PRIEBE, S., TU, J.H., ROWLEY, C.W. & MARTIN, M.P. 2016 Low-frequency dynamics in a shock-induced separated flow. *J. Fluid Mech.* **807**, 441–477.

ROBINET, J.-CH 2007 Bifurcations in shock-wave/laminar-boundary-layer interaction: global instability approach. *J. Fluid Mech.* **579**, 85–112.
I. Ben Hassan Saïdi, S. Wang, G. Fournier, C. Tenaud and J.-C. Robinet

SANSICA, A. 2015 Stability and unsteadiness of transitional shock-wave/boundary-layer interactions in supersonic flows. *PhD thesis*, University of Southampton, UK.

SANSICA, A., SANDHAM, N.D. & HU, Z. 2014 Forced response of a laminar shock-induced separation bubble. *Phys. Fluids* **26** (9), 093601.

SANSICA, A., SANDHAM, N.D. & HU, Z. 2016 Instability and low-frequency unsteadiness in a shock-induced laminar separation bubble. *J. Fluid Mech.* **798**, 5–26.

SARTOR, F., METTOT, C., BUR, R. & SIPP, D. 2015 Unsteadiness in transonic shock-wave/boundary-layer interactions: experimental investigation and global stability analysis. *J. Fluid Mech.* **781**, 550–577.

SCHMID, P.J. 2010 Dynamic mode decomposition of numerical and experimental data. *J. Fluid Mech.* **656**, 5–28.

SCHMIDT, O.T. 2020 Bispectral mode decomposition of nonlinear flows. *Nonlinear Dyn.* **102** (4), 2479–2501.

SCHMIDT, O.T. & COLONIUS, T. 2020 Guide to spectral proper orthogonal decomposition. *AIAA J.* **58** (3), 1023–1033.

SHINDE, V., MCNAMARA, J., GAITONDE, D., BARNES, C. & VISBAL, M. 2019 Transitional shock wave boundary layer interaction over a flexible panel. *J. Fluid. Struct.* **90**, 263–285.

SIROVICH, L. 1987 Turbulence and the dynamics of coherent structures. I. Coherent structures. *Q. Appl. Maths* **45** (3), 561–571.

SONG, Z. & HAO, J. 2023 Global instability of the interaction between an oblique shock and a laminar boundary layer. *Phys. Fluids* **35** (8), 084121.

SOUVEREIN, L.J., DUPONT, P., DEBIÈVE, J.-F., DUSSAUGE, J.-P., VAN OUDHEUSDEN, B.W. & SCARANO, F. 2010 Effect of interaction strength on unsteadiness in shock-wave-induced separations. *AIAA J.* **48** (7), 1480–1493.

THEOFILIS, V. 2003 Advances in global linear instability analysis of nonparallel and three-dimensional flows. *Prog. Aerosp. Sci.* **39** (4), 249–315.

THEOFILIS, V. 2011 Global linear instability. *Annu. Rev. Fluid Mech.* **43** (1), 319–352.

THOMPSON, K.W. 1987 Time dependent boundary conditions for hyperbolic systems. *J. Comput. Phys.* **68** (1), 1–24.

TOUBER, E. & SANDHAM, N.D. 2011 Low-order stochastic modelling of low-frequency motions in reflected shock-wave/boundary-layer interactions. *J. Fluid Mech.* **671**, 417–465.

TOWNE, A., SCHMIDT, O.T. & COLONIUS, T. 2018 Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis. *J. Fluid Mech.* **847**, 821–867.

VIRTANEN, P. et al. 2020 SciPy 1.0: fundamental algorithms for scientific computing in Python. *Nat. Meth.* **17** (3), 261–272.

WANG, Z.J. et al. 2013 High-order CFD methods: current status and perspective. *Intl J. Numer. Meth. Flow* **72** (8), 811–845.

WEISS, J., MOHAMMED-TAIFOUR, A. & SCHWAAB, Q. 2015 Unsteady behavior of a pressure-induced turbulent separation bubble. *AIAA J.* **53** (9), 2634–2645.

WELCH, P. 1967 The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms. *IEEE Trans. Audio Electroacoust.* **15** (2), 70–73.

WU, W., MENEVEAU, C. & MITTAL, R. 2020 Spatio-temporal dynamics of turbulent separation bubbles. *J. Fluid Mech.* **883**, A45.