1. Introduction
The fluid–structure interaction problem of vortex-induced vibration (VIV) of a circular cylinder and, in particular, the case where it is elastically mounted and constrained to only oscillate transversely to the free stream, has been well studied in the literature, as can be seen from many comprehensive reviews: Sarpkaya (Reference Sarpkaya1979), Bearman (Reference Bearman1984), Parkinson (Reference Parkinson1989), Sarpkaya (Reference Sarpkaya2004), Williamson & Govardhan (Reference Williamson and Govardhan2004), Williamson & Govardhan (Reference Williamson and Govardhan2008), Bearman (Reference Bearman2011), Wu, Ge & Hong (Reference Wu, Ge and Hong2012) and Soti et al. (Reference Soti, Zhao, Thompson, Sheridan and Bhardwaj2018), amongst others. Over the last few decades, there have been many two-dimensional (2-D) numerical investigations of VIV at low Reynolds numbers, which, of course, are only valid if the wake remains 2-D. Interestingly, other studies have indicated that the onset of wake three-dimensionality is delayed if the cylinder oscillates transverse to the free stream (Leontini, Thompson & Hourigan Reference Leontini, Thompson and Hourigan2007; Gioria et al. Reference Gioria, Jabardo, Carmo and Meneghini2009). Since the development of three-dimensional (3-D) flow is likely to strongly affect the VIV response, it is important to understand the Reynolds number limits of the 2-D VIV modelling.
For a stationary cylinder, the transition from a 2-D periodic von Kármán vortex street to a 3-D wake was observed by noticing the discontinuities of the Strouhal number vs Reynolds number curve (Williamson Reference Williamson1988). Two dominant transition modes can be observed in the 3-D wake, which are referred to as mode A and mode B. Stability analysis reveals that the onset of the first 3-D mode, mode A, occurs at 
${\textit {Re}}_{cr} \approx 190$. Mode A is a spanwise periodic 3-D flow structure with a spanwise wavelength, 
$\lambda$, of approximately four cylinder diameters (
$D$) (Barkley & Henderson Reference Barkley and Henderson1996). This is consistent with the experimental finding of Williamson (Reference Williamson1988), who observed a spanwise wavelength of between 
$3D$ and 
$4D$ in the saturated wake beyond transition. A second 3-D wake mode, mode B, predicted to occur at 
${\textit {Re}}_{cr}= 260$ from Floquet stability analysis (FSA), instead occurs at 
${\textit {Re}}_{cr} \approx 230$ in experiments because the mode A saturated wake substantially alters the otherwise 2-D base flow. It also has a spanwise periodic structure with a critical wavelength of 
$\lambda \approx 0.8 D$, again consistent with experimental measurements (Williamson Reference Williamson1988).
Stability analysis is concerned with the breaking of flow symmetries as flow parameters (in this case the Reynolds number) are changed. The 2-D Kármán vortex street behind the stationary cylinder possesses spanwise translational and reflectional symmetries. It also possesses a spatiotemporal symmetry, i.e. a half-period time-shift and spatial reflection about the wake centreline result in the same flow state. If the streamwise flow is directed along the 
$x$ axis and the base flow 
$(U,V$) is generated in the 
$(x,y)$ plane, this symmetry possessed by the 2-D base flow is defined by
\begin{equation} (U,V) (x,y,t) = (U,- V)(x,-y,t+T/2), \end{equation}
where 
$T$ is the base-flow period. Alternatively, expressed in term of spanwise vorticity, 
$\varOmega$, this translates to
\begin{equation} \varOmega (x,y,t) ={-}\varOmega(x,-y,t+T/2). \end{equation}
Barkley, Tuckerman & Golubitsky (Reference Barkley, Tuckerman and Golubitsky2000) observed that mode A breaks the translational symmetry of the 2-D flow (base flow) 
$\boldsymbol {U}(x,y,z,t) = \boldsymbol {U}(x,y,z+l,t)$ (where 
$l$ is an arbitrary real constant), but preserves the spatiotemporal symmetry. Mode B also breaks the translational symmetry, but unlike mode A, mode B also breaks the spatiotemporal symmetry.
Another instability mode, mode C, not found for a stationary cylinder, was (falsely) reported to occur in the wake of square section cylinder (Robichaux, Balachandar & Vanka Reference Robichaux, Balachandar and Vanka1999). However, it does occur for circular cylinder wakes perturbed by tripwires (Sheard, Thompson & Hourigan Reference Sheard, Thompson and Hourigan2003). Mode C has an intermediate wavelength and occurs when the 2-D time-dependent wake does not have the exact spatiotemporal symmetry typically seen in the Kármán vortex street of a circular cylinder. In general, mode C appears in a stationary cylinder wake if a symmetry-breaking mechanism is used, such as a tripwire, or, for example, by imposing transverse oscillation. Such a mode is seen in the FSA of a transversely oscillating cylinder (Leontini et al. Reference Leontini, Thompson and Hourigan2007).
Another 3-D mode shown to be theoretically possible in a stationary or oscillating cylinder wake is a quasi-periodic (QP) mode. This mode has a complex-conjugate pair of Floquet multipliers, allowing it to consist of travelling or standing waves (Blackburn & Lopez Reference Blackburn and Lopez2003). Importantly, this means that the perturbation mode and base-flow periods are not commensurate. FSA for a stationary circular cylinder wake shows that mode QP becomes unstable at 
${\textit {Re}}_{cr} \approx 377$. At this 
${\textit {Re}}$, the wake is already highly 3-D in experiments although showing identifiable dominant mode B spanwise structures but displaying spatiotemporal chaos. These unstable 3-D modes have an intermediate spanwise wavelength in the range 
$1.5D$ to 
$2D$. In fact, the falsely identified mode C instability observed for a square cylinder wake by Robichaux et al. (Reference Robichaux, Balachandar and Vanka1999) was actually a QP mode, but with a period that was very close to subharmonic.
The dominant physical mechanisms underlying mode A and mode B instabilities are different. Williamson (Reference Williamson1996) suggested that mode A instability was associated with an elliptic instability of the vortex cores, and that mode B instability was associated with an instability of the braid region (which includes the braid shear layer within the near-wake vortex formation region (Leweke & Williamson Reference Leweke and Williamson1998; Thompson, Leweke & Williamson Reference Thompson, Leweke and Williamson2001)). The physical mechanism of mode C appears to be a hybrid of modes A and B (Sheard, Thompson & Hourigan Reference Sheard, Thompson and Hourigan2005). However, in each of these cases, the instability should perhaps not be identified with a single physical mechanism; for example, there appears to be some evidence that mode B is also influenced by centrifugal instability of near-wake vortices while showing elliptic instability characteristics within the spanwise vortex cores (Ryan, Thompson & Hourigan Reference Ryan, Thompson and Hourigan2005).
Most of these referenced studies were undertaken for circular cylinders, either fixed or under forced oscillation. They focused on finding the dominant 3-D modes, the influence of Reynolds number, 
${\textit {Re}} = UD/\nu$, and amplitude of oscillation, 
$A^*=A/D$, on the observed modes, and the effect of loss of spatiotemporal symmetry of the base flow on mode selection. Here, 
$U$ is the uniform flow speed, 
$\nu$ is the kinematic viscosity and 
$A$ the fixed oscillation amplitude.
Typically, for 2-D numerical investigations of VIV, the Reynolds number is chosen to be low enough to ensure a pre-transition 2-D wake. However, potentially interesting VIV wake states can be suppressed if the Reynolds number is too low, and additionally the connection between 2-D modelling and typically higher-Reynolds-number experiments becomes increasingly decoupled. Thus, the question arises as to the Reynolds number limit to ensure the wake is guaranteed to remain 2-D, or alternatively when 3-D modelling must be undertaken. To date, it appears that the FSA of an elastically mounted circular cylinder at low mass ratios has received little attention in the literature. For the forced oscillation case, the variation of amplitude and possible flow state affect the critical 
${\textit {Re}}$ for different modes (Leontini et al. Reference Leontini, Thompson and Hourigan2007) and hence transition to a 3-D wake. The primary motivation for this study is to determine the 
${\textit {Re}}$ limits of 2-D VIV simulations, and in particular, whether changes to the standard von Kármán wake state strongly affect wake stability. In effect, by examining free vibration, this study extends the forced-oscillation transition study of Leontini et al. (Reference Leontini, Thompson and Hourigan2007) by examining the influence of oscillation frequency, and how deviations from sinusoidal oscillation, deviations from vibration periodicity and wake–frequency coupling affect wake transition. In this paper, we investigate the critical Reynolds number for the 2-D to 3-D transition in a VIV system through Floquet analysis based on high-fidelity numerical simulations. We also undertake some complementary 3-D direct simulations beyond 3-D transition to quantify the effect on the amplitude response and the wake state. In addition, the effect of mass ratio is considered.
The layout of this paper is as follows. In § 2, the theoretical background of FSA, used to determine the linear stability of periodic flow, is briefly provided, and the numerical approach to solve this coupled system is also briefly discussed. The simulation results are given in § 3 for VIV. In this section, first response curves for a circular cylinder wake under VIV are presented for 
${\textit {Re}}= 150, 200$ and 250 to show why 2-D simulations are inadequate for accurately modelling higher-Re VIV. Next, the instability of Kármán wake to 3-D perturbations at 
${\textit {Re}} = 260, 280$ and 300 are discussed, together with the instability modes. Furthermore, the 3-D transition is investigated in detail for 
${\textit {Re}} = 300$ as a function of reduced velocity, including comparing 2-D and 3-D simulation predictions. 3-D predictions at higher Reynolds numbers are also briefly discussed. Finally, the critical Reynolds number for elastically mounted VIV wakes is presented as a function of reduced velocity.
2. Methodology
2.1. Theoretical background
FSA is used to predict the critical Reynolds number at which a 2-D periodic flow will first become unstable to 3-D perturbations, similar to how linear stability analysis (LSA) is used for stationary flows. This approach has been adopted by researchers to establish the onset of 3-D instabilities for different bluff-body wakes. Barkley & Henderson (Reference Barkley and Henderson1996) examined a fixed cylinder wake to predict the critical Reynolds number, 
${\textit {Re}}_{cr}$, the corresponding perturbation mode structure and the wavelengths of modes A and B at onset. Previous instability analyses of other bluff-body wakes include: the square cylinder (Robichaux et al. Reference Robichaux, Balachandar and Vanka1999), elongated cylinders for different aspect ratios (Ryan et al. Reference Ryan, Thompson and Hourigan2005), toroids (Sheard et al. Reference Sheard, Thompson and Hourigan2003), transverse-oscillating cylinders (Leontini et al. Reference Leontini, Thompson and Hourigan2007), rotational oscillating cylinders (Jacono et al. Reference Jacono, Leontini, Thompson and Sheridan2010) and elliptical cylinders of different aspect ratios (Leontini, Jacono & Thompson Reference Leontini, Jacono and Thompson2015).
For mode A, the Floquet perturbation velocity field (
$\boldsymbol {{u}'} = ({u}',{v}',{w}'$)) shows the following spatiotemporal symmetry (Barkley & Henderson Reference Barkley and Henderson1996; Robichaux et al. Reference Robichaux, Balachandar and Vanka1999):
\begin{equation} \left.\begin{gathered} {u}'(x,y, z, t) = {u}'(x, - y, z, t+T/2), \\ {v}'(x,y, z, t) ={-} {v}'(x, - y, z, t+T/2),\\ {w}'(x,y, z, t) = {w}'(x, - y, z, t+T/2), \end{gathered}\right\} \end{equation}
which is the same symmetry present in the 2-D base-flow wake (with 
$w = 0$). For mode A, the velocity components, 
${u}'$ and 
${w}'$, have an even spatiotemporal symmetry, whereas the cross-stream velocity 
${v}'$ has an odd spatiotemporal symmetry. Combining these, the representation of the mode by its streamwise vorticity (
${\varOmega }_x'$) gives
\begin{equation} {\varOmega}_x'(x,y, z, t) ={-}{\varOmega}_x'(x, - y, z, t+T/2), \end{equation}which obeys an odd spatiotemporal symmetry.
The velocity components of mode B have a spatiotemporal symmetry that is just the opposite of mode A and can be written as
\begin{equation} \left.\begin{gathered} {u}'(x,y, z, t) ={-}{u}'(x, - y, z, t+T/2), \\ {v}'(x,y, z, t) = {v}'(x, - y, z, t+T/2),\\ {w}'(x,y, z, t) ={-}{w}'(x, - y, z, t+T/2). \end{gathered}\right\} \end{equation}
In terms of streamwise vorticity (
${\varOmega }_x'$) this displays an even spatiotemporal symmetry expressed as
\begin{equation} {\varOmega}_x' (x,y, z, t) = {\varOmega}_x'(x, - y, z, t+T/2). \end{equation}The subharmonic mode C has a Floquet multiplier that is real and negative, and shows the following temporal symmetry for the vorticity:
\begin{equation} {\varOmega}_x'(x,y, z, t) = {\varOmega}_x'(x, y, z, t+2T). \end{equation}Finally, the QP mode, because it is non-commensurate with the base-flow frequency, does not possess such a relationship between flow states at different times.
In terms of flow physics, during shedding, the newly forming braid shear layer lies in close proximity to the previously formed braid. Unlike mode A, for mode B the braid comprises an array of rolled-up streamwise vortex pairs that extend between the spanwise rollers. The disturbances due to the forming braid set the preferred locations and rotation directions of the new braid vortices, giving a specific symmetry that is different from mode A; mode B has an in-phase symmetry for the streamwise vortex pattern.
Blackburn & Lopez (Reference Blackburn and Lopez2003) observed the Floquet multiplier for (pure) mode C contained no imaginary component. In addition, the sign of the perturbation field swapped each period, i.e. the 3-D features alternate in sign from one shedding cycle to the next. It leads to a periodic state with a period of oscillation twice that of the fundamental one, i.e. 
$\textrm {period} = 2T$ (Sheard et al. Reference Sheard, Thompson and Hourigan2005). Superficially, the spatial structure of the mode C perturbation field is similar to that of modes A and B. For a torus, the near-wake of the mode C instability, within approximately three diameters of the ring, bears a strong resemblance to the mode B instability, with strong perturbation vorticity present in the braid region between base-flow vortices. Further downstream, the perturbation field resembles that of the mode A instability, with perturbation vorticity localised within the base-flow spanwise vortex cores.
2.2. Floquet stability analysis
As indicated previously, FSA is method to determine the linear stability of the periodic wake to 3-D perturbations. In this analysis, the growth of linear perturbations is evaluated with respect to a pre-calculated periodic base flow. We consider a periodic 2-D wake (base-flow) 
$\boldsymbol {U}(x,y,t)$ with period 
$T$, and analyse its stability with respect to an infinitesimal 3-D perturbation 
$\boldsymbol {u'}(x,y,z,t)$. The equations that governs the perturbation field are obtained by substituting 
$\boldsymbol {U}(x,y,t)+ \boldsymbol {u'}(x,y,z,t)$ for 
$\boldsymbol {u}(x,y,z,t)$; 
$P(x,y,t)+ p'(x,y,z,t)$ for 
$p$, in (2.11). The governing equations for perturbation growth are obtained by linearising to give
\begin{equation} \frac{\partial \boldsymbol{u'}}{\partial t} ={-} (\boldsymbol{U} \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u'} - (\boldsymbol{u'} \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{U} - \boldsymbol{\nabla} p' +\frac{1}{Re} (\nabla^2 \boldsymbol{u'}), \end{equation}and
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u'}=0. \end{equation} These are simplified further by writing the solution as a summation of Fourier terms in the spanwise direction, 
$z$. The individual Fourier terms decouple so that the analysis can be done for each Fourier mode (i.e. spanwise wavelength) individually. The detailed approach is described in Barkley & Henderson (Reference Barkley and Henderson1996).
The boundary conditions are shown in figure 1 and are described as follows. The fluid velocity is prescribed at the inlet, top and bottom boundaries, and is given by 
$(U,V) = U_\infty \textrm {i} - u_{cyl} \textrm {j}$. The perturbation velocity components at these far-field boundaries are 
$u' = v' = w' = 0$, where 
$u'$, 
$v'$ and 
$w'$ are the 
$x$, 
$y$ and 
$z$ perturbation velocity components, respectively. At the surface of the cylinder, the no-slip condition is imposed, i.e. 
$(U,V) = 0$ and 
$(u',v',w') = 0$. At the outlet, the boundary condition is the same as the base-flow condition, i.e. the normal velocity gradient is set to zero for both the base and perturbation velocity fields, and the pressure is fixed. At no-slip boundaries and at the far-field boundaries, higher-order boundary conditions are used for the normal pressure gradient for both the base and perturbation pressure, as described in Karniadakis, Israeli & Orszag (Reference Karniadakis, Israeli and Orszag1991).
Schematic of the computational domain and the boundary conditions (BC) for Floquet analysis used in the present study.
2.3. Perturbation kinetic energy
The transport equation for the perturbation kinetic energy can be derived from the stability equations above by taking the dot product with the perturbation velocity. In tensor notation, which leads to easier manipulation, this gives
\begin{equation} u'_i\frac{\partial u'_i}{\partial t} ={-} u'_i U_j \frac{\partial u'_i}{\partial x_j} - u'_i u'_j \frac{\partial U_i}{\partial x_j} - u'_i\frac{\partial p'}{\partial x_i} + \frac{1}{Re} u'_i\frac{\partial_{\phantom{i}}}{\partial x_j} \left(\frac{\partial u'_i}{\partial x_j}\right). \end{equation}This can be rewritten as
\begin{equation} \frac{\partial k'}{\partial t} + U_j \frac{\partial k'}{\partial x_j} ={-} u'_i u'_j \frac{\partial U_i}{\partial x_j} -\frac{\partial u_i'p'}{\partial x_i} + \frac{1}{Re} \left(\frac{\partial^2 k'}{\partial x_j^2}\right), \end{equation}
where 
$k' = \frac {1}{2}u'_iu'_i$ is the perturbation kinetic energy. The left-hand side is just the time derivative of the kinetic energy following an advected fluid parcel whereas the last term represents viscous diffusion of kinetic energy. Of the remaining two terms, the first is the only one referencing the base flow and, hence, represents transfer of energy to and from the base flow, whereas the second represents work done by the perturbation pressure on the perturbation velocity field. Thus, examining the first term on the right-hand side provides insight into the generation of perturbation kinetic energy. Further, by expanding this term into Fourier components in the spanwise direction, i.e. 
$\boldsymbol {u'} = (\hat {u}(x,y,t)\cos (kz), \hat {v}(x,y,t)\cos (kz), \hat {w}(x,y,t)\sin (kz)),\ p' = \hat {p}\cos (kz)$, and integrating over a spanwise wavelength, gives the spanwise average approximation per unit span as
\begin{equation} \dot{K}' = \frac{1}{\lambda}\int_0^\lambda \left({-}u'_i u'_j \frac{\partial U_i}{\partial x_j}\right){\rm d} z ={-}\frac{1}{2}\left(\hat{u}\hat{u} \frac{\partial U}{\partial x} + \hat{v}\hat{v} \frac{\partial V}{\partial y} + \hat{u}\hat{v} \left(\frac{\partial U}{\partial y} + \frac{\partial V}{\partial x}\right)\right), \end{equation}
where 
$\lambda = 2{\rm \pi} /k$ is the wavelength corresponding to wavenumber 
$k$. This term can be used to examine which parts of the wake contribute to a particular instability mode.
2.4. The spectral-element method
Fluid flow is modelled in the moving reference frame attached to the cylinder so that remeshing is not required. The governing equations are then the non-dimensional Navier–Stokes equations in an accelerated frame of reference (Thompson, Hourigan & Sheridan Reference Thompson, Hourigan and Sheridan1996; Thompson et al. Reference Thompson, Hourigan, Cheung and Leweke2006)
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} ={-} (\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u} - \boldsymbol{\nabla} p + \frac{1}{Re} (\nabla^2 \boldsymbol{u}) - \boldsymbol{a}, \end{equation}and
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{equation}
where 
$\boldsymbol {u}$ and 
$p$ are the dimensionless fluid velocity and kinematic pressure, respectively, and 
$\boldsymbol {a}$ is the acceleration of the reference frame attached to the cylinder, 
$\boldsymbol {a} = {\textrm {d}\boldsymbol {u_{cyl}}}/{\textrm {d} t}$, where 
$\boldsymbol {u_{cyl}}$ is the cylinder velocity vector.
These equations are discretised in space using the spectral-element method (Thompson et al. Reference Thompson, Hourigan and Sheridan1996, Reference Thompson, Hourigan, Cheung and Leweke2006; Karniadakis & Sherwin Reference Karniadakis and Sherwin2013). This is essentially a high-order finite-element approach, where the fluid domain is broken up into quadrilateral elements or cells, and within each cell the function is represented by a tensor product of high-order Lagrange polynomial interpolants based on the node points of Gauss–Legendre–Lobatto quadrature. This allows efficient evaluation of the integrals resulting from the application of the weighted residual method to the governing fluid equations. Importantly, the method converges exponentially (or spectrally) as the number of nodes/elements is increased. The method is described in detail in the above references.
In terms of time discretisation, a second-order three-step time-splitting method is used (Karniadakis et al. Reference Karniadakis, Israeli and Orszag1991). This treats the effects of advection, mass conservation (and pressure) and diffusion separately. The advection step is treated explicitly using a third-order Adams–Bashforth approximation, whereas the linear pressure and diffusion substeps are treated implicitly resulting in linear matrix problems. This means that these substeps of the time integration are reduced to matrix–vector multiplies after the computing the lower–upper decomposition of the matrices in a pre-processing step. Further details can be found in Thompson et al. (Reference Thompson, Hourigan, Cheung and Leweke2006). This computational approach is used for both the fluid (2.11), (2.12) and the stability analysis (2.6) equations.
The dimensionless structural differential equation, which provides the acceleration of the cylinder, is given by
\begin{equation} \ddot{Y}^* + 4 {\rm \pi}\zeta f_n^* + 4 {\rm \pi}^2 {f_{n}^*}^2 {Y}^* = \frac{2}{{\rm \pi} m^*} {C_L}. \end{equation}
Here, 
$\ddot {Y}^*=\ddot {y}D/U^2$, 
$\dot {Y}^*=\dot {y}/U$, 
$Y^* = y/D$ and 
$C_L = F_y/(\frac {1}{2}\rho _f U^2 D L)$ correspond to the non-dimensional acceleration, velocity, displacement and lift force, respectively, with 
$\rho _f$ the fluid density and 
$\zeta$ the damping ratio. Here 
$L$ represents the spanwise length of the cylinder and is taken as 1. In addition, 
$f_n^* = Df_n/U$ is the non-dimensional natural frequency of the system in a vacuum ( 
$f_n = ({1}/{2{\rm \pi} })\sqrt {k/m}$), and 
$m^* = \rho _c/\rho _f$ is mass ratio, with 
$\rho _c$ the cylinder density. In general, the oscillation response is recorded as a function of the reduced velocity, 
$U_r = 1/f_n^*$, although often the theoretical system oscillation frequency in the working fluid, which takes account of the added mass contribution, is used to define an alternative reduced velocity instead, i.e. 
$U^* = 1/f_w^*$, where 
$f_w^* = Df_w/U$, with 
$f_w = ({1}/{2{\rm \pi} })\sqrt {k/(m+m_A)}$. Here, 
$m_A = C_A\rho _c {\rm \pi}D^2/4$ is the added mass with 
$C_A = 1$ for a circular cylinder. There are also two other important frequencies that need to be considered. These are 
$f_v$, the shedding frequency of a stationary cylinder, and 
$f$ or 
$f_y$, the actual measured frequency of the cylinder response when the oscillations are periodic. The latter is, of course, different but related to the other frequencies.
These equations are discretised in time following the approach described in Rajamuni, Thompson & Hourigan (Reference Rajamuni, Thompson and Hourigan2018), and are integrated simultaneously in time with the fluid equations to provide the frame acceleration as a function of time. Typically, we take the damping ratio 
$\zeta$ to be zero.
For the stability analysis undertaken for this study, the fluid (2.11), (2.12) and structural (2.13) equations are integrated in time until a periodic state is reached. At that stage, the velocity field is recorded at discrete integrals (typically 50–100 snapshots) over a complete period. These fields are used to provide the velocity field at each timestep for the linear stability equations (2.6) using quadratic interpolation based on the three closest saved times.
2.5. Steady to unsteady transition
The initial Hopf bifurcation (i.e. steady to unsteady transition) for the transition to 2-D periodic flow is also examined for the elastically mounted case. Previously, Mittal & Singh (Reference Mittal and Singh2005) found that a possible resonance between the shedding and system frequencies could lower the transition Reynolds number considerably from the fixed cylinder case. In this paper, we reexamine this transition especially focusing on the effect of mass ratio.
In terms of stability theory, the transition is from a stationary cylinder to coupled motion, which means we can expand the fluid and body variables about the steady base state as
\begin{equation} \boldsymbol{u} \rightarrow \boldsymbol{U}(x,y) + \boldsymbol{u}'(x,y,t),\quad p \rightarrow P(x,y) + p'(x,y,t),\quad Y^*{} \rightarrow Y^*{}'(t). \end{equation}Thus, after linearising and subtracting the equations of the base flow, the governing stability equations for this case are
\begin{gather} \frac{\partial \boldsymbol{u'}}{\partial t} ={-} (\boldsymbol{U} \boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u'} - (\boldsymbol{u'} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{U} - \boldsymbol{\nabla} p' + \frac{1}{Re} (\nabla^2 \boldsymbol{u'})- \frac{{\rm d}^2Y^*{}'}{{\rm d} t^2}, \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u'}=0, \end{gather}together with
\begin{equation} {\ddot{Y}}^*{}' + 4 {\rm \pi}\zeta f_n^* {\dot{Y}}^*{}'+ 4 {\rm \pi}^2 {f_{n}^*}^2 {Y^*}{}' = \frac{2}{{\rm \pi} m^*} {C_L'}. \end{equation}
Only the zero damping case is considered here, i.e. 
$\zeta = 0$.
The dominant eigenmodes can be extracted using the same approach as for the Floquet analysis, i.e. by integrating forward in time from an initial random velocity field until only the most dominant modes remain and then using Arnoldi decomposition on a set of regularly saved snapshots to extract the most dominant modes.
2.6. 3-D simulations
Subsequent to undertaking the linear and FSA, we also undertook some fully 3-D flow modelling to examine the cylinder response and the wake state post-3-D transition. This was done by extending the 2-D spectral-element code using a Fourier expansion in the spanwise direction as described in Karniadakis & Triantafyllou (Reference Karniadakis and Triantafyllou1992), Thompson et al. (Reference Thompson, Hourigan and Sheridan1996) and Gupta et al. (Reference Gupta, Zhao, Sharma, Amit, Hourigan and Thompson2023). Typically, 128 Fourier planes were used to cover a span length of 
$12D$. The latter was chosen as approximately three times the wavelength of mode A. In addition, this gives a spanwise (
$z$) spatial resolution of 
${\sim }0.1D$, which is much smaller than the smallest instability wavelength (mode B, wavelength 
${\sim }D$). It was also necessary to increase the cross-plane (
$x$–
$y$) resolution for the higher Reynolds number simulations: relative to the 2-D simulations, the intra-element nodal resolution was increased from 
$8 \times 8$ to 
$10 \times 10$ for the 
$Re=500$ simulations.
Most simulations at 
$Re=1000$ were undertaken by switching on spectral-vanishing viscosity (SVV) in the spectral/spectral-element code (Kirby & Sherwin Reference Kirby and Sherwin2006). This increases the viscosity experienced by shorter length-scale features, which stabilises the code at higher Reynolds numbers. An important point is that the range of length scales (wavelengths) affected decrease as the resolution is increased. Sometimes this is referred to as implicit large-eddy simulation (LES). The implementation and details follow its application and recommended use in the open-source software Nektar++ (described in Kirby & Sherwin Reference Kirby and Sherwin2006), which has been successfully used to simulate flow past a circular cylinder at 
$Re=3900$ (Jiang & Cheng Reference Jiang and Cheng2021), noting an excellent match with PIV velocity fields obtained by Parnaudeau et al. (Reference Parnaudeau, Carlier, Heitz and Lamballais2008) in that case. The implementation of SVV in current code has been verified using the same benchmark 
$Re=3900$ non-oscillating cylinder flow, giving statistically identical predictions to the published predictions of Jiang & Cheng (Reference Jiang and Cheng2021).
A resolution study was undertaken to ensure that SVV was minimally influencing statistical maximum response amplitude predictions. Simulations were run for 
$m^*=2.456$ using meshes with 518 spectral elements and 128 Fourier planes with both 
$8 \times 8$ and 
$9 \times 9$ noded elements over the entire 
$U^*$ range. The mode cutoff recommendations provided by Nektar++ documentation, i.e. the lower 75 % of spectral-element (polynomial) modes remained undamped and the SVV diffusion coefficient was set at 0.1 (Nektar++ User Guide, also see Jiang & Cheng (Reference Jiang and Cheng2021) for the 
$Re=3900$ cylinder flow predictions and the effect of the choice of SVV parameters). The timestep for these cases was set at 
${\rm \Delta} t = 0.005 D/U_\infty$. The maximum variation in the amplitude prediction between these two cases across the entire 
$U_r$ range was less than 2 %, which is consistent with the statistical sampling error given finite integration times. In addition, higher-resolution simulations were undertaken for a reduced velocity corresponding to an upper-branch response (
$U^*=5.0$ or 
$U_r=4.2$) using a mesh with 2072 
$7 \times 7$ elements and 128 or 256 Fourier planes. The maximum amplitude response predictions between these four cases were in the range 0.766–0.770, i.e. a variation of approximately 0.5 %. Finally, simulations were undertaken using the 2072 
$7 \times 7$ 128 Fourier plane mesh without SVV stabilisation. These were stable for most 
$U^*$ cases outside the upper branch although a smaller timestep was required for stability. The predicted maximum amplitude responses for the stable simulations were in statistical agreement with the SVV stabilised predictions. In the upper branch, presumably the spatial resolution was not quite adequate to resolve the transfer of energy down to the smallest scales, noting this is a turbulent flow. Although higher resolution could have been used to fully capture scales down to the Kolmorogov scale, this would have been computationally prohibitive, but also noting that this higher-Reynolds-number case is not a major focus of the paper. Of course, the response curve for 
$Re=1000$ shown in figure 21 is in good agreement with the predictions of Zhao et al. (Reference Zhao, Cheng, An and Lu2014) obtained using a stabilised finite-element code for almost the same mass ratio.
2.7. 2-D resolution studies
The implementation of the spectral-element approach used in the current study has been validated extensively in previous studies, see, e.g., Hourigan, Thompson & Tan (Reference Hourigan, Thompson and Tan2001), Sheard et al. (Reference Sheard, Thompson and Hourigan2003), Rao et al. (Reference Rao, Radi, Leontini, Thompson, Sheridan and Hourigan2015) and Leontini et al. (Reference Leontini, Thompson and Hourigan2007) and references therein. Resolution studies for base flow and stability calculations are considered separately. For the base-flow simulations, the domain and element distribution is based on the previous study of VIV of a circular cylinder, where resolution and other validation studies were reported (Mishra et al. Reference Mishra, Soti, Bhardwaj, Kulkarni and Thompson2020). For the stability calculations, care has been taken during the construction of the mesh to ensure adequate resolution, as verified through the base-flow study. However, perturbation fields can have smaller length scales than the base flow, especially as the spanwise wavelength is reduced, so we investigated the convergence with respect to the polynomial order of the tensor-product intra-element interpolating polynomials. Three sets of element resolutions, 
$6\times 6$, 
$8\times 8$ and 
$10\times 10$ (nodes/element) for a parameter set that leads to near-maximal displacement: 
${\textit {Re}} = 300$, 
$\zeta = 0.0$, 
$m^* = 2.546$ and natural frequency 
$f_n^* = 0.2$ have been used. Figure 2 indicates that increasing the polynomial order from 7 (
$8 \times 8$) to 9 (
$10 \times 10$) changed the value of the Floquet multiplier less than 0.3 % over the range of spanwise wavelengths investigated, providing confidence that the results presented in this paper are adequately converged.
Floquet multiplier modulus, 
$|\mu |$, vs spanwise wavelength, 
$\lambda /D$, of Floquet modes as a function of number of nodes per spectral element, for a cylinder undergoing VIV at 
${\textit {Re}} = 300$, 
$\zeta = 0.0$, 
$f_n^* = 0.2$ and 
$m^* = 2.546$: (a) influence of element-order on the mode B Floquet multiplier for the mesh used for the bulk of the simulations; (b) influence of domain size.
2.8. Influence of computational domain size
We also investigated the influence of the locations of the outer computational boundaries relative to the centre of the cylinder. Two different domain sizes were examined. The first domain had dimensions of 
$15D$ upstream and 
$25D$ downstream from the cylinder centre, as depicted in figure 1, and was discretised with 518 macro-elements. The second domain extended 
$50D$ upstream and 
$50D$ downstream from the cylinder centre, resulting in a mesh consisting of 1134 elements. Figure 3 displays the amplitude of oscillation (
$A^*$) and the root mean square of the lift coefficient (
$C_{L,rms}$) for 
${\textit {Re}} = 150$, 
$m^* = 2.546$ and 
$\zeta = 0$. This indicates only minor changes to the response curve over the entire reduced velocity range, with only slightly larger variation in the lift coefficient, which is more sensitive. In addition, the effect on the Floquet multiplier for mode B is shown in figure 2(b), which shows a difference of less than 0.2 % near the peak. Thus, the smaller domain size was selected for the remainder of the simulations because of the reduced computational resource requirement for this large parameter study.
Effect of domain size of amplitude 
$A^*$ and root mean square of the lift coefficient (
$C_{L,rms}$) for 
${\textit {Re}} = 150$, 
$m^* = 2.546$ and 
$\zeta = 0$.
2.9. Validation
We have undertaken a validation study to ensure that the results produced are in agreement with published results for standard test cases. Comparisons of base-flow VIV predictions for the current set-up against literature results were reported in Mishra et al. (Reference Mishra, Bhardwaj, Kulkarni and Thompson2021). We conducted simulations at 
${\textit {Re}} = 150$ and 
$m^* = 2.546$ across a range of reduced velocities. Figure 4 illustrates the variations of the maximum displacement of the cylinder, mean lift coefficient and root mean square of drag and lift coefficients with reduced velocity. The comparison reveals that our results closely align with those presented in Bao et al. (Reference Bao, Huang, Zhou, Tu and Han2012) and Zhao (Reference Zhao2013). Any minor differences observed at specific reduced velocities may be attributed to variations in blockage ratios considered.
Comparisons of the predicted amplitude response and mean and fluctuating force coefficients against simulation predictions from the literature. Here, 
$\zeta = 0$, 
$m^* = 2.546$ and 
${\textit {Re}} = 150$. See the text for further details.
Validation of stability predictions under VIV is through a comparison against forced oscillation predictions of Leontini et al. (Reference Leontini, Thompson and Hourigan2007), while noting that Gioria et al. (Reference Gioria, Jabardo, Carmo and Meneghini2009) has produced similar predictions for a different forcing frequency ratio. In particular, Leontini et al. (Reference Leontini, Thompson and Hourigan2007) presented a transition map showing where different Floquet modes become unstable as a function of forced oscillation amplitude and Reynolds number for a fixed forcing frequency of 
$f^* = f(D/U) = 0.2$, noting that this is very close to the vortex shedding frequency at 
$Re=200$ and for higher Reynolds numbers after the flow has become 3-D. For the current problem, the oscillation amplitude is a function of the reduced velocity, so by running simulations over a range of Reynolds numbers, the transition Reynolds number can be determined and the oscillation amplitude noted for a particular system frequency. For 
$f^* = 0.2$, the transition Reynolds number can then be compared with the published result of Leontini et al. (Reference Leontini, Thompson and Hourigan2007), under the assumption that the displacement variation is close to sinusoidal in the freely oscillating case.
The validation has been done by extracting the critical Reynolds number corresponding to the observed amplitude for each transition mode from the forced oscillation study of Leontini et al. (Reference Leontini, Thompson and Hourigan2007, figure 5). Those predictions show that for an undamped forced oscillation amplitude of 
$A^* =0.5$ for 
$f^* = 0.2$, mode A becomes unstable at 
${\textit {Re}}_{c,A} \sim 305$ and mode B at 
${\textit {Re}}_{c,B} \simeq 277$. Thus, simulations were run for 
${\textit {Re}} = 277$ and 305 for 
$\zeta = 0$ (undamped system) and 
$m^* = 2.546$, corresponding to the standard VIV spring–damper system. Since 
$f_n^* = 0.2$ (i.e. 
$U_r = 5$) does not give a response oscillation frequency of 
$f^* = 0.2$, it was first necessary to determine the 
$U_r$ corresponding to 
$f^*=0.2$ for each Reynolds number. For both Reynolds numbers, the saturation amplitude is close to 
$A \simeq 0.5$. Figure 5 shows the power spectrum of the displacement signal for the 
$Re=305$ case, indicating that the response is very close to a pure sine wave with only minor higher-order frequency content at the odd harmonics. Figure 6 shows that at 
${\textit {Re}} = 277$, the maximum amplitude multiplier is close to unity for the solution branch corresponding to mode B, whereas for 
${\textit {Re}}=305$, the maximum amplitude of the mode A branch also reaches a value close to unity. Thus, both of these predictions are in close accord with those from the forced oscillation study, with the same critical Reynolds numbers of 
${\textit {Re}}_{cr} = 277$ and 305. This provides validation of the implementation of FSA for the elastically mounted spring–damper VIV system. Further validation for the free-oscillation case was provided by examining the initial growth phase of 3-D perturbations for 3-D simulations at post-critical Reynolds numbers, which showed that the measured growth rates in the linear growth regime matched the predictions of FSA. This comparison is not shown because of space considerations.
Power spectrum of the displacement signal for 
$f^* = 0.2$, 
$m^* = 2.546$ at 
$Re=305$. This indicates that the displacement signal for this reduced velocity (
$U_r = 5.194$) the oscillatory response is very close to a pure sinusoid, in accordance with the assumption required for a direct comparison of free and forced vibration cases. See the text for further details.
Floquet multiplier vs spanwise wavelength at 
${\textit {Re}} = 277$ and 305, for 
$\zeta = 0$, 
$f^* = 0.20$, 
$m^* = 2.546$, corresponding to a standard spring–damper VIV system. Mode B becomes critical first at 
${\textit {Re}} \simeq 277$ and then mode A for 
${\textit {Re}} = 305$. In this figure, 
$A^*$ is the amplitude of oscillation. A comparison with Leontini et al. (Reference Leontini, Thompson and Hourigan2007) shows that for a undamped forced oscillation amplitude of 
$A^* =0.5$ and oscillation frequency = 0.2, and for purely sinusoidal oscillation, mode A becomes unstable at 
${\textit {Re}}_{c,A} \sim 305$ and mode B at 
${\textit {Re}}_{c,B} \simeq 277$.
3. Results
In this section, the stability of the wake to 3-D transition for the VIV for a circular cylinder under elastic support is analysed. There are three independent parameters in this study: mass ratio (
$m^*$), reduced velocity (
$U_r$) and Reynolds number (
${\textit {Re}}$). Recall that the reduced velocity is defined as 
$U_r = U / (f_n D)$, where 
$U$ is the free-stream velocity, 
$f_n$ is the natural frequency of the structure in a vacuum. The limit of two-dimensionality for the oscillating cylinder is determined by varying three independent parameters, 
${\textit {Re}}$, 
$U_r$, i.e. the non-dimensional oscillation period, and 
$m^*$.
We initially analyse an undamped elastically-mounted cylinder for 
$m^* = 2.546$ at 
$f_n^* = 0.2$, noting that the natural frequency is close to the natural vortex-shedding frequency of a fixed cylinder. The occurrence of modes A, B and QP is examined as 
${\textit {Re}}$ is increased. We determine the critical values of 
${\textit {Re}}$ and non-dimensional spanwise wavelength, 
$\lambda /D$, for which the 2-D wake becomes linearly unstable for the periodic oscillation of a cylinder under VIV.
3.1. Critical 
${\textit {Re}}$ for the onset of vortex shedding
Before investigating the loss of stability of the Kármán wake to 3-D perturbations, it is useful to determine the underlying periodic base-flow states. Although this dependence of the transition on reduced velocity has been examined previous by Mittal & Singh (Reference Mittal and Singh2005) for the two-degree-of-freedom case, they only considered mass ratios 
$m^* \gtrsim 5$. For the current simulations, this range is extended to lower mass ratios where the effect on transition is amplified.
For a stationary circular cylinder, the wake undergoes a transition from a stable to periodic Kármán wake at 
${\textit {Re}} \sim 47$ (e.g. Dusek, Gal & Fraunie Reference Dusek, Gal and Fraunie1994; Williamson Reference Williamson1996; Le Gal, Nadim & Thompson Reference Le Gal, Nadim and Thompson2001). This instability leads to unsteady lift and drag forces that cause an elastically mounted cylinder to undergo VIV. The onset of this periodic shedding occurs at different Reynolds numbers for a cylinder undergoing VIV. The value of this critical 
${\textit {Re}}$ depends on the reduced velocity.
Based on the stability analysis of § 2.5, figure 7 provides growth rate contours for mass ratios 
$m^*=1$, 
$2.546$ and 
$10$. The zero growth rate curves are marked, which indicate the transition Reynolds number as a function of reduced velocity.
Growth rate contour maps for the steady to unsteady transition of an elastically mounted cylinder showing the effect of mass ratio on transition as a function of reduced velocity. The rows correspond to (a) 
$m^* = 1$, (b) 
$m^*=2.546$ and (c) 
$m^*=10$. The zero growth contour is marked indicating the variation of the transition Reynolds number as 
$U^*$ is varied.
As 
$U_r$ is increased from very small values, the critical Reynolds number, 
$Re_{co}$, initially rises from a value consistent with the stationary cylinder transition Reynolds number of 
$Re_{co} \simeq 46$ to reach a peak, prior to rapidly dropping after a different instability mode becomes dominant (corresponding to the system frequency rather than the unsaturated shedding frequency). This shift corresponds to the kinks in the contour lines of figure 7. For the lightest mass ratio considered, 
$m^*=1$, the critical Reynolds number is increased to approximately 
$Re_{co}=51$ at 
$U_r=4.6$. For larger mass ratios, the effect is less; for 
$m^*=10$, the increase is only to 
$Re_{co} \simeq 47$ at 
$U_r=5.7$ before again dropping rapidly. In all cases, beyond the initial peak, the critical Reynolds number reaches a value of slightly less than 20, noting that the minimum occurs at different reduced velocities for the different mass ratios. This result is consistent with the findings of Mittal & Singh (Reference Mittal and Singh2005), who also found a minimum critical Reynolds number of 
${\sim }20$.
3.2. Variation of the mode growth close to 3-D transition
In this section, for the system undergoing VIV, loss of stability of the Kármán wake to 3-D perturbations is analysed for a mass ratio 
$m^* = 2.546$, as a function of Reynolds number. The Reynolds number range chosen covers the range where analogues of both modes A and B become unstable. Figure 8 shows the magnitude of the predicted Floquet multiplier, 
$|\mu |$, vs wavelength, 
$\lambda /D$, for undamped free vibration. The natural frequency is fixed at 
$f_n^* = 0.2$, close to shedding frequency of the fixed cylinder; this corresponds to 
$U_r = 5.0$. Figure 8 shows that for 
${\textit {Re}}= 300$ and 320, the flow is unstable to 3-D perturbations, whereas at 
${\textit {Re}} = 280$ and 
$260$ the 2-D wake flow remains stable (
$|\mu | < 1$), with no Floquet modes showing positive growth. The amplitude of oscillation increases slightly with decreasing 
${\textit {Re}}$, with 
$A^* = 0.49, 50$, 
$0.50$ and 
$0.51$, for 
${\textit {Re}}= 320$, 
$300$, 280 and 260, respectively. The figure shows that for 
${\textit {Re}} \geq 280$, two clear peaks can be observed although both remain stable at 
$Re=280$, whereas for 
${\textit {Re}} = 260$ only the mode B peak has become visible. The modes corresponding to the two branches are analogous to modes A and B for the stationary cylinder. We observe that mode B becomes critical first then mode A; a similar result was found for 
$A^* > 0.3$ under forced cylinder oscillation by Leontini et al. (Reference Leontini, Thompson and Hourigan2007).
Floquet multiplier vs wavelength at 
${\textit {Re}} = 320$, 
$300$, 
$280$ and 
$260$; 
$\zeta = 0$ (undamped system), 
$f_n^* = 0.20$ and 
$m^* = 2.546$.
In the case of 
${\textit {Re}}= 300$, only mode B is unstable, with the maximum growth rate occurring at 
$\lambda \approx 0.9 D$. The multiplier branch corresponding to Mode QP, which has a complex Floquet multiplier, remains stable but is dominant for 
$1.4D < \lambda < 1.9 D$. At 
${\textit {Re}} = 300$, mode A is also stable with a maximum 
$|\mu | \approx 0.9$ occurring at 
$\lambda \approx 2.5D$. In fact, for 
$\lambda > 3.7 D$, the dominant mode has a complex Floquet multiplier corresponding to a QP mode. Since it is observed at a longer wavelength than mode A, to make a clear distinction from the previously observed QP mode, it is designated as mode QPL, named for long-wavelength QP. For 
${\textit {Re}} = 280$, the flow is near to unstable to mode B, again with a wavelength centred around 
$\lambda \approx 0.9D$, whereas the branch corresponding to mode QP remains stable and exists for 
$\lambda > 1.4 D$ although it is not shown explicitly here. Thus, modes A and B become critical at lower Reynolds numbers than those of modes QP and QPL. This latter behaviour is similar as that reported by Blackburn, Marques & Lopez (Reference Blackburn, Marques and Lopez2005) for a stationary cylinder, who estimated the critical 
${\textit {Re}}$ for mode QP to be about 377 with 
$\lambda _{cr} \approx 1.8 D$. Mode QP has no local maximum over this Reynolds number range, unlike modes A and B; however, we have not examined the stability for 
${\textit {Re}} > 300$.
Figure 9 shows the spanwise and streamwise perturbation vorticity fields at 
${\textit {Re}} = 300$ in order to indicate the similarities and differences in the structures of modes A and B. The solid lines in the figure outline the base-flow vorticity with two single vortices shed per oscillation cycle (referred to as a 
$2S$ wake mode). This mode occurs at low amplitudes of oscillation, in particular for 
$A^* \lesssim 0.5$ at 
${\textit {Re}} =300$ (Leontini et al. Reference Leontini, Stewart, Thompson and Hourigan2006a). Mode B in figure 9(a,b) is the first occurring unstable mode showing a maximum Floquet multiplier at 
$\lambda = 0.9D$. Mode B obeys the spatiotemporal symmetry of the base flow about the wake centre line, with spanwise perturbation vorticity repeated (on reflection through the centreline) every half base period. The successive spanwise vortex structures from each side of the wake have the same sign. The still stable mode A correspondingly shows a maximum 
$|\mu |$ at 
$\lambda = 2.5D$. This mode is shown in figure 9(c,d). The primary vortex at the rear of the cylinder has both positive and negative regions of perturbation vorticity; further downstream, the perturbation becomes concentrated in the primary vortex cores. Mode A has spatiotemporal symmetry opposite to that of mode B. The perturbation field structures for modes A and B are found to be similar to those of modes A and B for a fixed circular cylinder (Barkley & Henderson Reference Barkley and Henderson1996), although the oscillation certainly has a strong effect on the underlying base flow. In particular, the near-wake base-flow vortex structures are quite different from those for a stationary cylinder in that they are stretched and elongated as the cylinder moves up and down rather than more elliptical in character for the stationary cylinder. Since mode A has been at least partially attributed to elliptic instability (Leweke & Williamson Reference Leweke and Williamson1998; Thompson et al. Reference Thompson, Leweke and Williamson2001), this may help to explain why mode A remains stable until higher Reynolds numbers under oscillation. Since mode B is likely to depend more on the shear of the braid regions between Kármán wake vortices (Leweke & Williamson Reference Leweke and Williamson1998), which of course persist in the oscillatory case, this is consistent with the smaller effect on the onset of mode B. However, this question is further considered in the following through an examination of the perturbation kinetic energy generation.
Spanwise perturbation vorticity depicting (a,b) mode B and (c,d) mode A at 
${\textit {Re}} = 300$, 
$\zeta = 0$, 
$f_n^* = 0.20$ and 
$m^* = 2.546$. Red (blue) represents positive (negative) spanwise vorticity. The solid lines outline base-flow vorticity.
The effect of mass ratio on 3-D stability under VIV at 
${\textit {Re}}= 300$ is examined for mass ratios 
$m^* = 1$, 2.546 and 10 in figure 10. The reduced velocity is adjusted to result in 
$f^* = 0.20$ in each case. The response amplitude is close to 0.49 for each mass ratio. Indeed, the dominant Floquet multiplier curves are virtually identical as the mass ratio is varied by a factor of 10, again showing a collapse against the actual system response frequency (rather than 
$U_r$ or 
$U^*$), and consistent with figure 6.
Floquet multiplier vs wavelength for (a) 
$m^* = 10$, 2.546 and 1, at 
$Re=300$ for 
$f^* = 0.2$.
3.3. Relationship of the 3-D instability modes to stationary cylinder modes A, B and QP
Figure 11 plots the spatial distributions of the perturbation kinetic energy (
$k'$) and generation rate term (
$\dot {K}'$) given by (2.10) for the different instability modes. This latter term (
$\dot {K}'$) provides the rate of transfer of perturbation kinetic energy to or from the base flow, thus providing insight into which parts of the wake contribute to the development of the instability modes. This analysis was undertaken for modes A, B and QP for the stationary circular cylinder and the analogous A, B and QP instability modes for a cylinder undergoing VIV.
Perturbation kinetic energy, 
$k'$, (a,c,e,g,i,k), and kinetic energy generation rate, 
$\dot {K}'$, (b,d, f,h, j,l), from (2.10), which indicates the generation of perturbation kinetic energy from the base flow: (a–d) mode A comparison for the stationary and VIV cases; (e–h) same for mode B, bottom two rows; (i–l) same for mode QP. See the text for further details.
Each pair of rows of this figure shows side-by-side comparisons of the dominant stationary-cylinder and VIV oscillating-cylinder instability modes close to onset. The left column shows perturbation kinetic energy and the right column the perturbation energy generation rate. The top, middle and bottom pair of rows show comparisons for modes A, B and QP, respectively. Clearly, for mode A of the stationary cylinder, energy input into the mode occurs primarily in the wake vortex cores, consistent with the interpretation of elliptic instability assisting the maintenance of the instability (Leweke & Williamson Reference Leweke and Williamson1998; Thompson et al. Reference Thompson, Leweke and Williamson2001). On the other hand, for mode B perturbation energy transport from the base flow is primarily in the sheared region between the near-wake vortices and extending into the braids between the row of vortices further downstream. This is consistent with the interpretation that the instability is supported mostly by an instability of the braids (Leweke & Williamson Reference Leweke and Williamson1998). Mode QP appears to be supported through a more equal combination of elliptic and hyperbolic instabilities (Sheard et al. Reference Sheard, Thompson and Hourigan2003).
In terms of the VIV instability modes, mode A shows distinct differences from the 
$k'$ and 
$\dot {K}'$ distributions of the stationary cylinder modes. Under oscillation, there appears to be little generation associated with the vortex cores and more generation in the braids. In addition, further downstream, the instability mode appears supported by energy transfer in the vortex cores, while for the stationary cylinder mode A the generation remains small. The study by Leontini et al. (Reference Leontini, Thompson and Hourigan2007) indicates that the mode A seen for a stationary cylinder wake is continuously transformed into the mode A seen under high-amplitude oscillations as the forced oscillation amplitude is steadily increased. Thus, an interpretation of the perturbation kinetic energy analysis is that the stretching of vortex structures in the near-wake caused by the body oscillation suppresses the presence of near-wake elliptic instability, which causes the mode A instability to occur at considerably higher Reynolds numbers (
$Re_c\simeq 305$ at 
$U_r=5$) than for the stationary case (
$Re_c \simeq 190$).
For mode B similar spatial distributions of 
$k'$ and 
$\dot {K}'$ are seen for both the stationary and oscillating cylinder modes. Given this, it is not surprising that the instability mode onset is similar in both cases: 
$Re_c = 259$ for the stationary cylinder and 
$Re_c=277$ at 
$U_r=5$ for the oscillating cylinder. Thus, mode B is relatively minimally affected by large amplitude oscillations and the non-trivial change to the wake structure, noting a significantly increased wake width.
As shown later in § 3.6, a QP mode is the first occurring VIV instability mode at 
$U_r=7$. The comparison of the stationary and VIV QP modes in figure 11 indicates that the VIV QP mode more closely resembles mode B than the stationary cylinder QP mode. Hence, it is unsurprising that the critical Reynolds number for these modes is so different (
$Re_c\simeq 377$ for the stationary cylinder mode and 
$Re_c\simeq 220$ for the VIV mode).
Further discussion on how these modes change with reduced velocity is given in § 3.6.
3.4. Response curves for a circular cylinder under VIV at ${\textit {Re}} = 150$, 200 and 250
The VIV response of a circular cylinder for an undamped system is shown in figure 12 at 
$m^* = 2.546$ for 
${\textit {Re}} = 150$, 200 and 250. The following variables, the maximum oscillation amplitude (
$A^*$), peak lift coefficient (
$C_{L,max}$), normalised frequency ( 
$f^*$), total phase difference (
$\phi _{tot}$) and vortex phase difference (
$\phi _{vor}$), are plotted as a function of reduced velocity (
$U_r)$ in figure 12(a–d), respectively. Normalised frequency is defined as 
$f^* = f_y/f_n$, where 
$f_y$ is the frequency corresponding to the highest power in the power spectral density of displacement, and 
$f_n$ is the natural frequency of the structural system. The total phase 
$(\phi _{tot})$ and the vortex phase (
$\phi _{vor}$) are the phase differences between fluid force and displacement signals, obtained using a Hilbert transform (Khalak & Williamson Reference Khalak and Williamson1999). In particular, 
$\phi _{tot}$ and 
$\phi _{vor}$ are the mean phase difference between lift force and cylinder displacement, and the vortex force (total 
$-$ added mass force) and cylinder displacement, respectively.
VIV responses of the cylinder at 
$m^* = 2.546$ for 
${\textit {Re}} = 150$, 200 and 250: (a) maximum oscillation amplitude, 
$A^*$; (b) peak lift coefficient, 
$C_{L,max}$; (c) normalised frequency, 
$f^*\ (= f/f_n)$ (green dotted line represents the vortex shedding frequency for a stationary cylinder); (d) mean phase difference between lift force and cylinder displacement, 
$\phi _{tot}$; (e) mean phase difference between vortex force and cylinder displacement, 
$\phi _{vor}$.
We observe a typical three-branch response, the initial (red colour), upper (blue colour) and lower (black colour) branches, followed by a desynchronisation (magenta colour) regime. These branches are not necessarily one-to-one related to those observed for high-Reynolds number VIV at a low mass-damping ratio. Based on the observation by Govardhan & Williamson (Reference Govardhan and Williamson2000), Leontini, Thompson & Hourigan (Reference Leontini, Thompson and Hourigan2006b), Soti et al. (Reference Soti, Zhao, Thompson, Sheridan and Bhardwaj2018), Mishra et al. (Reference Mishra, Soti, Bhardwaj, Kulkarni and Thompson2020) and Sharma, Garg & Bhardwaj (Reference Sharma, Garg and Bhardwaj2022), the criteria to classify different VIV branches are as follows. The initial branch is demarcated with monotonically increasing 
$A^*$ and 
$C_{L,max}$, multiple frequencies in 
$f_y$, the vibration frequency overlaps the vortex shedding frequency of a stationary cylinder ( 
$f^* \approx f_v/f_s$), 
$\phi _{tot} \approx 0$ and 
$\phi _{vor} \approx 0$. The range of the upper branch is identified with large 
$A^*$, constant 
$f^*$ slightly lower than 1 ( 
$f^*<1$), large 
$C_{L, max}$, 
$\phi _{tot} \approx 0$ and transition of 
$\phi _{vor}$ from 0 to 
${\rm \pi}$. Beyond this, conditions for the existence of the lower branch is large 
$A^*$, 
$f^* =1$, monotonic decrease in 
$C_{L, max}$, transition of 
$\phi _{vor}$ from 0 to 
${\rm \pi}$ and 
$\phi _{vor} \approx {\rm \pi}$. Finally, desynchronisation is characterised by negligible 
$A^*$, 
$f^*$ parallel to 
$f_v/f_s$, 
$\phi _{tot} \approx {\rm \pi}$ and 
$\phi _{vor} \approx {\rm \pi}$.
As evident from figure 12(a) 
$A^*$ is a continuous function of 
$U_r$, except there are sudden jumps at 
$U_r = 3.3$, 3.1 and 3.0 for 
${\textit {Re}} = 150$, 200 and 250, respectively, thus allowing the identification of the transition between the initial and the upper branches. Although displacement amplitude is continuous, figure 12(b) shows a sudden jump in the lift-force amplitude, thus demarcating the upper and lower branches at 
$U_r = 4.0$ (
${\textit {Re}} = 150$), 3.9 (
${\textit {Re}} = 200$) and 3.8 (
${\textit {Re}} = 250$). Figure 12(c–e) shows the branch-to-branch jumps in the oscillation frequency and phase difference. We further quantify the branch ranges for 
${\textit {Re}} = 150$, 200 and 250 in table 1, together with the criteria used to distinguish the branches.
Range of 
$U_r$ and criteria distinguishing the different VIV branches for 
${\textit {Re}} = 150$, 200 and 250.
Figure 12(a) shows that the peak oscillation amplitude at all three Reynolds numbers is 
${\simeq }0.58D$, and the overall response curves are similar for 
${\textit {Re}} = 150$, 200 and 250. These findings agree with Govardhan & Williamson (Reference Govardhan and Williamson2006), who showed that peak oscillation amplitude is almost constant for low-Reynolds-number flows. However, Williamson & Govardhan (Reference Williamson and Govardhan2004) has noted that such low-
${\textit {Re}}$ simulations, in general, predict peak oscillation amplitudes around 
$0.6D$, whereas experiments for higher 
${\textit {Re}}$ produce amplitudes up to 1
$D$. Figure 13 shows wake vorticity contours at a reduced velocity, 
$U_r = 4.0$, corresponding to the maximum oscillation amplitude. At 
${\textit {Re}} = 150$, the wake vortices assemble into a double-row, 
$2S$ (two single vortices per shedding period), configuration consistent with the high amplitude oscillation (figure 13a). At 
${\textit {Re}} = 200$ and 250 (figure 13b,c), the newly shed vortices show signs of a 
$2P$ (two vortex pairs per shedding period) structure close to the rear of the cylinder, although the second vortex of each pair is very weak. However, we certainly do not claim that the wakes have a 
$2P$ structure for any of these Reynolds numbers, since the experimentally observed 
$2P$ structures show approximately similar strength opposite-signed vortices in each pair.
Vorticity contours at 
$U_r = 4.0$ for Reynolds number (a) 
${\textit {Re}} = 150$, (b) 
${\textit {Re}} = 200$ and (c) 
${\textit {Re}} = 250$. The red (blue) colours represents clockwise (anticlockwise) spanwise vorticity.
On this point, Blackburn, Govardhan & Williamson (Reference Blackburn, Govardhan and Williamson2001) compared low Reynolds number experimental results (
$Re \sim O(1000)$) with 2-D and 3-D numerical predictions. The 3-D simulations matched the experiments well. Furthermore, demarcation of the upper and the lower branches based on amplitude was not observed for the 2-D simulations, noting that they were clearly evident by the gaps between branches in the 3-D simulation. Therefore, and as has been recognised for a long time, 2-D simulations are inadequate for accurately modelling higher-
${\textit {Re}}$ VIV, when the flow is 3-D and turbulent.
3.5. Effect of reduced velocity on different modes
In § 3.2, stability analysis for VIV is presented for an oscillation frequency of 
$f_n^* = 0.2$ at 
${\textit {Re}} = 300$, corresponding to the strongly resonant case. In the current section, results for fixed 
${\textit {Re}} = 300$ for a standard spring–damper VIV system are presented as a function of reduced velocity. Figure 14 shows the maximum amplitude of vibration vs reduced velocity at 
${\textit {Re}} = 150$ and for 
$m^*=2.546$ depicting the periodic and QP nature of the oscillations. The situation is similar at other Reynolds numbers. Since FSA is only applicable to a periodic base flow, only the reduced velocity range 
$U_r \in [4.0, 8.0]$ is considered, noting periodic flow together with a strong resonant response.
Variation of amplitude (
$A^*$) vs reduced velocity (
$U_r$) for undamped VIV in the range 
$U_r \in [2.0, 8.0 ]$, depicting the periodic and QP range.
Figure 15 presents 
$|\mu |$ as a function of spanwise wavelength for 
$U_r \in [4.0- 8.0]$. The 
$U_r = 4.0$ case has a complex Floquet multiplier corresponding to a stable 
$2S$ wake mode, as shown in figure 16. At 
$U_r = 4.0$ the oscillation amplitude 
$A^* = 0.56$, corresponds to the highest amplitude during VIV. Interestingly, Leontini et al. (Reference Leontini, Thompson and Hourigan2007) observed a 
$P+S$ (base-flow) wake for forced oscillations at 
$A^* = 0.55$ and 
$f_s^* = 0.2$. That case corresponds to a forcing frequency of 
$St = 0.2$ or 
$U_r \simeq 5.0$. They observed the transition from a 
$2S$ to 
$P+S$ wake configuration at 
${\textit {Re}} = 280$. Thus, the system frequency, which controls the shedding frequency over the resonant range, does have a non-negligible effect on the base-flow structure, with increasing the forcing/shedding frequency leading to a higher-Reynolds-number transition from the 
$2S$ to the 
$P+S$ state.
Floquet multiplier vs wavelength, at 
${\textit {Re}} = 300, m^*=2.546$, 
$\zeta = 0$ for: (a) 
$U_r = 4.0$, (b) 
$U_r = 4.5$, (c) 
$U_r = 5.0$, (d) 
$U_r = 5.5$, (e) 
$U_r = 6.0$, ( f) 
$U_r = 6.5$, (g) 
$U_r = 7.0$, (h) 
$U_r = 7.5$ and (i) 
$U_r = 8.0$.
Base-flow vorticity contours at 
${\textit {Re}} = 300$ for reduced velocity 
$U_r = 4.0$. The red (blue) represents clockwise (anticlockwise) spanwise vorticity.
For 
$U_r \geq 4.5$, modes A, B, QP, QPL (a long-wavelength QP mode) and SH are observed. Of interest, 
$U_r = 4.5$ shows that a subharmonic mode (SH) is the first to become unstable. This occurs for a wavelength range typically dominated by Mode A. Modes B, QP and QPL are not unstable, whereas mode SH is unstable and grows fastest at 
$\lambda \approx 1.4D$. Therefore, at 
$U_r = 4.5$, mode SH is responsible for 3-D transition. It is important to note that the Floquet multiplier for the SH mode is always real and negative consistent with its subharmonic nature. The spanwise perturbation vorticity field is shown in figure 17. The subharmonic mode nature is highlighted by comparing images at each base-flow time period (
$T$), verifying that the perturbation field is repeated over two periods (
$2T$) and the structure of the mode is the same at the end of a single period except for a sign change (Sheard et al. Reference Sheard, Thompson and Hourigan2005).
Spanwise perturbation vorticity indicating mode SH for time (a) 
$T$, (b) 
$2T$, (c) 
$3T$ and (d) 
$4T$ at 
$\lambda = 1.4 D$, 
$U_r = 4.5$, 
$Re = 300$ and 
$m^* = 2.546$. Here, 
$T$ is the base-flow field period.
For the 
$\lambda < 5 D$, mode QPL is only observed at reduced velocity 
$U_r = 4.5$ and 5.0, whereas mode QPL is absent for other reduced velocities 
$U_r \in [5.5, 8.0]$. It is also observed that as 
$U_r$ increases, the stability of the modes are affected. For 
$U_r = 4.0$ the wake is stable, whereas for 
$U_r = 4.5$ mode SH is just critical and modes B, QP and QPL can be observed. For 
$U_r=5.0$, mode A is observed for the first time but remains stable, however mode B is unstable. Furthermore, for 
$U_r \geq 5.0$, mode B is unstable (
$|\mu | > 1$) and 
$|\mu |$ for this mode increases up to 
$U_r = 6.5$. The maximum 
$|\mu | \approx 1.7$ occurs at 
$\lambda \approx 0.9D$ for 
$U_r = 5.0$, increasing to a maximum 
$|\mu | \approx 4.0$ at 
$\lambda \approx 1.0D$ for 
$U_r = 6.5$. The growth rate for mode B further decreases at 
$U_r= 7.0$, with a maximum of 
$|\mu | \approx 3.5$ at 
$\lambda \approx 1.0D$, which remains constant on further increasing 
$U_r$. The maximum 
$|\mu | \approx 3.5$ occurs at 
$\lambda \approx 1.0D$ for 
$U_r = 7.0$, 7.5 and 
$8.0$.
Similar behaviour can be observed for mode A. At 
$U_r = 5.5$, the maximum 
$|\mu | \approx 1.03$ occurs at 
$\lambda \approx 2.7D$, indicating that the mode is just supercritical. The growth rate increases on further increasing 
$U_r$ up to 
$U_r = 7.5$ with a maximum 
$|\mu | \approx 1.8$ at 
$\lambda \approx 3.6D$. Beyond 
$U_r \geq 7.5$, the growth rate remains constant. Mode QP is seen to be stable for 
$U_r \in [4.0, 6.0]$ and unstable for 
$U_r \geq 6.0$. The maximum 
$|\mu |$ increases up to 
$U_r = 7.0$, with a highest growth rate of 
$|\mu | \approx 1.9$ at 
$\lambda \approx 1.9 D$. On increasing 
$U_r$ beyond 7.5, the maximum 
$|\mu |$ slightly decreases with a maximum of 
$|\mu | \approx 1.5$ and 
$|\mu | \approx 1.3$ at 
$U_r = 7.5$ and 8.0, respectively, centred around 
$\lambda = 1.9 D$. The flow is observed to be more stable at 
$U_r = 4.0$, 4.5 and 5.0, the highly resonant region. The synchronisation between cylinder vibration and vortex shedding stabilises the flow and hence it is more stable over this range.
3.6. Critical Reynolds number and wavelength variation with reduced velocity
This section presents the maximum 
${\textit {Re}}$ for which the wake remains 2-D for the undamped VIV (standard spring–damper) system as a function of reduced velocity, 
$U_r$. We have focused our analysis on the lower branch of VIV, where the base flow exhibits periodic wake oscillations. In the desynchronisation branch, the amplitude is relatively small and periodic with the oscillation frequency close to but lower than the natural shedding frequency. On the other hand, for the initial and (pseudo-)upper branch, the oscillation is not periodic, so FSA, which requires a purely periodic signal, cannot be applied. Thus, we have limited our consideration to the resonant reduced-velocity range, 
$U_r \in [4.5, 8]$, where the base-flow amplitude variation is periodic, noting that at the upper end of this range corresponds to the desynchonisation branch.
Figure 18 shows the critical Reynolds number, 
${\textit {Re}}_{cr}$, and corresponding critical wavelength, 
$\lambda _{cr}$, as a function of reduced velocity, 
$U_r$, for an undamped VIV system. Here, 
${\textit {Re}}_{cr}$ is the value of the Reynolds number at which a given mode first becomes unstable, i.e. 
$\mu$ first exceeds 
$|\mu | = 1$. In addition, 
$\lambda _{cr}$ is the wavelength at which 
${\textit {Re}}_{cr}$ occurs. To find an accurate 
${\textit {Re}}_{cr}$ at each 
$U_r$, simulations were run for each 
${\textit {Re}}$ for a series of 
$\lambda$ in steps of 0.1. For each 
${\textit {Re}}$, 
$\mu$ locally varies quadratically about each local maximum (as shown in previous 
$\mu$ vs 
$\lambda$ plots). We determine the coordinate (
$|\mu |, \lambda$) for the local maximum. The value of 
${\textit {Re}}$ is then incremented in steps of five, noting down the maximum 
$|\mu |$ and the corresponding value of 
$\lambda$ for each 
${\textit {Re}}$. This process is repeated until resolving values that give a maximum 
$|\mu | = 1$. The values of 
${\textit {Re}}_{cr}$ are then found by interpolating 
${\textit {Re}}$ and corresponding maximum 
$|\mu |$ values using quadratic interpolation. The same process is repeated for each reduced velocity.
(a) Critical Reynolds number, 
${\textit {Re}}_{cr}$, and (b) critical wavelength, 
$\lambda _{cr}$, against reduced velocity for the undamped system at 
$U_r \in [4.5, 8.0]$, and 
$m^* = 2.546$.
Figure 18 shows the transition boundary for modes A, B, SH and QP, corresponding to reduced velocity 
$[4.5, 8.0]$. The critical Reynolds number, 
${\textit {Re}}_{cr}$, decreases with reduced velocity up to 
$U_r=7.25$. At 
$U_r = 4.5$, mode SH is the only mode that becomes critical, which occurs at 
${\textit {Re}} \approx 299$. For 
$U_r \in [5.0, 6.0]$, modes A and B become critical, with mode B being the leading 3-D mode. At 
$U_r = 6.25$, mode QP is critical first, with the order of inception of modes being mode QP, mode B and mode A. It is observed that for 
$U_r = 7.0$ and 7.5, there are two distinct QP modes (referred to as QP and QP
$_0$) together with mode A that become critical, with mode QP
$_0$ becoming critical first. Interestingly, mode QP
$_0$ does not become critical at 
$U_r=6.25$. The perturbation kinetic energy analysis (see figure 11 and the discussion in § 3.3) suggests that the physical triggering mechanism for this mode has more in common with mode B rather than mode QP. It also occurs for the same wavelength band as for mode B, however, the transition Reynolds number is much lower than observed for mode B. Also for this 
$U_r$ range, at 
$U_r=7.25$ and 
$Re>270$, the wake develops asymmetry about the centreline and the wake period doubles. This part of parameter space is where the transition to mode B would be expected to occur. Transition to the desynchronisation branch occurs between 
$U_r = 7.5$ and 8.0, resulting in the relatively sudden shifting of shedding and cylinder oscillation to the be close to the natural shedding frequency, and a significant decrease in oscillation amplitude. At 
$U_r = 8.0$, mode A precedes the occurrence of modes B and QP, with modes B and QP transitioning at approximately the same 
${\textit {Re}}$. Interestingly, the transition Reynolds numbers and wavelengths of modes A and QP at 
$U_r=8.0$ are consistent with the trends at lower 
$U_r$ values, despite the significant change in oscillation frequency and amplitude. Mode B, on the other hand, shows a significantly lower transition Reynolds number for this case.
The transition boundary, corresponding wavelength and amplitude for modes A, B, QP and SH, and for 
$U_r \in [4.5, 8.0]$, are summarised in table 2. These findings are in accordance with Leontini et al. (Reference Leontini, Thompson and Hourigan2007) for the forced oscillation case: for 
$A^* < 0.3$ mode inception order (with increasing 
${\textit {Re}}$) is A, B and QP; whereas for 
$0.3< A^*<0.55$, the order of mode inception is B then A; and for 
$A^*>0.55$ mode SH is the only unstable mode. Table 2 indicates that for 
$U_r \in [4.5, 7.5]$, 
$A^*$ decreases approximately from 0.52 to 0.42 with 
$U_r$, whereas for 
$U_r = 8.0$, 
$A^* \approx 0.11$. Mode SH occurs at 
$A^* \approx 0.52$ for 
$U_r =4.5$, slightly varying from the prediction of Leontini et al. (Reference Leontini, Thompson and Hourigan2007), who observed mode SH at 
$A^* >0.55$. For 
$U_r \in [5.0, 6.0]$ with 
$0.42 \leq A^* \leq 0.50$ results are in agreement with Leontini et al. (Reference Leontini, Thompson and Hourigan2007). It is also observed that for 
$U_r = 6.5$ and 7.5, 
$0.40 \leq A^* \leq 0.44$, mode QP first becomes critical, which is different from Leontini et al. (Reference Leontini, Thompson and Hourigan2007); however, of course, the system frequency is then significantly different from the forcing frequency of 
$St = 0.2$ used in that study. On the other hand, at 
$U_r = 8.0$ with 
$A^* \simeq 0.11$ and the shedding frequency is close to 0.2, modes A and B undergo transition at close to the values seen for the forced oscillation case: 
$Re_{cr}=243$ at 
$St=0.18$ (VIV) compared with 
$Re_{cr}=255$ at 
$St=0.20$ (forced) for mode B; and 
$Re_{cr}=212$ at 
$St=0.18$ (VIV) compared with 
$Re_{cr}=230$ at 
$St=0.20$ (forced) for mode A. Overall, there appear to be non-negligible effects of the oscillation frequency and amplitude on transition.
Transition boundary characteristics: corresponding critical Reynolds numbers, wavelengths and amplitudes for modes A, B, SH and QP. The post-simulation oscillation frequency is shown in the last column corresponding to the Reynolds number of the first transition.
The envelope of the transition lines in figure 18(a) shows the variation of critical Reynolds number with reduced velocity for which the flows remains 2-D. The critical Reynolds number decreases with reduced velocity over the entire range of 
$U_r$ considered. Specifically, 
${\textit {Re}}_{cr} \approx 299$ at 
$U_r = 4.5$, decreasing monotonically to 
${\textit {Re}}_{cr} \approx 208$ at 
$U_r = 7.5$.
Figure 18(b) shows the variation of 
$\lambda _{cr}$ with reduced velocity. In terms of the first modes to become unstable, the figure shows that 
$\lambda _{cr}\approx 1.4 D$ for mode SH at 
$U_r = 4.5$. For 
$U_r = 5.0$ to 6.5, 
$\lambda _{cr}\approx 0.9D$ for mode B, however, for mode A there is an increase in 
$\lambda _{cr}$ from 
$2.4D$ to 
$4.3D$. On further increasing 
$U_r$, 
$\lambda _{cr}$ for mode A decreases towards 
$3.9D$ at 
$U_r = 8.0$. Mode QP, which appears for 
$U_r \geq 6.5$, decreases from 
$\lambda _{cr}\approx 2.0D$ to 
$\lambda _{cr}\approx 1.0D$ over the interval 
$6.5 \leq U_r \leq 7.5$. It then increases to 
$\lambda _{cr}\approx 2.0D$ at 
$U_r = 8.0$. The critical values for 
$U_r = 8.0$, where the oscillation amplitude is only 0.11, are close to those for the non-vibrating case (Barkley & Henderson Reference Barkley and Henderson1996) as should be expected.
Overall, the synchronisation between cylinder vibration and vortex shedding has a significant role in maintaining the two-dimensionality of the flow. Clearly, high-amplitude oscillations delay onset of 3-D transition but there is also a distinct dependence on reduced velocity. In summary, for the standard elastically mounted VIV system, the maximum 
${\textit {Re}}$ for which the wake remains 2-D over the resonant and periodic reduced-velocity range 
$U_r \in [4, 8]$ is 208 and occurs at 
$U_r = 7.5$, corresponding to a critical wavelength of 
$\lambda _{cr} \approx 1.2D$. The maximum delay in the transition from 2-D to 3-D flow occurs for 
$U_r = 4.5$, with the wake remaining 2-D up to 
${\textit {Re}}_{cr} \approx 299$.
3.7. Influence of mass ratio on the amplitude response and transition
Figure 19 shows the influence of mass ratio on the amplitude response under VIV as a function of Reynolds number. The subfigures in the first column correspond to (a) 
$Re=150$ (where the wake is certainly 2-D), (b) 
$Re=200$ (with the wake on the verge of transition) and (c) 
$Re=300$ (where the wake is mostly 3-D). The response curves have a similar shape but with the onset reduced velocity of resonance shifted to the right as the mass ratio is increased. In addition, desynchronisation occurs at smaller 
$U_r$ values for higher mass ratios. The right-hand column provides the equivalent amplitude curves but plotting against the ratio of actual system oscillation frequency (and vortex shedding frequency under VIV) to vortex shedding frequency for a stationary cylinder. The displacement response is only plotted for the resonant 
$U_r$ range where the oscillation is periodic. Clearly, for each Reynolds number, there is an excellent collapse of the amplitude response for the wide range of mass ratios considered. As already commented, higher mass ratios result in the resonance dropping out at higher oscillation amplitudes as 
$U_r$ is increased; however, the shape of the curves remains the same. The only slight variation where this does not occur is in the neighbourhood of the maximum amplitude response for 
$Re=300$, where there are differences in the maximum response amplitude for the different mass ratios. The upshot of this collapse suggests that the transition scenario discussed in the previous section should hold for other mass ratios by using the post-simulation measured oscillation frequency as the scaling parameter rather than 
$U_r$. To this end, table 2 also provides the oscillation frequency for the transitions for 
$m^*=2.546$. Thus, to find the equivalent transition for, say, 
$m^*=1.0$, that occurs at 
$U_r=6$ for 
$m^{*}=2.546$ (initially to mode B at 
$Re=261.3$), simulations can be undertaken at 
$m^*=1.0$ for 
$Re=261.3$ varying 
$U_r$ to find when 
$f(D/U)= 0.1729$. For 
$m^*=1$, this occurs for 
$U_r=6.37$, and indeed determining the Floquet multiplier for this 
$Re$–
$U_r$ combination for 
$\lambda /D = 1$ gives 
$\mu = 1.0$ in agreement with the 
$m^*=2.546$ result.
(a–c) Amplitude response for different Reynolds numbers and mass ratios. (d–f) Corresponding amplitude response against 
$f/f_v$ for the same cases in the left column.
3.8. Differences in 2-D and 3-D amplitude responses
In order to show the differences between 3-D and 2-D VIV simulations, even at relatively low Reynolds numbers, direct 3-D numerical simulations were undertaken. Figure 20 illustrates the comparison between 2-D and 3-D simulations at 
${\textit {Re}} = 300$ for 
$m^* = 2.546$. This particular Reynolds number choice is interesting because the previous transition analysis shows that the wake should be 3-D except for the highest-amplitude cases. Indeed, the predicted maximum amplitudes for the 3-D simulations centred in the neighbourhood of 
$U_r=4$ agree for predictions from 2-D simulations. A check on these cases does indicate that the wake remains 2-D there. However, there are significant differences between the amplitude response at smaller and larger reduced velocities. In particular, the predicted 3-D amplitudes between 
$4.5 \lesssim U_r \lesssim 8$ are significantly larger than their 2-D counterparts and the response branch is beginning to show the typical response observed in high-Reynolds-number experiments (where the lower branch maintains an almost constant amplitude until finally dropping at higher reduced velocity after desynchronisation (e.g. Khalak & Williamson Reference Khalak and Williamson1997). Although less noticeable, the initial branch amplitude is also over-predicted in the 2-D simulation. Both the 2-D and 3-D simulations fail to predict the upper branch, where the response amplitude is largest.
Amplitude response for 2-D and 3-D simulations with respect to 
$U_r$ at 
${\textit {Re}} = 300$ for 
$m^* = 2.546$.
3.9. Amplitude response as the wake becomes increasingly turbulent
We also used 3-D simulations to examine the response curves as the flow becomes increasingly turbulent, especially to document the appearance of the upper branch and to further determine differences from 2-D simulation predictions. The response curves are shown in figure 21 for 
$Re=300$, 500 and 1000. These predictions are compared with the experimental response curve for 
$m^*=2.4$ of Khalak & Williamson (Reference Khalak and Williamson1997) for 
$2000 < Re < 12\,000$, and numerical simulation results from Zhao et al. (Reference Zhao, Cheng, An and Lu2014) at a slightly smaller mass ratio of 
$2$. There is excellent agreement between the current predictions and those of Zhao et al. (Reference Zhao, Cheng, An and Lu2014), noting the effect of the smaller mass ratio is to extend the reduced velocity range of the lower branch slightly. Both sets of simulations predict a noticeable upper branch. The amplitude in the upper branch is below that seen in experiments, which is consistent with the previous observed dependence of the upper branch amplitude with Reynolds number for VIV (Govardhan & Williamson Reference Govardhan and Williamson2006; Rajamuni et al. Reference Rajamuni, Thompson and Hourigan2018). On reducing the Reynolds number to 
$Re=500$, the upper branch disappears. This was also observed by Blackburn et al. (Reference Blackburn, Govardhan and Williamson2001) for simulations for 
$Re = O(500)$. That study focused on demonstrating that shedding in the lower branch was of type 
$2P$, consistent with experiments unlike the typical 
$2S$ wake seen over most of the resonant range in lower-
$Re$ simulations.
Amplitude response for 2-D and 3-D simulations with respect to 
$U_r$ at 
${\textit {Re}} = 300$ for 
$m^* = 2.546$.
In terms of the lower branch wake state observed in the current set of simulations, figure 22 provides a comparison of wake vortices at 
$Re=300$ for the 2-D (a) and 3-D (b) simulations for 
$U_r=7.5$. Thus, development of 3-D flow changes the wake from 
$2S$ (for 2-D) to 
$P+S$ (for 3-D), which can be associated with the higher lower-branch response amplitude in the latter case. On increasing the Reynolds number to 
$Re=500$, the lower-branch wake changes to a 
$2P$ mode, as shown in figure 23 and consistent with Blackburn et al. (Reference Blackburn, Govardhan and Williamson2001). The top two images display spanwise-averaged wake vorticity contours for 
$Re=300$ (left) and 
$Re=500$ (right). The bottom row of images shows a top down view of these wakes visualised using isosurfaces of constant 
$Q$-criterion (
$Q=0.1$), which highlights rotational vorticity rather than shear. Clearly, the 
$Re=300$ wake is much less complex, although there is still significant distortion of the spanwise rollers. At 
$Re=500$, the wake is effectively turbulent while still showing the remnants of the instability modes examined in this paper.
Effect of 3-D transition on the lower branch. (a) Spanwise vorticity contours at 
$Re=300$, 
$m^*=2.546$, 
$U_r = 7.5$ for 2-D flow. (b) Same for 3-D flow at a representative spanwise cross-section.
Mode selection for the lower branch as Reynolds number is increased. (a) Span-averaged spanwise vorticity contours at 
$m^*=2.546$, 
$U_r = 7.5$ and 
$Re=300$ for 3-D flow. (b) 
$P+S$ wake mode for this case depicted using isosurfaces of 
$Q$-criterion (
$Q = 0.1$). (c) Same as (a) for 
$Re=500$, showing the mode has switched to 
$2P$. (d) Corresponding 
$Q$-criterion contours for this case.
4. Conclusions
FSA of flow past a circular cylinder undergoing VIV subject to standard elastic support has been undertaken. It has been found that for 
${\textit {Re}}= 300$, all three previously identified stationary cylinder wake modes, A, B and QP, can be observed. In this case, 3-D transition is through mode B for a spanwise wavelength of 
$\lambda \approx 0.9D$, close to the value observed for a stationary cylinder, whereas modes A have Floquet multipliers less than unity and so remain stable at this transition 
${\textit {Re}}$. Over the resonant reduced-velocity range 
$U_r \in [4.5, 8]$ corresponding to purely periodic shedding, the 3-D onset mode changes from mode B to mode A via mode QP as 
$U_r$ is increased. Interestingly, for 
$U_r =4.5$, a subharmonic (SH) mode is observed and is found to be the critical mode. In all other cases, mode A is observed instead of mode SH. This may be due to the higher amplitude and lower frequency of oscillation as compared to the other reduced-velocity cases examined. The critical 
${\textit {Re}}$ for 3-D transition has been observed to decrease as 
$U_r$ increases from 4.5 to 8, with the critical mode changing from mode B to mode A.
Thus, each of the modes A, B, QP and SH have been shown to drive the wake transition from 2-D to 3-D flow depending on the amplitude of oscillation and reduced velocity considered. The maximum 
${\textit {Re}}$ for which the wake is guaranteed to remain 2-D for reduced velocities covering the resonant response range, 
$U_r \in [4.5, 8.0]$, is 202. The 3-D simulations at 
$Re=300$ show that the amplitude response curve is distinctly different from that obtained using 2-D simulations. This study improves our understanding of 3-D wake transition for a circular cylinder undergoing VIV under standard elastic mounting, documenting how response oscillation amplitude and reduced velocity effect transition.
Funding
This research was made possible by a computer time grant (NCMAS d71) through the National Computing Infrastructure (NCI) initiative. We would like to acknowledge these computer resources from NCI and the Pawsey Supercomputing Centre that have been generously supported by the Australian Government and the Government of Western Australia. We also acknowledge support from an Australian Research Council Discovery Grant DP190103388.
Declaration of interests
The authors report no conflict of interest.
Open access
