Bifurcation scenario in the two-dimensional laminar flow past a rotating cylinder

Abstract The aim of this paper is to provide a complete description of the bifurcation scenario of a uniform flow past a rotating circular cylinder up to $Re = 200$. Linear stability theory is used to depict the neutral curves and analyse the arising unstable global modes. Three codimension-two bifurcation points are identified, namely a Takens–Bogdanov, a cusp and generalised Hopf, which are closely related to qualitative changes in orbit dynamics. The occurrence of the cusp and Takens–Bogdanov bifurcations for very close parameters (corresponding to an imperfect codimension-three bifurcation) is shown to be responsible for the existence of multiple steady states, as already observed in previous studies. Two bistability regions are identified, the first with two stable fixed points and the second with a fixed point and a cycle. The presence of homoclinic and heteroclinic orbits, which are classical in the presence of Takens–Bogdanov bifurcations, is confirmed by direct numerical simulations. Finally, a weakly nonlinear analysis is performed in the neighbourhood of the generalised Hopf, showing that above this point the Hopf bifurcation is subcritical, leading to a third range of bistability characterised by both a stable fixed point and a stable cycle.


Introduction
The flow past a circular cylinder is a classical configuration which has been widely adopted in the fluid dynamics community as a canonical model to investigate vortex shedding behind bluff bodies. In the case of a fixed cylinder, i.e. without rotation, the dynamics and the corresponding bifurcations are well known (Williamson 1996). The case of a rotating cylinder, which has implications for flow control using wall motion (Modi 1997;el Hak 2000), has recently received attention. A number of numerical studies in a two-dimensional framework have been conducted (Kang, Choi & Lee 1999;Stojković, Breuer & Durst 2002, 2003Mittal 2004) and have revealed the existence of several steady and unsteady regimes. Linear stability approaches (Pralits, Brandt & Giannetti 2010;Pralits, Giannetti & Brandt 2013) have shown the existence of two separated regions of instability in the (Re, α) plane, where α is the dimensionless rotation rate and Re is the Reynolds number. The so called Mode I becomes unstable via a supercritical Hopf bifurcation and it is present for 0 ≤ α ≤ 2. This mode is the one associated with the classical Bénard-von-Kármán vortex street, and characterised by the alternate shedding of vortices of opposite sign. At higher rotation rates, around 4.5 ≤ α < 6 another unsteady mode exists, denoted as Mode II. The physical mechanism driving this mode is rather different, as it corresponds to a slow-frequency shedding of vortices with the same vorticity sign. Its onset is less well characterised than Mode I from the point of view of bifurcation theory: the fact that the frequency is very low suggests a more complex bifurcation scenario and its supercritical or subcritical nature is still unclear. The full characterisation of Mode II is complicated by the fact that, in approximately the same range of (Re, α) parameter space, a region where three steady-state solutions coexist has been evidenced (Pralits et al. 2010;Rao et al. 2013a). A more thorough characterisation of this phenomenon has been carried out by Thompson et al. (2014) who observed that the region of existence of multiple steady-state solutions grows with the Reynolds number. Note also that the picture is further complicated by the existence of three-dimensional (3-D) instabilities in this range. This point is outside of the range of the present paper which restricts to 2-D dynamics, but a brief review on 3-D stability properties of this flow can be found in appendix E.
To explain the existence of multiple steady states, Rao et al. (2013a) conjectured that they emerge from a cusp bifurcation point. Indeed, a cusp correctly explains the change in the number of steady states from one to three. However, a cusp is not generally associated with the existence of a Hopf bifurcation in the same range of parameters, so it cannot explain, alone, all the features discussed above. The fact that the frequency of Mode II is very small is an indicator of a second kind of codimension-two bifurcation, namely a 0 2 or Takens-Bogdanov bifurcation (Kuznetsov 2013, chapter 8, p. 314) This bifurcation typically occurs when the frequency of a limit cycle vanishes. However, in the vicinity of a standard Takens-Bogdanov bifurcation, only two steady states generally exist, not three. This combination of features suggests that the picture could hide a codimension-three bifurcation point, also known as a generalised Hopf bifurcation. The unfolding of this generalised Takens-Bogdanov bifurcation has been studied by Dumortier et al. (2006) and Kuznetsov (2005) from a mathematical point of view, but to our knowledge such a feature has not yet been evidenced in a fluid dynamics system such as the one considered here.
The main purpose of the present work is to review the classification of the possible 2-D states in the (Re, α) ∈ [0, 200] × [0, 10] parameter plane with the point of view of dynamical system theory. Firstly, we will characterise the nature of the codimension-one bifurcation curves (Hopf or saddle nodes). We give a cartography of the regions where multiple steady states exist and give a detailed description of these multiple states as well as their stability properties. We further identify three codimension-two points, namely a Takens-Bogdanov (TB) bifurcation, a cusp and a generalised Hopf (GH) bifurcation. We show that the two first are effectively located very close to each other and that the whole dynamics in this range of parameters is effectively described by the unfolding of a codimension-three bifurcation point.
The article is organised as follows: in § 2 the formulation of the problem is discussed together with the methodology adopted in the present analysis. Section 3 begins with a characterisation of the multiple steady states. A complete bifurcation diagram covering the range (Re, α) ∈ [0, 200] × [0, 10] is then presented. The next subsections aim at clarifying the picture in the vicinity of the identified codimension-two points.

Geometrical configuration and general equations
The two-dimensional flow past a rotating circular cylinder is controlled by two parameters: the Reynolds number Re = U ∞ D/ν and the rotation rate α = ΩD/2U ∞ . Here, Ω is the dimensional cylinder angular velocity, U ∞ is the free stream velocity, D the diameter of the cylinder and ν the dynamic viscosity of the fluid. The fluid motion inside the domain is governed by the two-dimensional incompressible Navier-Stokes equations, where U is the velocity vector whose components are (U, V), P is the reduced pressure and the viscous stress tensor τ (U) can be expressed as ν(∇U + ∇U T ). The incompressible Navier-Stokes equations (2.1) are complemented with the following boundary conditions: on the cylinder surface, no-slip boundary conditions are set by U · t = ΩD/2 and U · n = 0, where (t, n) are the director vectors of the surface in the plane (x, y); in the far field, uniform boundary conditions are set U → (U ∞ , 0) when r → ∞, where r is the distance to the cylinder centre (see figure 1). In the discussion we consider clockwise rotation of the cylinder surface (α > 0). In the following, Navier-Stokes equations (2.1) and the associated boundary conditions will be written symbolically under the form B(∂Q/∂t) = N S(Q), where Q = (U, P) is the state vector and B is a linear projection operator, meaning that the time derivatives apply only on the velocity components.

Linear stability analysis
Under the framework of linear stability analysis, we first need to identify base-flow solutions defined as the steady solutions Q b of the (two-dimensional) Navier-Stokes equations, namely the solutions of N S(Q b ) = 0. We then characterise the dynamics of small-amplitude perturbations around this base flow by expanding them over the basis of linear eigenmodes, i.e. (2.2) Here, is a small parameter, λ j the eigenvalues andq j the eigenmodes. The eigenpairs [λ j ,q j ] have to be determined as the solutions of the following eigenvalue problem: Which will be written in the following under the symbolic form λ j Bq j + LN Sq j = 0. In the following we consider that eigenmodesq(x, y) have been normalised, see appendix C for further details. Note that in (2.2), to fully represent the dynamics, the summation over eigenmodes may involve a continuous sum over the spectrum, i.e. the discrete and the continuous or essential spectra of the operator (see Kapitula & Promislow (2013) for a rigorous discussion). However, to determine global stability we only need to consider a limited number of eigenmodes, so we keep the summation as a discrete sum indexed by j.
Owing to the eigenvalues, two cases can be distinguished: (i) If all eigenvalues λ j have negative real part the considered base flow is a stable solution. (ii) If n eigenvalues have positive real part, the considered base flow will be referred to as a n-unstable solution. Note that 1-unstable solutions are commonly referred to as saddle points because a projection of their dynamics in a 2-D plane (phase portrait) has an attractive direction and another repulsing one, while 2-unstable solutions are either unstable nodes or unstable foci depending if the leading eigenvalues are both real or complex conjugates.
The transition from stable to unstable (or from n-unstable to n + 1-unstable) is called a local bifurcation. The simplest bifurcations (such as saddle nodes and Hopf) are said to be codimension-one and occur along given curves in the parameter plane (Re, α). The intersection of two such curves tangentially is called a codimension-two bifurcation and generally leads to a rich dynamics in the vicinity of the intersection point.

Notions of bifurcation theory
From the viewpoint of dynamical system theory, the expression (2.2) can be generalised as a decomposition of the perturbations over the leading modes of the system (2.4) Then, the problem can be reduced to a low-dimensional system governing the amplitudes where (NL) represent the nonlinear interactions between modes. Investigation of these nonlinear terms allows us to predict the dynamics in the vicinity of bifurcation points. Systematic methods exist to compute these nonlinear terms (such as weakly nonlinear expansions, centre manifold reduction or Lyapunov-Schmidt reduction). However, restricting ourselves to a qualitative point of view (up to a continuous change of coordinates with continuous inverse), it is also possible to predict a number of features by examining the generic normal form of the bifurcation, namely, a standard form to which the dynamical system can be reduced by a series of elementary manipulations (see Wiggins (2003) for details). Particular forms of codimension-two bifurcations encountered in the rotating cylinder are discussed in § § 3.5 and 3.6.

Numerical methodology
In the present manuscript, we adopt the same numerical methodology used in Fabre et al. (2020) and described in Fabre et al. (2019). The computation of the steady-state solutions, the resolution of the linear problems and the time stepping techniques are implemented using the open-source finite element software FreeFem++. Parametric studies and generation of figures are performed using Octave/Matlab thanks to the generic drivers of the StabFem project (see a presentation of these functionalities in Fabre et al. 2019). According to the philosophy of this project, codes reproducing parts of the results of the present paper are available from the StabFem website (https://gitlab.com/stabfem/ StabFem). On a standard laptop, all the computations discussed below can be obtained in a few hours, except time stepping simulations which take longer. Results presented in § 3 are obtained with a computational domain L x = 120 and L y = 80 in the streamwise and cross-stream directions, respectively. The cylinder centre is located 40 diameters downstream of the inlet, symmetrically between the top and bottom boundaries. Numerical convergence issues are discussed in appendix D by a meticulous comparison between results obtained with different meshes, where domain dimension and grid density were varied.
Steady nonlinear Navier-Stokes equations are solved by a Newton method. In the degenerated cases, pseudo-arc length continuation is performed to be able to compute multiple steady-state solutions, as described in appendix A. The generalised eigenvalue problem (2.3) is solved by the Arnoldi method or by a simple inverse iteration algorithm. Finally, nonlinear unsteady Navier-Stokes equations are integrated forward in time with a second-order time scheme (Jallas, Marquet & Fabre 2017).

Characterisation of multiple steady-state solutions
To introduce the existence of multiple steady states, we first characterise them by plotting in figure 2 the associated lift as function of the rotation rate α, for four different values of α. In these plots, stable solutions are indicated by continuous lines and unstable ones by dashed lines, following the usual convention in dynamical systems theory.
For Re = 60, as illustrated in figure 2(a), only one steady state exists for all values of α, for Re = 60. This state is stable except in the ranges α 2 (corresponding to the existence of Mode I), and 5.2 α 5.5 (corresponding to the existence of Mode II).
For higher Reynolds numbers, a small region of multiple solutions arises in a small-scale interval around α ≈ 5. This phenomenon is illustrated in figure 2(b) for Re = 100 and is associated with an 's' shape of the curve, featuring two successive folds. Note that, before the first fold, the steady solution is 2-unstable (focus type); at the first fold it turns into 1-unstable (saddle type) and at the second fold it turns into stable. To detect these folds, pseudo-arc length continuation is carried out with α as a parameter and the horizontal force exerted on the cylinder surface F x as a monitor to track and distinguish multiple steady states (see appendix A for a more detailed discussion).
For larger values of the Reynolds number, as illustrated in figure 2(c) for Re = 170, the interval of existence of multiple states for α ≈ 5 expands to α ∈ [4.75, 5.12]. In addition, we observe a second range displaying multiple states for α > 5.87. This second interval is associated with a fold bifurcation at α = 5.87, giving rise to two additional and disconnected steady solutions. Note that both these solutions are unstable, respectively of node and saddle types.
Finally, for Re = 200, as illustrated in figure 2(d), we observe that the two ranges of multiple steady states are merged into a single one. In this case there is a single saddle-node bifurcation around α = 4.75 leading to two branches of steady states which are disconnected from the branch existing for lower values of α. Here, one of these branches is stable and the second is unstable (saddle type).

Topological description of steady-state solutions
We now illustrate the spatial structure of some steady-state solutions, with emphasis on the topological structure of the corresponding flows. We restrict ourselves to the case Re = 200, as previously considered in figure 2(d). Figure 3(a) corresponds to α = 1.8, the value at which Mode I is re-stabilised. The corresponding flow is characterised by a stagnation point located beneath the cylinder axis, on the left side of the y-cylinder axis. Compared to the steady flow in the non-rotating case, which is characterised by a symmetrical recirculation region, the upper recirculating bubble is reduced whereas the lower one is moved downwards.
Further increasing the rotation speed, both recirculation bubbles shrink and eventually vanish. At α = 4.35 (figure 3b) corresponding to the lower threshold for the existence of Mode II, recirculating bubbles have already disappeared and the vorticity wraps the cylinder. Stagnation point is located on the opposite side but downstream the cylinder vertical axis.   Compared to the previous state, the flow is topologically different as no stagnation point is observed along the wall of the cylinder. On the other hand, two stagnation points are observed within the flow. One of them is elliptic and located at the centre of the detached recirculation bubble. The other is hyperbolic and located along the streamline bounding the recirculation bubble. Figure 3(d-f ) displays the three coexisting steady states at α = 5.25 and Re = 200. The topology of the streamlines of unstable and stable steady states differs. In the stable case (panel d) there is a single recirculation region encircling the cylinder and bounded by a hyperbolic stagnation point, as in the classical potential solution existing in this range of rotation rates. On the other hand, for both unstable states, the topology is similar to the case of figure 3(c). The recirculation region is detached from the cylinder and contains an elliptic stagnation point located approximately in the midpoint between the hyperbolic point and the bottom point of the cylinder surface. In the unstable steady state, the recirculating region is more stretched, as it can be seen in figure 3(d-f ).
We highlight that even though topological changes in the streamlines of the steady states and bifurcations of the velocity field are in general independent events (see Brøns 2007), in some cases these two events occur in a small neighbourhood of the space of parameters (see Heil et al. 2017). In the current situation it has been confirmed that there is not a one-to-one relation between both phenomena. For instance, the transition between detached recirculation bubble (as in panel c) and recirculation bubble encircling the cylinder (as in panel d) along the stable branch occurs at some value of α in the range [4.75-5.25] where no dynamical bifurcation occurs. Yet, for larger Reynolds numbers, i.e. Re 190, successive creation and destruction of vortices seems to be relevant in the preservation of the disconnected branch of steady states.

Analysis of the spatial structure of direct and adjoint eigenmodes
To explain why the steady state displayed in figure 3( f ) is unstable, the two corresponding unstable modes (both associated with real eigenvalues) are displayed in figure 4 for Re = 200 and α = 5.25. Direct modes are characterised by two recirculating regions of opposite vorticity. Vorticity is stronger and more localised in Mode IIa while Mode IIb displays a larger region with non-zero vorticity. Adjoint eigenvectorsq † for Mode IIa and Mode IIb are also displayed in figure 4. Adjoint fields (Luchini & Bottaro 2014) can be interpreted as a kind of Green's function for the receptivity of the global mode. Scalar product of the adjoint field with a forcing function or an initial condition provides the amplitude of the instability mode (see Giannetti & Luchini 2007). Therefore, Mode IIa is highly receptive in the upper right side of the near wake of the cylinder. The region of maximum receptivity extends from the close upper right region of the cylinder to a larger region at the bottom right of the cylinder and it is weaker than Mode IIa. Both modes present weak sensitivity to forcing upstream of the cylinder.

Bifurcation diagram in the parameter plane (Re, α)
The bifurcation curves detected in the α < 10, Re < 200 range by linear stability analysis of all steady-state solutions are depicted in figure 5.
Three Hopf bifurcation curves are detected and plotted with full lines. The first one encircles the range of existence of unsteady Mode I. The second one delimits the range of existence of unsteady Mode II in its lower and left parts, but not on its upper part. The third one (in grey) occurs along a steady state which is already unstable, and hence is not likely to be related to a bifurcation observable in DNS or experiments.
In addition, we have identified two bifurcation curves associated with saddle nodes or 'folds', here denoted F + and F − . These curves delimit the range of existence of multiple two-dimensional steady states, displayed as a grey region in figure 5. Note that the extension of this region explains the difference between the cases Re = 170 and Re = 200 discussed in the previous paragraph; according to the figure a single interval of α is found for Re 190.
In figure 5, the two fold curves seem to merge with the Hopf curve existing for lower Re at a point with coordinates Re ≈ 75, α ≈ 5.4. Inspection shows that there are actually both a 0 2 or TB bifurcation and a cusp (C) bifurcation in very close vicinity in this range of parameters. This region will be studied in § 3.5. Additionally, in another range of parameters located at the lower threshold of existence of the Mode II, we have identified the existence of a Bautin or GH bifurcation which splits the Hopf curve into supercritical (Re < Re GH ) and subcritical (Re > Re GH ). This region will be studied in § 3.6.

Qualitative study of the normal form
The transition occurring for Re ≈ 75 and α ≈ 5.4 is characterised by the end of the Hopf curve (H − ) at a fold curve (F + ) (characteristic of a Takens-Bogdanov bifurcation), and a transition between one and three steady states (characteristic of a cusp). This suggests that the present situation is actually very close to a codimension-three bifurcation. The dynamical behaviour of the system can thus be expected to be well predicted using the normal form describing the universal unfolding of the codimension-three planar bifurcation, also called a generalised TB bifurcation. This normal form has been studied by both Dumortier et al. (2006) and Kuznetsov (2013, chapter 8.3). The normal form can be written as follows: where β 1 , β 2 and β 3 are unfolding parameters (mapped from the physical parameters (Re, α)), c 1 , (which can be rescaled to ±1) are fixed coefficients which depend on the nonlinear terms of the underlying system. Note that this normal form generalises both the normal form of the standard TB bifurcation (which is recovered for β 1 (Re, α) = 0) and the one of the fold bifurcations (which is recovered for β 3 (Re, α) = 0). The occurrence of both these codimension-two conditions for very close values of the parameters is characteristic of an imperfect codimension-three bifurcation and justifies the relevance of the associated normal form.
The dynamics of the normal form (3.1) has been explored by Dumortier et al. (2006) who classified the possible phase portraits and the associated bifurcation diagrams as functions of the unfolding parameters (β 1 , β 2 , β 3 ) along a spherical surface. They showed that all possible bifurcation diagrams fall into three possible categories, called focus, saddle and node according to the values of the coefficients c 1 and . The situation 0 < c 1 < 2 √ 2 and = −1 corresponds to the stable focus case and is found to lead to a bifurcation diagram consistent with the present situation, so we concentrate on this case. Figure 6 illustrates all the possible behaviours of the dynamical system, sketched by sample phase portraits, along with their range of existence in the (β 1 , β 2 ) plane. This figure corresponds to a subset of the complete diagram displayed in Dumortier et al. (2006, chapter 1, pp. 6-8), restricted to a range of parameters which is sufficient to explain all the dynamical features of the present problem. The bifurcation diagram displays two codimension-two points, a cusp C and a TB. These codimension-two points result from the tangential intersection of two codimension-one curves: the cusp point C occurs when the two fold curves F + and F − collide, while the TB point arises from the intersection of the supercritical H − Hopf curve and the F + fold. In addition, the bifurcation diagram predicts a homoclinic global bifurcation along a curve H ∞ originating from the TB point and terminating along the F − fold on a point denoted SNL (for saddle-node loop). Left from this point, the F − curve corresponds to a local saddle node while right from this point it corresponds to a homoclinic saddle-node bifurcation (appearance of two fixed points along a previously existing cycle). Note that the SNL point and the intersection of H − and F − are formally not codimension-two points (see Dumortier et al. 2006).
Phase portraits obtained in the various regions delimited by bifurcation boundaries are displayed in the panels of figure 6. One of the most interesting predictions is the existence of two regions characterised by the existence of two stable states, a bistability phenomenon. The first region (3), in the vicinity of the cusp, is characterised by two stable  (1), (2), (3), (4),  (4) is characterised by both a stable steady state and a stable cycle. In all other regions, there is a single stable solution which is either a steady state (in regions 1 and 5) or a cycle (in region 2). Note that in these phase portraits nodes and foci are not distinguished. Distinguishing between these cases (Dumortier et al. 2006) leads to consideration of a larger number of subcases (for instance region 1 could be split in two subregions corresponding to a stable node and a stable focus . . . ) but the transitions between these subcases are not associated with bifurcations.

Numerical results in the C-TB region
In order to check the predictions of the normal form approach, we have conducted an accurate exploration of the range of parameters corresponding to the C-TB region. The exploration allowed us to confirm the existence of both a cusp and a Takens-Bogdanov point. The locations in the (α, Re) plane are given in table 1. Figure 7 displays 'zooms' of the full bifurcation diagram (figure 5) in two narrow ranges centred on the C and TB codimension-two points. The bifurcation curves and the regions  are numbered with the same convention as in figure 6. Although it is not possible to present all results in a single figure because the curves are very steep and close to each other, the numerical results fully confirm the predictions of the normal form. In particular, the numerical results allow us to confirm the coexistence of two stable states (in regions 3) and of a stable cycle and a stable state (in region 4). However, a precise mapping of the curve H ∞ bounding the region 4 could not be achieved, but the occurrence of a global homoclinic bifurcation was confirmed (see § 3.5.3).

Homoclinic bifurcation
As explained in § 3.5, the normal form predicts a homoclinic curve H ∞ and a homoclinic saddle-node bifurcation along the F − curve, right from the SNL point, corresponding to the appearance of two steady solutions along a previously existing cycle. A generic feature of the imminent presence of a homoclinic saddle-node bifurcation is the divergence of the period of the limit cycle on which the saddle node appears. More precisely, the period is expected to scale as ∝ 1/ √ α SN − α as α → α SN (see Gasull, Mañosa & Villadelprat 2005). To check this prediction, time stepping simulations were conducted for Re = 170 and values of α just below the F − curve. As shown in figure 8 the period of the limit cycle effectively diverges as one approaches the bifurcation following the theoretical behaviour. Dynamics near the threshold can be perfectly understood in a two-dimensional manifold. Phase portraits of the bifurcation are displayed in figure 9. These phase portraits were computed with an initial guess generated by a small linear perturbation to a steady state in the direction of its corresponding eigenmode. The initial guess is then integrated in time until it reaches its limit set, i.e. a periodic, homoclinic orbit or another steady state. Below the bifurcation threshold (figure 9a) a stable limit cycle exists, represented by a thick solid line. At the bifurcation threshold, a saddle node arises along this cycle, which ceases to exist, giving rise to a homoclinic connection (an approximation of this orbit is delineated by a thick solid line in figure 9b). Beyond the saddle-node bifurcation, the saddle node splits into two fixed points. Hence, three steady states exist, including a stable one (see figure 9c). There exist four stable heteroclinic connections, two between unstable-stable steady states represented by a dashed line in figure 9(c) and other two between saddle-stable steady states denoted by a solid line. This sequence of events is fully consistent to the sequence connecting phase portraits (2), (SNL) and (4) in figure 6.

Normal form analysis
Bautin bifurcation or GH is a codimension-two bifurcation where the equilibrium has purely imaginary eigenvalues λ 1,2 = ±ω 0 with ω 0 > 0, and the third-order coefficient of the normal form vanishes. Generalised Hopf bifurcation is thus a degenerate case of the generic Hopf bifurcation, where the cubic normal form is not sufficient to determine the nonlinear stability of the system. To unravel the dynamics near the Bautin bifurcation point consider the normal form (3.2) Three curves are of special interest: (i) System (3.2) undergoes a supercritical Hopf bifurcation in the half-line H + = {(β 1 , β 2 )|β 2 > 0, β 1 = 0}. This curve separates a region containing a stable focus to a region containing an unstable focus plus a stable limit cycle. (ii) System (3.2) undergoes a subcritical Hopf bifurcation in the half-line H − = {(β 1 , β 2 )|β 2 < 0, β 1 = 0}. This curve separates a region containing an unstable focus, from one containing a stable focus and two limit cycles (one being stable and the other one being unstable). (iii) System (3.2) undergoes a fold cycle bifurcation on the curve F LC = {(β 1 , β 2 )|β 2 1 + 4β 2 = 0, β 1 < 0}. This curve separates a region containing two limit cycles from one which does not contain any limit cycle (a stable fixed point also exists in both regions).
The most notable feature of this bifurcation is the existence of a bistability region characterised by two stable states (a fixed point and a cycle). Therefore, hysteretic behaviour is expected as one successively crosses curves H − and F LC . The bistability range is also characterised by the existence of an unstable limit cycle constituting the 'edge state' bounding the basins of attraction of the two stable states (figure 10).

Weakly nonlinear analysis
Unstable limit cycles are not easy to track, since they require stabilisation techniques, such as BoostConv  or edge-state tracking (Bengana et al. 2019), or the use of continuation techniques, such as harmonic balance (Fabre et al. 2019). Alternatively, we have performed a multiple-scale analysis up to fifth order (see appendix C). This method was previously used to study thermoacoustic bifurcations in the Rijke tube (Orchini, Rigas & Juniper 2016), displaying a good match with time stepping simulations with a much lower computational cost. By performing a weakly nonlinear analysis up to fifth order it is possible to determine a complex amplitude equation for the amplitude A of the critical linear modeq. Here, the critical linear mode is normalised so that its L 2 B-norm (see appendix C), i.e. its kinetic energy, is unity, which corresponds to the same normalisation as in Mantič-Lugo, Arratia & Gallaire (2014 We remark that (3.3) is equivalent to (3.2) if separating real and imaginary parts. Searching for a solution under the form A = |A| e iωt , and injecting into (3.3) leads to where ν 1 = ν 1,0 + 2 ν 1,1 , λ = 2 λ 0 + 4 λ 1 , ν 2 = ν 2,0 and subscripts r, i denote real and imaginary parts respectively. It turns out that ν 2,r is always negative while ν 1,r changes sign at (Re, α) = (Re GH , α GH ). One can deduce the following consequences: (i) If Re < Re GH (i.e. ν 2r < 0), (3.4) has a single solution |A| for λ r > 0 (i.e. Re > Re c ) and none for λ r < 0 (i.e. Re < Re c ). In this case, the Hopf bifurcation is supercritical. (ii) If Re > Re GH , (i.e. ν 2r > 0), (3.4) has a single solution |A| for λ r > 0 (i.e. Re > Re c ), two solutions if λ c < λ r < 0 with λ c = ν 2 1,r /4ν 2,r and no solution if λ r < λ c . In this case, the Hopf bifurcation is subcritical. The condition λ r = λ c defines a curve in the (Re, α) plane which corresponds to the fold cycle bifurcation associated with the emergence of the two limit cycles. Figure 11 represents the amplitude and frequency of the limit cycles predicted by (3.4) for three values of Re. According to these results, the fold curve is predicted to be very close to the Hopf curve, i.e. within a few tenths of Re up to Re = 250. This behaviour allows us to clarify the transition occurring at the GH point in figure 6. For Re < Re GH , when increasing Re for fixed α (or increasing α with fixed Re), the transition occurs via a supercritical Hopf bifurcation. On the other hand, for Re > Re GH , the transition is predicted to be subcritical, involving the existence of a band where both steady state and Mode II coexist. Note that the width of the bistability band predicted by the weakly nonlinear analysis is very narrow, and could thus be difficult to evidence using direct numerical simulations.

Conclusion and discussion
The present study allowed us to clarify the bifurcation scenario in the two-dimensional flow past a rotating cylinder, especially concerning the range of parameters corresponding to the onset of the 'Mode II' unsteady vortex shedding mode. Using steady-state calculation involving arclength continuation and linear stability analysis, we have been able to draw all bifurcation curves existing in the range of parameters corresponding to Re < 200 and α < 5. Three codimension-two bifurcations have been identified along the border of the range of existence of this mode, namely a Takens-Bogdanov, a cusp and a generalised Hopf. The first two are located in close vicinity, in such a way that the whole dynamics can be understood using the normal form of the codimension-three bifurcation (for a generalised Takens-Bogdanov bifurcation). The analysis also allowed us to identify three ranges of parameters characterised by bistability, two of them located in the vicinity of the Takens-Bogdanov and cusp points, the third one emanating from the generalised Hopf point. Time stepping simulations and a weakly nonlinear analysis have confirmed these findings, and have also allowed us to characterise the homoclinic and heteroclinic orbits connecting the fixed points, in full accordance with the predictions of the normal form theory. The most surprising result of the study is the existence of an almost perfect codimension-three bifurcation in a problem characterised by only two control parameters. Such a feature suggests that the problem could be quite sensible to any small perturbations in a way such that small perturbations could completely change the scenario. We have checked that the scenario is robust with respect to numerical discretisation issues (see appendix D). The dependency with respect to additional physical parameters is more interesting. The effect of compressibility is an interesting question which we expect to investigate in future studies. Preliminary results have shown that for a Mach number of order 0.1, the dynamics in the region of the near-codimension-three point is effectively greatly modified. Other additional parameters, such as for instance shear or confinement, could be added. Finally, one may question the relevance of the present findings for three-dimensional flows. A short review of three-dimensional stability properties of the rotating cylinder flow is given in appendix E. The discussion confirms that the most important results of the present study occur in range of parameters where no three-dimensional instabilities are present.

Declaration of interests
The authors report no conflict of interest.

Appendix A. Pseudo arc-length continuation
Arc-length continuation is a standard technique in dynamical systems theory. It allows for the continuation of a given solution branch through a turning or fold point. At the turning point the Jacobian of the system is singular; therefore, any iterative method based on the Jacobian is doomed to failure. To prevent the stall in the convergence of the Newton's method, an extra condition needs to be added to the system of equations. In the current study we have chosen a pseudo arc-length methodology, which is based in a predictor-corrector strategy. The extended system adds an extra equation which ensures the tangency to the branch of the solution. For that purpose, a parameter is chosen, here either Re or α, and a monitor of the variation, either the horizontal force acting on the cylinder surface F x or the vertical force F y . The parameter and the monitor are parametrised by the length of the branch, here indicated by the parameter s. The current solution is varied by a given step Δs tangent to the solution branch and later corrected by a orthogonal correction. Let us denote by the subscript j the arc-length iteration and by the superscript n the Newton iteration of the corrector step, where N is used to denote the last step. In the description below, let us consider without loss of generality we have fixed the parameter α and the monitor F x .

A.1. Predictor
The predictor step consists in the determination of a initial guess α 0 j for the iteration j of the arc length. The initial guess is determined from a tangent extrapolation of the solution branch.
In (A 1), dα N j−1 /ds is the slope of the tangent in the α direction and dq N j−1 /ds in the direction of the vector field. The tangent is computed from the differentiation of the stationary Navier-Stokes equations (2.1) where we have used the notation NS q N j−1 = 0 to denote the steady incompressible Navier-Stokes equation whose solution is q N j−1 . The tangent is completed with a normalisation condition in the arc length A.2. Corrector This step consists in an orthogonal correction of the tangent guess. To do so one needs to solve the following system of equations ⎡ , (A 4) where the last equation of (A 4) comes from the differentiation of the normalisation condition (A 3) and considering that Δα j = Δα n+1 (see Gallaire et al. 2016) and to determine the validity of stability analysis on the mean flow (see Sipp & Lebedev 2007). In this article WNL analysis is used to determine the existence of a generalised Hopf bifurcation (see § 3.6). The starting point of the weakly nonlinear method is the decomposition of the flow field into multiple scales Q = Q b + A wnlq e iω 0 t + c.c. + 2 q 2,0 + |A wnl |q |A wnl | 2,0 + (A 2 wnl q 2,2 e 2iω 0 t + c.c.) If we take into account the definition of the slow time scale τ = 2 t, the fact that up to leading order O( ) we have dA wnl /dt = iω 0 A wnl and we define a new amplitude which depends on as A = A wnl we can rewrite (C 2) as dA dt = (iω 0 + 2 λ 0 + 4 λ 1 )A + (ν 1,0 + 2 ν 1,1 )|A| 2 A + ν 2,0 |A| 4 A.
In the following we consider that the eigenmodeq and its adjointq † have been normalised so that ||q|| 2 B = q, Bq = û,û = 1 and q † , Bq = û † ,û = 1. This normalisation is the same as that one used in the self-consistent methodology (see Mantič-Lugo et al. 2014): with this choice, A is a real constant representing the amplitude of the linear mode with respect to its L 2 norm. In the following we will use the notation LN S iω q = iωBq − LN Sq to denote the application of the linearised operator at a specific frequency ω.

Appendix D. Mesh convergence
Mesh independence of the solutions has been verified systematically. First, we have considered a given mesh refinement and varied the physical size of the domain, see table 2. We have observed that for a domain length of 80 diameters downstream the cylinder centre, 40 diameters upstream the cylinder centre and 40 in the cross-stream direction the solution is not affected by the imposition of boundary conditions. Secondly, we have looked at the effect of the mesh refinement on the properties of the solution. For that purpose a parametric study of eigenvalues, Hopf WNL coefficients and global monitors of a given steady-state solution have been carried out, see (table 3). The sensitivity to mesh convergence of cusp and Takens-Bogdanov bifurcation points has been also tested. Results show that each of them is found within ΔRe c < 0.2. Every mesh is computed by Delaunay triangulation. Mesh M 1 has been generated by blocks, as it is generally done with structured meshes; M 2 and M 3 have been computed following the mesh adaption procedure described in Fabre et al. (2019, appendix A), with respect to base flow only and with respect to base flow and direct mode structure; M 4 and M 5 are the consequence of successive division of each triangle edge by two and four respectively, with respect to mesh M 3 . The mesh selected for this study is M 1 which provides results within the one per cent of relative error with respect to the finest mesh. One of the reasons that led us not to use mesh adaptation is the fact that the structure of the mode greatly changes within  Appendix E: Three-dimensional stability of steady-state solutions In this section, we review three-dimensional stability studies carried out by Pralits et al. (2013), Rao et al. (2013a,b), Radi et al. (2013) and Rao et al. (2015).
It is now well known the secondary three-dimensional transition from a two-dimensional unsteady flow towards a three-dimensional flow at Re ≈ 190 and α = 0, see Williamson (1996). Vortices in the wake of the fixed cylinder, i.e. α = 0, develop spanwise waviness whose wavelength is approximatelyfour cylinder diameters. The rotation of the cylinder surface on this linear steady mode, denoted as Mode A in Rao et al. (2015), has a stabilising effect for rotation rates α < 1, see figure 12.
Instead, if we consider the stability of an infinitesimal spanwise perturbation on a steady-state solution, the flow displays spanwise waviness at a much lower Reynolds number Re ≈ 100 and α = 0. The onset of instability of this stationary mode, denoted as Mode E in Rao et al. (2015), is shown in figure 12 as a function of (Re, α).
In the same region of existence of the unsteady two-dimensional Mode II, experimental evidence has shown the presence of a three-dimensional mode, see Linh (2011). A steady three-dimensional mode, here denoted as Mode II-3D, extends to lower Reynolds values than the two-dimensional threshold of the non-rotating cylinder, and for a larger interval in α than the two-dimensional Mode II. The instability mechanism of Mode II-3D is of hyperbolic nature, see Pralits et al. (2013). Finally, note that the occurrence of two unstable modes has also been documented in the flow past rotating spheres (Citro et al. 2016;Fabre et al. 2017). However, the spatial structure of the direct and adjoint modes for our geometrical configuration is very different with respect to the case of the rotating sphere flow.