# From low-frequency oscillations to Markovian bistable stall dynamics

Ivan Kharsansky Atallah, Luc Pastur, Romain Monchaux, Laurent Zimmer

► To cite this version:

Ivan Kharsansky Atallah, Luc Pastur, Romain Monchaux, Laurent Zimmer. From low-frequency oscillations to Markovian bistable stall dynamics. Physical Review Fluids, 2024, 9 (6), pp.063902. 10.1103/PhysRevFluids.9.063902. hal-04727344

HAL Id: hal-04727344

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

Submitted on 9 Oct 2024

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
# From Low-Frequency Oscillations to Markovian Bistable Stall Dynamics

Ivan Kharsansky Atallah,¹,² Luc Pastur,¹ Romain Monchaux,¹ and Laurent Zimmer²

¹Fluid Mechanics Department, ENSTA Paris, Institut Polytechnique de Paris, F-91120 Palaiseau, France

²EM2C Laboratory, CNRS, CentralSupelec, Université Paris-Saclay,
8-10 rue Joliot-Curie, 91190 Gif-sur-Yvette, France

(Dated: October 9, 2024)

We study the dynamics of a fixed-wing at stall in a wind tunnel by measuring the aerodynamic forces. We report experimental evidence of a critical Reynolds number from which low frequency oscillations in the force are replaced by random bistable dynamics. In this new regime, the flow explores each state intermittently with long residence times. This stochastic process can be modeled as a continuous Markov chain, and equivalently, it shows a super-exponential scaling for the mean residence times. Furthermore, the probability density function of the lift coefficient exhibits the characteristic heavy tail of extreme events. Extreme minima and maxima are at the origin of the transitions. We analyzed the evolution of these tails using extreme value theory to identify the bifurcation points of the associated dynamical system. The results are in good agreement with a user-defined threshold method, the advantage being the unambiguity in the computation.

## I. INTRODUCTION

Following the pioneering work of Prandtl [1], one calls stall regime any flow pattern that involves a massive detachment of the boundary layer from a wall. Stall has been studied in various geometries like backward-facing steps, cavities, and airfoils, among others; for a nice review on the subject, see e.g. [2]. Stall in airfoils appears as the angle of attack $α$ increases and can be classified as one or a combination of three main types [3] that take into account the behavior of the separation point and the appearance or not of a laminar separation bubble. The dynamics of the bubble and hence, the laminar to turbulent transition of the boundary layer, is believed to be at the origin of complex stall phenomena in airfoils, such as low-frequency oscillations [4–8] or hysteresis [6, 9–12]. These phenomena, in turn, may have a dramatic effect on profile efficiency (wind turbine, propulsion...) or be a source of mechanical damage. Low-frequency oscillations are global oscillations of the flow with a defined frequency much lower than the vortex shedding frequency [4], while hysteretic loops in the lift force usually occur when two flow states coexist over a finite range of $α$, resulting in bi-stability [9]. Parameters like velocity, geometry, free stream turbulence, and acoustic disturbances can either amplify, delay or suppress the appearance of low-frequency oscillations and hysteresis [4, 13]. It has been shown [6, 14] that large-scale modes dominate the dynamics at stall, involving saddle-node and Hopf bifurcations in the case of low-frequency oscillations [6]. Multi-stability has been encountered in several fluid mechanic systems ranging from small-scale flows [15–19] to large-scale geophysical flows [20–22] and in the transition to turbulence [23–25]. They all share common statistical features such as exponential statistics of the residence times, defined as the time-lapse in each state, and their mean values orders of magnitude larger than any hydrodynamic timescale of the flows [25–27]. All these characteristics make this process fall into the theory of Markov chains in which the transitions between states are considered rare events, namely events in which the probability of occurrence is low, and they are approximately exponentially distributed [28]. Rare events are at the origin of extreme events [29] as they result from the synergistic action between a low-probability event and the intrinsic dynamics of the system [30]. Extreme events are large deviations from normality exhibited by specific observables of a system. For a comprehensive review of extreme events in fluid dynamics systems, see [30]. Under this framework, extreme rare events can either be analyzed by using tools of residence time statistics or extreme value laws as both are equivalent [31–33]. In this article, we experimentally study the stall dynamics of a thin symmetric airfoil (NACA0012) for Reynolds numbers ranging from 50 000 to 110 000 (see definition in the next section). We confirm that low-frequency oscillations can actually be observed on the lowest range of the Reynolds numbers, as reported in the numerical and experimental studies of [5, 34, 35]. We discovered that beyond a critical value of the Reynolds number, the airfoil dynamics transits into random memory-less switches between the fully attached and the fully detached flow states, replacing the low-frequency oscillations. In this new regime, the expected hysteresis loop is not observed due to the intermittent dynamics, which connects the two branches of solutions for a fixed angle of attack. To the best of our knowledge, this intermittent bi-stable dynamics has never been reported in airfoil stall dynamics yet and we show that it can be described by the theory of continuous Markov chains and extreme rare events.
<figure><img src="image_2.png" /><figcaption>FIG. 1: Experimental setup showing the wing of chord c and span s, the rotation plate, which controls the angle of attack α, on which the force balance is mounted, and the hot wire placed in the wake.</figcaption></figure>

## II. EXPERIMENTAL SETUP

All the experiments were carried out in a closed loop wind tunnel with a 450 x 450 mm² test section. The turbulence intensity is constant for all the velocities tested and equal to 0.4%. The model analyzed is a 3D printed NACA0012 wing with a chord length c = 120 mm and a span s = 450 mm (FIG. 1). The wing is flushed on one side to an endplate which in turn is connected to a balance-rotating platform (FIG. 1). This system allows to rotate the model with a resolution of 0.02° and to measure lift (L) and drag forces with a resolution of 1 mN. We also carried out hot-wire (HW) measurements to differentiate between purely aerodynamic effects and airfoil vibrations in the spectra of the forces. The probe model was P15 used together with a Mini CTA 54T30 from DANTEC, see FIG. 1 for its position in the flow. The temperature in the test section was measured via a thermocouple with an accuracy of 0.1 °C. Six Reynolds numbers were analyzed $Re = U_\infty c/\nu \in [5.31; 6.54; 7.84; 9.03; 10.3; 11.60] \times 10^4$, where $U_\infty$ is the velocity at the inlet of the test section and $\nu$ the air viscosity at the flow temperature. In order to accurately characterize phenomena around stall, the same rigorous measuring protocol was used for each of them. First, the stall zone is identified by doing quick measurements with large $\Delta\alpha$ steps; once the region was identified, small $\Delta\alpha$ steps and long measurements were performed. For the bistable case, the values adopted for the acquisition time and $\Delta\alpha$ were based on previous explorations that revealed the sensitivity of this phenomenon to small changes in the angle of attack and the long characteristic times associated. A moving average filter of 1 s was applied in order to highlight the bistability. For all measurements, the sampling frequency was set at 1 kHz. Most of the analysis shown in this article was done over the lift coefficient $C_L = 2L/(\rho U_\infty^2 c s)$, where L is the lift force, although drag and hot wire data were also analyzed and highlighted too the phenomena around stall.

## III. DYNAMICS AND BIFURCATIONS

### A. Dynamic regimes

As already stated, two different regimes were observed around stall when increasing $Re$. In order to analyze their evolution FIG. 2 was constructed by finding the range of existence of each phenomenon in the $\alpha-Re$ plane. The phase diagram is divided into two main regions: (i) below $Re_{bs} \approx 9 \times 10^4$, the flow exhibits a low frequency regime (LFO), for a certain range of angle of attack which depends on $Re$, and (ii) above $Re_{bs}$, the LFO regime disappears and a bi-stable regime appears over a range of angles of attack which also depends on $Re$. Below $Re_{bs}$, LFOs are observed for angles of attack between the two red-circles dashed lines, while above $Re_{bs}$, bistability is observed for angles of attack between the two magenta squares-dashed-lines. Stall angles defined as the values from which the lift decreases are also added as a reference. Low-frequency oscillations are observed when increasing $Re$ above a certain threshold $Re_{lfo} \approx 3 \times 10^4$ (not shown in the figure). The range of existence in $\alpha$ diminishes as $Re$ increases until the critical Reynolds number $Re_{bs}$ is reached, from where the phenomenon is not visible anymore. The characteristic signature
of LFO is found by looking at the power spectral density (PSD) of any of the measured variables. In FIG. 3(b), the hot wire spectrum is shown along with the corresponding time series FIG. 3(a), for a typical *Re* in this regime. A clearly defined peak is visible around $St = 0.007$ where $St = fc\sin(\alpha)/U_\infty$ is the Strouhal number based on the projected cross-section of the wing and *f* the frequency in Hz. This low-frequency large band peak is characteristic of wings showing thin airfoil stall [4, 5, 8]. Increasing *Re* in this zone only shifts the peak to a slightly larger value of the Strouhal number, up to $St \approx 0.01$ at $Re_{bs}$. This low-frequency peak is far from the wing-free vibration frequency. The latter can be computed with the following equation [36], $f_1 = (1.875^2/2\pi)\sqrt{EI/\rhoBs^4}$ where *E* is the Young's Modulus, *I* the second order area moment, $\rho$ the density, *B* the cross-section and *s* the span of the wing. With the dimensions and material of the airfoil, it yields $f_1 \approx 15$ Hz and $St \approx 0.045$, which is about an order of magnitude larger than that of the LFO. Beyond $Re_{bs}$, low-frequency oscillations are replaced by intermittent bistable dynamics distinguished by two different states that are intermittently explored in time. Only at $Re_{bs}$ both regimes could be observed. For a given $Re$ in the bi-stable regime, hot wire time series (FIG. 4(a)) shows a random switching between two states and the spectra (FIG. 4(b)) exhibits the characteristic $1/f^2$ behavior at low frequencies [15] with no sign of the low-frequency peak. Another remarkable feature of each phenomenon is the small range of their existence in $\alpha$, going as small as $\Delta\alpha = 0.18^\circ$. As the focus of this article is on the intermittent dynamics, from now on, all the results shown will be for a typical $Re$ in this regime, particularly $Re = 1.16 \times 10^5$. Similar results have been obtained for the other two $Re$ tested in the bistable zone.

We also analyzed the influence of vibrations, surface roughness, boundary conditions, and installation effects on the bistable regime. When both ends of the wing are fixed, the bistability threshold is shifted towards larger values of the Reynolds numbers, but the dynamics remain qualitatively the same. We interpret the delay of the transition as a consequence of the attenuation of the vibrations and, consequently, a reduction in the available energy to promote random jumps. Wind tunnel characteristic frequencies were also ruled out as a possible source of excitation for bistability. As stated in section II, the wind tunnel is a closed loop so for a fan with *n* blades, the blade pass frequency is $f_b = nr$ where *r* is the rotating speed of the fan wheel. In our experiment $f_b = 100$ Hz which is way larger than the timescales associated with the bistability. Preliminary tests done over an aluminum NACA 0012 profile with a smooth surface showed the intermittent bistable behavior around stall for a larger $Re_{bs}$, showing that the phenomenon is robust against surface roughness. This is similar to the influence of rugosity in the drag crisis on a sphere, which is universal in some aspects, but still strongly affected by some effects, notably rugosity [37, 38]. Thicker rugosity of 3D printed wing promotes turbulence and hence an earlier transition to bistability, whereas the smooth surface of an aluminum wing delays the transition to larger values of the Reynolds number. Regarding simulations, intermittent bursts in aerodynamic coefficients were observed in a NACA 4412 far from stall [39]. All of these results show the universality and robustness of this phenomenon, but also the non-universality of the transition threshold.

## B. Random switching dynamics

The mean behavior of the lift coefficient versus $\alpha$ is shown in FIG. 5. Samples inside the bistable region were acquired for 30 minutes and in steps of $\Delta\alpha = 0.04^\circ$ while outside this region, the time was 2 minutes and $\Delta\alpha = 1^\circ$. The mean lift ($C_{L_m}$) evolution is typical of thin airfoils [9, 40], showing a small change of slope at small angles due to the formation of the recirculation bubble, a soft decrease during stall followed by an increase at high angles. The inset in FIG. 5 follows closely the progression of the bistability, revealing two bifurcation points at $\alpha = \alpha_A \approx 10.90^\circ$ and $\alpha_D \approx 10.72^\circ$, between which two stable branches of high and low lift, respectively associated with an attached (A) or detached (D) flow regime, both coexist. The values corresponding to each branch i.e. $C_{L_A}$, $C_{L_D}$ are the two most likely values of $C_L$ in the bimodal probability mass functions (PMF). The bifurcation angles are determined by looking at the ratio between the probability associated with $C_{L_A}$, $C_{L_D}$. When this ratio goes below $10^{-3}$, we consider the dynamics mono-stable. It is important to note that no hysteresis of the mean lift nor of the onset of bistability $\alpha_{bs}$ was observed in decreasing and increasing ramps of the angle of attack. An interesting feature is that once bifurcated, each state remains almost the same throughout the bifurcation with only a small shift in value, and the overall average $C_{L_m}$ depends mostly on the time spent in each state. This fact is well understood while seeing FIG. 6 in where the PMF and the time series for three different angles are shown. At lower $\alpha$ (FIG. 6(a)) one state, A, prevails even though some attempts to explore D are visible. When $\alpha$ increases, D occurs more frequently (FIG. 6(b)), until a mono-stable state is reached when D prevails over A (FIG. 6(c)). Another salient feature is the standard deviation associated with each state. As it is clearly seen in the time series of FIG. 5, the fluctuations of the lift force in state D are much larger than those in state A, which is reflected in the width of the PMF peaks. High lift and small fluctuations (state A) are related to a mostly attached boundary layer state while low lift and large fluctuations (state D) are associated with a mostly separated boundary layer, as previously reported in [6, 41].
<figure><img src="image_3.png" /><figcaption>FIG. 2: Phase diagram $(\alpha, Re)$ showing the two observed regimes. Green shaded area corresponds to the range of existence of the low frequency regime (LFO) while the cyan shaded area corresponds to bistability. The abrupt transition between the LFO regime and bistability is marked by the vertical line at $Re = Re_{bs}$. The blue triangles mark the start of the stall zone, the red circles the extreme angles between which the LFOs are observed, and the magenta squares the extreme angles of the bistable dynamics.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3: Hot wire (a) time signal and (b) power spectral density (PSD), in the low frequency oscillation regime at $Re = 5.31 \times 10^4$. The broadband spectrum of oscillations is centred on $St \approx 6 \times 10^{-2}$.</figcaption></figure>

IV. STATISTICS OF THE RANDOM SWITCHING DYNAMICS

In this section, the statistics of the bistability are analyzed with the aim of understanding the nature of the stochastic
process and its evolution with the control parameters. Measurements with a duration of 1 hour were performed for
several angles in the region of interest as they showed fairly converged statistics and a reasonable overall time for a
given set. Nevertheless, isolated 8-hour acquisitions were also done to verify the statistical convergence.
<figure><img src="image_5.png" /><figcaption>FIG. 4: Hot wire (a) time signal and (b) power spectral density (PSD), in the bi-stable regime at $Re = 1.16 \times 10^5$. The characteristic slope in $St^{-2}$ is clearly visible on the Strouhal range $St \le 10^{-3}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 5: $C_L$ vs $\alpha$. Inset: zoom around stall exhibiting the bistability and its evolution. In the bistable range, the blue circles mark the high $C_L$ associated with flow state A, the red squares the low $C_L$ associated with flow state D. The black circles represent, for each value of $\alpha$, the time-averaged value of $C_L$.</figcaption></figure>

A. Continuous Markov chains and superexponential scaling

After a temporal filtering of 1 s, the signal is binarized to extract the residence times $\tau_i$, $i \in [A, D]$ of each state. The complementary cumulative distribution function (CCDF) of the residence times exhibits an exponential trend, $P(\tau_i, \alpha) \propto \exp[-\tau_i/\bar{\tau}_i(\alpha)]$, where $\bar{\tau}_i(\alpha)$ is the mean residence time, as shown in Fig. 7 for several angles of attack. The exponential distribution is typical of a memoryless random process, which can thus be described as a continuous Markov chain. The state space of the Markov process contains the two states $S = \{A, D\}$. In addition, as in bistable systems cited in the introduction, the mean residence times in each state are much larger than any typical time scales of the system. For example, compared to the convective time $c/U_\infty$ the residence time is 4 orders of magnitude higher. The transitions between states can hence be considered as rare events, e.g., the probability of seeing a transition is low.
<figure><img src="image_7.png" /><figcaption>FIG. 6: Time series (left) and probability mass functions (PMF) (right) of $C_L$ for (a) $\alpha = 10.70^{\circ}$, (b) $\alpha = 10.82^{\circ}$ and (c) $\alpha = 10.94^{\circ}$. The PMF are stitched on a high-lift state $C_{L_A}$ (a) or a low-lift state $C_{L_D}$ (c), or both (b).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIG. 7: Complementary cumulative distribution function (CCDF) of the residence times $\tau_A$ in state A (top) and $\tau_D$ in state D (bottom) for increasing or decreasing angle of attacks $\alpha$ (as marked by the arrow on the right figures). Dashed lines represent the exponential distribution fitted on the data.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 8: Mean residence time $\bar{\tau}$ evolution with the angle of attack $\alpha$ for states A (blue circles) and D (red squares). Solid lines are the double-exponential fits of the associated residence time distributions.</figcaption></figure>

The evolution of the mean residence time $\bar{\tau}_i(\alpha)$ (inverse of slopes in Fig. 7) with our control parameter $\alpha$ is presented in Fig. 8 in a double log-linear scale. The linear behavior is compatible with a double-exponential functional shape similar to other sub-critical systems like in the transition to turbulence [23, 42, 43] and the flow around a pendulum [16]. The crossing between lines has a particular interpretation which will be discussed in the next section. Note that because the range of angles of incidence is small (not even a decade), any super-exponential law could be a good fit for the data, as was already mentioned in Hof et al. [24] on the transition to turbulence. This super-exponential scaling directly results from the previously stated equivalence between continuous Markov chains and extreme value theory.

## B. Extreme rare events analysis

As already established, the transitions between states can be considered as rare events. Additionally, we showed the rise of an extreme value law as a consequence of the memory-less property of the process. All these facts lead us to consider extreme events as responsible for the transitions. To support this theory, the probability density function (PDF) of $C_L$ was plotted for two different angles of attack (see Fig. 9): one at the beginning of the bi-stable range (Fig. 9(a)) and another one at the end (Fig. 9(b)). The comparison is made against a Gaussian fit of the core of the PDF. The PDFs show the characteristic heavy tail typical of systems showing extreme events in both cases. We can also look at the amplitudes of these events in terms of the standard deviation of the Gaussian fit. They span from 10 to 20 times this value, showing the large deviation from normality expected from these types of events. We note that in Fig. 9(a) the extreme events are smaller values of $C_L$ e.g. minima, and in Fig. 9(b) larger values e.g. maxima. In the attached state A, the system undergoes extreme minima in $C_L$, triggering the transitions to D. In state D, extreme maxima bring back the flow to state A. We can now see the interest in knowing more about the distribution of maxima, $C_{Lmax}$ and minima $C_{Lmin}$. Extreme value theory states that in the asymptotic limit i.e when the sample size $n \to \infty$, these distributions converge to the generalized extreme value distribution (GEV) [32] with the following cumulative density function:

$$P = \exp\left\{-\left[1 + \xi\left(\frac{x - \mu}{\sigma}\right)\right]^{-1/\xi}\right\} \quad (1)$$

where $\mu$ is the location parameter, $\sigma$ is the scale parameter, and $\xi$ is the shape parameter. The latter is related to the kind of tail from the parent distribution, and consequently, it will determine the type of extreme value law. We call parent distribution, the probability distribution of $C_L$. When $\xi < 0$, the extremes follow a Weibull distribution, and the tail of the parent distribution is bounded. When $\xi > 0$, the extreme events follow a Frechet distribution. The tail
<figure><img src="image_10.png" /><figcaption>FIG. 9: Probability density function of $C_L$ showing the characteristic heavy tails of extreme events, when flow state (a) A or (b) D dominates in the flow dynamics. The dashed line is a Gaussian fit of the dominant peak.</figcaption></figure>

corresponds to a heavy one with a power decay law, as in Fig. 8. When $\xi = 0$, the extreme value law corresponds to a Gumbel distribution and an exponential decay tail from the parent distribution. Therefore, by following the evolution of the sign of $\xi$ with $\alpha$ we have insight into the behavior of extreme events in the bistable regime and hence into the probability of transition. This statistical procedure was applied in [32, 44] as a tool to detect the critical value of a control parameter in turbulent flows. Here, "critical value" is understood as the one producing a qualitative change in the dynamics, e.g. a bifurcation. This threshold corresponds to a change of sign in $\xi$, meaning the transition from a Weibull to a Frechet distribution or vice-versa, depending on the type of extreme events analyzed (maximum or minimum). We apply this method to our time series following the steps and recommendations detailed in [32, 44] with special care in the convergence of the GEV fit. All results are summarized in Fig. 10 along with representative PDFs of the parent distribution $C_L$ to help the understanding. For the maxima (blue circles), negative values of $\xi$ are observed until $\alpha \approx 10.78^\circ$, then a zero crossing between $\alpha \approx 10.78^\circ$ and $\alpha \approx 10.82^\circ$ with positive values up to $\alpha \approx 10.94^\circ$ where $\xi$ becomes negative again. The first zero-crossing can be interpreted by looking at Fig. 8. The angle of equal mean residence time for states A and D is between $\alpha \approx 10.78^\circ$ and $\alpha \approx 10.82^\circ$; before that angle, the system spends more time in A, and afterward, in D. This reversal of states coincides with the first change of sign in Fig. 10. By looking at PDFs inset, we can see the aforementioned reversal of probabilities. The second zero-crossing corresponds to the bifurcation angle $\alpha_A$ determined in section IIIB by means of a peak ratio criteria. Here, the PDF right before the transition (see Fig. 10(IV)) shows the characteristic heavy tail (Frechet), and after this angle, this tail disappears and is replaced by a bounded one (Weibull) like Fig. 10(V). Doing the same analysis for the minima (red squares) we see that $\xi$ goes from negative to positive around $\alpha \approx 10.70^\circ$ this value is close to the one previously defined for $\alpha_D$ i.e. the beginning of the bistable range. The system goes from bounded fluctuations (Fig. 10(I)) to sporadic bursts of extreme minima (Fig. 10(II)). Then $\xi$ continues to be positive until it crosses zero again between $\alpha \approx 10.78^\circ$ and $\alpha \approx 10.82^\circ$ corresponding to the reversal of states. We can point out that $\alpha_A$ is only identified when looking at the maximum and $\alpha_D$ at the minimum. This is because, in the first case, we look at the right tail of the distribution (state A) and, in the second case, the left tail of the distribution (state D). It is important to note that GEV are usually applied to distributions that are unimodal. Here we apply the GEV in an original way to distributions showing bimodality, i.e. in the bistable range. This use is subsequently validated by the smooth evolution of the shape parameter over the whole range of $\alpha$ for the maximum and minimum distributions. The GEV, when fitted to our bimodal data, captures well the dominant peak and the tail of interest (maxima or minima). Had we split the signal into two states, the result would have always been a Gumbel distribution ($\xi = 0$), missing the reversal of states.

V. DISCUSSION AND CONCLUSIONS

The experiments revealed an unexpectedly rich dynamics with two very different regimes (Fig. 2): at low Reynolds numbers, a deterministic oscillatory regime with a well-defined frequency, while at large Reynolds numbers, the
<figure><img src="image_11.png" /><figcaption>FIG. 10: Shape parameter $\xi$ vs $\alpha$ for the maxima (blue circles) and for the minima (red squares). Insets in black showing lift coefficient PDF (parent distribution) evolution with $\alpha$. The insets show the evolution of the pdf of the lift coefficient when the angle of attack increases, from the single flow state with high lift (I) to the single flow state with low lift (V), via the bistable regime (III) and intermediate regimes (II) and (IV).</figcaption></figure>

dynamics is stochastic with no sign of oscillations. This represents a transition from a mono-stable system with a limit cycle to a bistable system with random transitions between states. The experimental results and the study of [6] suggest that the two deterministic states A and D appear and disappear through two saddle-node bifurcations, one direct for $\alpha = \alpha_D(Re)$, the second indirect for $\alpha = \alpha_A(Re)$. This indicates that an unstable flow regime, presumably transiently observed during the transitions between states A and D, also co-exists with the A and D branches, which should collide with the A branch in $\alpha_A$ and with the D branch in $\alpha_D$, as detailed in [6]. In contrast to [6, 41] where A and D were absorbing states, the novel result here is that both states are explored intermittently over time, all other parameters being fixed. The hysteresis loop usually observed at stall could not be observed due to the intermittent switches between the two states, favored by the turbulent background fluctuations in the flow, analogously to the effect of thermal noise or inlet boundary conditions in other bistable systems [45, 46]. By analyzing the bistability statistics, we have shown that the process can be modeled as a continuous Markov chain. Former studies have already made analogies between such continuous Markov chains and extreme value laws in order to compute transition thresholds for the order parameter [16]. Here, we highlighted the mathematical link between continuous Markov chains and extreme value laws expressed by the super-exponential behavior of the mean residence times. Based on the hypothesis that extreme rare events are responsible for the transitions, we modeled the tail of the observable using extreme value theory. We used this model to compute the bifurcation points of our system in an unambiguous manner without the need to define thresholds by hand. Following [47] and although sometime found in literature, no critical dynamical value of the order parameter is sought, beyond which the system would irremediably change of state, since these bi-stable systems are likely to transit from one state to the other regardless of their current state. The current state of the system only affects the probability it has to transit from one state to the other.

Although outside of the scope of this article, a stochastically forced model of the double saddle-node bifurcation, together with a particle image velocimetry-based investigation of the flow structures are currently under investigation.
ACKNOWLEDGMENTS

We would like to thank T. Pichon, L. Cherfa, and E. Jean-Bart for their contribution to setting up the experiment.
We also thank Davide Faranda for his constant help in the analysis of extreme values theory

[1] O. Tietjens, Fundamentals of hydro-and aeromechanics, based on lectures by l. prandtl, Bull. Amer. Math. Soc **41**, 173 (1935).

[2] S. Kline, On the nature of stall, J Basic Eng-t Asme **81**, 305 (1959).

[3] G. B. Mccullough and D. E. Gault, *Examples of three representative types of airfoil-section stall at low speed*, Tech. Rep. (1951).

[4] K. B. M. Q. Zaman, D. J. Mckinzie, and C. L. Rumsey, A natural low-frequency oscillation of the flow over an airfoil near stalling conditions, J. Fluid Mech. **202**, 403–442 (1989).

[5] H. Tanaka, Flow visualization and piv measurements of laminar separation bubble oscillating at low frequency on an airfoil near stall, in *24th international congress of the aeronautical sciences* (2004) pp. 1–15.

[6] D. Busquet, *Study of a high Reynolds number flow around a two dimensional airfoil at stall; an approach coupling a RANS framework and bifurcation theory*, Ph.D. thesis, Institut polytechnique de Paris (2020).

[7] A. Broeren and M. Bragg, Low-frequency flowfield unsteadiness during airfoil stall and the influence of stall type, in *16th AIAA Applied Aerodynamics Conference*.

[8] P. J. Ansell and M. Bragg, Characterization of ice-induced low-frequency flowfield oscillations and their effect on airfoil performance, in *31st AIAA Applied Aerodynamics Conference*.

[9] T. Mueller, Unsteady stalling characteristics of thin airfoils at low reynolds number, in *Fixed and Flapping Wing Aerodynamics for Micro Air Vehicle Applications* (American Institute of Aeronautics and Astronautics, 2001) pp. 191–213.

[10] Z. Yang, H. Igarashi, M. Martin, and H. Hu, An experimental investigation on aerodynamic hysteresis of a low-reynolds number airfoil, in *46th AIAA Aerospace Sciences Meeting and Exhibit*.

[11] L. Pohlen and T. Mueller, Boundary layer characteristics of the miley airfoil at low reynolds numbers, in *Applied Aerodynamics Conference*.

[12] S. Mittal and P. Saxena, Prediction of hysteresis associated with the static stall of an airfoil, AIAA J. **38**, 933 (2000).

[13] J. A. Hoffmann, Effects of freestream turbulence on the performance characteristics of an airfoil, AIAA J. **29**, 1353 (1991).

[14] W. He, R. S. Gioria, J. M. Pérez, and V. Theofilis, Linear instability of low reynolds number massively separated flow around three naca airfoils, J. Fluid Mech. **811**, 701–741 (2017).

[15] M. Grandemange, M. Gohlke, and O. Cadot, Turbulent wake past a three-dimensional blunt body. part 1. global modes and bi-stability, J. Fluid Mech. **722**, 51–84 (2013).

[16] A. Gayout, M. Bourgoin, and N. Plihon, Rare event-triggered transitions in aerodynamic bifurcation, Phys. Rev. Lett. **126**, 104501 (2021).

[17] F. Ravelet, L. Marié, A. Chiffaudel, and F. Daviaud, Multistability and memory effect in a highly turbulent flow: Experimental evidence for a global bifurcation, Physical review letters **93**, 164501 (2004).

[18] R. Volpe, P. Devinant, and A. Kourta, Experimental characterization of the unsteady natural wake of the full-scale square back ahmed body: flow bi-stability and spectral analysis, Exp Fluids **56**, 99 (2015).

[19] G. Bonnavion and O. Cadot, Boat-tail effects on the global wake dynamics of a flat-backed body with rectangular section, J Fluid Struct **89**, 61 (2019).

[20] K. R. Sreenivasan, A. Bershadskii, and J. J. Niemela, Mean wind and its reversal in thermal convection, Phys. Rev. E **65**, 056306 (2002).

[21] Y. Tian, E. R. Weeks, K. Ide, J. Urbach, C. N. Baroud, M. Ghil, and H. L. Swinney, Experimental and numerical studies of an eastward jet over topography, Journal of Fluid Mechanics **438**, 129 (2001).

[22] M. Berhanu, R. Monchaux, S. Fauve, N. Mordant, F. Pétrélis, A. Chiffaudel, F. Daviaud, B. Dubrulle, L. Marié, F. Ravelet, M. Bourgoin, P. Odier, J.-F. Pinton, and R. Volk, Magnetic field reversals in an experimental turbulent dynamo, Europhys. Lett. **77**, 59001 (2007).

[23] N. Goldenfeld, N. Guttenberg, and G. Gioia, Extreme fluctuations and the finite lifetime of the turbulent state, Phys. Rev. E **81**, 035304 (2010).

[24] B. Hof, A. de Lozar, D. J. Kuik, and J. Westerweel, Repeller or attractor? selecting the dynamical model for the onset of turbulence in pipe flow, Phys. Rev. Lett. **101**, 214501 (2008).

[25] D. Barkley, Theoretical perspective on the route to turbulence in a pipe, J. Fluid Mech. **803**, P1 (2016).

[26] G. Bonnavion and O. Cadot, Unstable wake dynamics of rectangular flat-backed bluff bodies with inclination and ground proximity, J. Fluid Mech. **854**, 196–232 (2018).

[27] V. Parezanović, R. Monchaux, and O. Cadot, Characterization of the turbulent bistable flow regime of a 2d bluff body wake disturbed by a small control cylinder, Exp Fluids **56**, 12 (2015).

[28] O. Benois, C. Landim, and M. Mourragui, Hitting times of rare events in markov chains, J. Stat. Phys. **153**, 967 (2013).

[29] M. Farazmand and T. P. Sapsis, Extreme events: Mechanisms and prediction, Applied Mechanics Reviews **71**, 050801 (2019).

[30] T. P. Sapsis, Statistics of extreme events in fluid flows and waves, Annual Review of Fluid Mechanics **53**, 85 (2021).
[31] A. C. M. Freitas, J. M. Freitas, and M. Todd, Hitting time statistics and extreme value theory, Probab. Theory Related Fields 147, 675 (2010).

[32] V. Lucarini, D. Faranda, J. M. M. de Freitas, M. Holland, T. Kuna, M. Nicol, M. Todd, S. Vaienti, et al., *Extremes and recurrence in dynamical systems* (John Wiley & Sons, 2016).

[33] V. Araújo and H. Aytaç, Decay of correlations and laws of rare events for transitive random maps, *Nonlinearity* 30, 1834 (2017).

[34] E. ElJack, High-fidelity numerical simulation of the flow field around a naca-0012 aerofoil from the laminar separation bubble to a full stall, *Int. J. Comput. Fluid Dyn.* 31, 230 (2017).

[35] Y. A. ElAwad and E. M. ElJack, Numerical investigation of the low-frequency flow oscillation over a naca-0012 aerofoil at the inception of stall, *International Journal of Micro Air Vehicles* 11, 1756829319833687 (2019).

[36] L. Meirovitch, Analytical methods in vibrations., (1967).

[37] E. Achenbach, Experiments on the flow past spheres at very high reynolds numbers, *Journal of fluid mechanics* 54, 565 (1972).

[38] E. Achenbach, The effects of surface roughness and tunnel blockage on the flow past spheres, *Journal of fluid mechanics* 65, 113 (1974).

[39] S. H. Rudy and T. P. Sapsis, Prediction of intermittent fluctuations from surface pressure measurements on a turbulent airfoil, *AIAA Journal* 60, 4174 (2022).

[40] T. Ohtake, Y. Nakae, and T. Motohashi, Nonlinearity of the aerodynamic characteristics of naca0012 aerofoil at low reynolds numbers, *Japan Society of Aeronautical Space Sciences* 55, 439 (2007).

[41] G. Hristov and P. J. Ansell, Poststall hysteresis and flowfield unsteadiness on a naca 0012 airfoil, *AIAA J.* 56, 2528 (2018).

[42] T. Nemoto and A. Alexakis, Do extreme events trigger turbulence decay? – a numerical study of turbulence decay time in pipe flows, *J. Fluid Mech.* 912, A38 (2021).

[43] S. Gomé, L. S. Tuckerman, and D. Barkley, Extreme events in transitional turbulence, *Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci.* 380, 20210036 (2022).

[44] D. Faranda, M. Bourgoin, S. Miralles, P. Odier, J.-F. Pinton, N. Plihon, F. Daviaud, and B. Dubrulle, Robust estimate of dynamo thresholds in the von kármán sodium experiment using the extreme value theory, *New Journal of Physics* 16, 083001 (2014).

[45] T. Dessup, C. Coste, and M. Saint Jean, Hysteretic and intermittent regimes in the subcritical bifurcation of a quasi-one-dimensional system of interacting particles, *Phys. Rev. E* 93, 012105 (2016).

[46] T. Schikarski, H. Trzenschiok, W. Peukert, and M. Avila, Inflow boundary conditions determine t-mixer efficiency, *Reaction Chemistry & Engineering* 4, 559 (2019).

[47] D. Faranda, (private communication), .