1. Introduction
Finite-length cylinders and their variants are ubiquitous in human lives and engineering, including the cylinder with one free end, e.g. chimneys, cylindrical tall buildings, etc. and the cylinder with two free ends, e.g. submarine-like shape (Tezuka & Suzuki Reference Tezuka and Suzuki2006; Sheard, Thompson & Hourigan Reference Sheard, Thompson and Hourigan2008), torpedoes-like shape (Schouveiler & Provansal Reference Schouveiler and Provansal2001), wheels (Zdravkovich et al. Reference Zdravkovich, Brand, Mathew and Weston1989) and the short cylindrical bluff bodies (Prosser & Smith Reference Prosser and Smith2016; Yang et al. Reference Yang, Guo, Liu, Chen, Xing and Zhao2021) etc. However, compared with the extensively studied flow past an infinite cylinder, there are relatively fewer works on the flow past a finite-length cylinder with two free ends and, thus, the current understanding of the flow dynamics past a finite-length cylinder is insufficient. The works on the flows past finite-length cylinders in the literature also appear scattered, mainly because of the many different configurations that are possible for the cylinder. For example, the cylinders can be classified as (1) having two flat free ends vs two non-flat free ends; (2) having two free ends vs one free end and one fixed end (connected to a ground plane); (3) being axially perpendicular to the incoming flow vs axially parallel to the incoming flow, etc. Plus, important parameters, such as Reynolds number (
$Re$), aspect ratio (
${\small \text{AR}}$), yaw angle, etc., are all different in each work, further enriching but at the same time complicating the current literature. In this work we will focus on the flow past a cylinder with two free ends and with an axis perpendicular to the incoming flow with 
$0.5\leq {\small \text{AR}} \leq 2$ and 
$Re\le 1000$. In the following, we will summarise the relevant works in the literature on the studies of the flows past a finite-length cylinder and on stability analyses based on a time-mean flow. At the end of this section, we will clarify the position of the current work by identifying the research gap we are going to fill.
1.1. Flows past a finite cylinder with two free ends
In an oil shadow experiment to visualise the flow past a finite-length cylinder, Zdravkovich et al. (Reference Zdravkovich, Brand, Mathew and Weston1989) found that periodical vortices can be observed at the aspect ratio of 
$2<{\small \text{AR}}<8$ and 
$6\times 10^4< Re<2.6\times 10^5$ (where 
${\small \text{AR}}=L/D$ with 
$L$ being the length or height and 
$D$ the diameter of the cylinder), and when 
${\small \text{AR}}$ is reduced to 3, the eye-like low pressure area near the free ends gradually disappears. Moreover, when 
${\small \text{AR}}<3$, the pressure distribution on the curved surface is no longer symmetric about the centre plane, and this asymmetrical flow produces yaw and roll moments. In general, the drag coefficient decreases as the aspect ratio decreases. However, Zdravkovich et al. (Reference Zdravkovich, Brand, Mathew and Weston1989) later conducted experiments on the flow around a coin-like cylinder with 
$0< {\small \text{AR}}\leq 1$ and showed that in the range of 
$2\times 10^5< Re<6\times 10^5$, the drag coefficient and the aspect ratio satisfy the relationship 
$\bar {C}_d=0.024/{\small \text{AR}}+0.633, (2\times 10^5< Re<6\times 10^5)$, that is, when 
${\small \text{AR}}$ decreases, 
$\bar {C}_d$ increases. This trend contradicted their previous measurements (in Zdravkovich et al. Reference Zdravkovich, Brand, Mathew and Weston1989), which the authors attributed to the inappropriate reference area used in this range of 
${\small \text{AR}}$ (once the side area 
$({D^2}/{4}){\rm \pi}$ was used as the reference area, instead of the projected area 
$LD$, the previous trend that 
$\bar {C}_d$ decreases as 
${\small \text{AR}}$ decreases held valid for small values of 
${\small \text{AR}}$). Furthermore, the topology of the flow field has also been depicted: there are two horseshoe-shaped vortices on the two free ends, which are separated at an angular position of approximately 
$90^{\circ }$, and the free vortex forms two counter-rotating vortex pairs in the leeward zone. Based on this observation, Zdravkovich analysed and proposed a drag-reducing strategy of rounded sharp edges with a drag reduction effect up to 33.6 %.
More recently, direct numerical simulations (DNS) have been adopted in this field of research; for example, Inoue & Sakuragi (Reference Inoue and Sakuragi2008) studied the flow around a finite cylinder with two free ends with 
$0.5\leq {\small \text{AR}}\leq 100$ and 
$40\leq Re\leq 300$. Their results showed that the 
${\small \text{AR}}$ and 
$Re$ have a great influence on the vortex shedding pattern, and five types of vortex shedding patterns were identified as: (I) periodic oblique vortex shedding (large 
${\small \text{AR}}$, 
$Re>Re_c$); (II) quasi-periodic oblique vortex shedding (large 
${\small \text{AR}}$, 
$Re< Re_c$); (III) hairpin vortex that periodically falls off (moderate 
${\small \text{AR}}$); (IV) two stable counter-rotating vortex pairs (small 
${\small \text{AR}}$ and 
$Re$); (V) the counter-rotating vortex pairs alternately shed from the free ends (small 
${\small \text{AR}}$, high 
$Re$). Prosser & Smith (Reference Prosser and Smith2016) conducted large-eddy simulations coupled with the unsteady Reynolds-averaged Navier–Stokes (NS) equations using a finite-volume method to study the flow around bluff bodies (finite cylinder or prism) with much higher 
$Re$ in the range of 
$[10^5, 10^6]$ and 
${\small \text{AR}}=1$ and 
$2$. The influences of the geometric features and attitude of the bluff body on the reattachment distance, pressure coefficient and stagnation point position were investigated and modelled empirically. For example, they found that through proper normalisation, on the flat ends of the bluff body, the pressure coefficient and the stagnation point position are single-valued functions of the incident angle; but on the curved face of the bluff body, these two variables are affected by the aspect ratio. Gao et al. (Reference Gao, Nelias, Liu and Lyu2018) numerically studied the flow pattern of a cylinder with two free ends. Unlike the flow field of a cylinder with one free end (where there is a large horseshoe vortex surrounding the fixed end of the cylinder), there is no horseshoe vortex that appears in the flow past the cylinder with two free ends, and a new relationship between 
$\bar {C}_d$ and 
$Re$ has been proposed within the range of 
$5\leq Re \leq 5\times 10^6$, 
$1\leq {\small \text{AR}} \leq 6$. In addition, Pierson et al. (Reference Pierson, Auguste, Hammouti and Wachs2019) used the finite-volume–fictitious-domain method to study the axial flow (i.e. the cylinder axis is parallel to the streamwise direction) around a finite-length three-dimensional (3-D) cylinder (
${\small \text{AR}}=1, 20< Re<460$) with two free ends. They identified two bifurcation points both without hysteresis, one is the regular bifurcation at 
$Re\approx 278$ and the other is the Hopf bifurcation at 
$Re\approx 355$. Unfortunately, they did not explain how the aspect ratio affects these two bifurcation points.
In addition, finite-length cylinders with non-flat ends have also been studied by some groups. These works are reviewed here because some of them are inspiring and comparisons have been made to these works in our result section. Zdravkovich et al. (Reference Zdravkovich, Brand, Mathew and Weston1989), Schouveiler & Provansal (Reference Schouveiler and Provansal2001) and Sheard, Thompson & Hourigan (Reference Sheard, Thompson and Hourigan2005) used wind tunnel experiment and numerical simulations, respectively to study flow around a finite cylinder with two hemispherical ends. The results show that the hemispherical ends can effectively reduce drag coefficient. Since the flow field past this type of cylinder has no fixed boundary layer (due to the rounded ends), there is no horseshoe vortex in the flow field either. When the two hemispherical ends connect each other, a sphere is formed. The wake pattern, transition and instability of the flow around a sphere have been studied for a long time, especially regarding the critical Reynolds number, beyond which the steady flow becomes unsteady. This critical Reynolds number has been reported as 
$Re_c=130$ in early experimental research (Taneda Reference Taneda1956), as 
$Re_c=175$ in linear stability analysis (LSA) (Kim & Pearlstein Reference Kim and Pearlstein1990), as 
$Re_c=270$ in experiments of Magarvey & MacLatchy (Reference Magarvey and MacLatchy1965), Wu & Faeth (Reference Wu and Faeth1993) and Johnson & Patel (Reference Johnson and Patel1999), as 
$Re_c=277.5$ in LSA of Natarajan & Acrivos (Reference Natarajan and Acrivos1993) and finally as 
$270< Re_c<285$ in the DNS using the spectral element method (SEM) by Tomboulides & Orszag (Reference Tomboulides and Orszag2000). In summary, the value of Hopf bifurcation 
$Re_c=277.5$ obtained by LSA has been generally accepted. This is the second transition in the sphere flow; and the first transition occurs around 
$Re=210$ (Natarajan & Acrivos Reference Natarajan and Acrivos1993) or 
$Re=212$ (Tomboulides & Orszag Reference Tomboulides and Orszag2000). In this work we will focus on the wake patterns and transition of the flow around a short finite cylinder (
${\small \text{AR}}\leq 2$). As we will later show, the flow past a short finite-length cylinder (
${\small \text{AR}} \leq 1.7$) shows a transition process similar to that of a sphere, but also exhibits some unique characteristics.
1.2. Stability analyses of a mean flow
As we will conduct global stability analysis of the flow around the finite cylinder, it is pertinent to review works on the application of this analysis to the infinite-length cylinder, pioneered by Jackson (Reference Jackson1987), Zebib (Reference Zebib1987) and Noack & Eckelmann (Reference Noack and Eckelmann1994). Especially, we will review the stability analyses applied to the two-dimensional (2-D) wake flows based on a time-mean flow.
For the wake flow behind an infinite-length (2-D) cylinder, Hammond & Redekopp (Reference Hammond and Redekopp1997) and Pier (Reference Pier2002) noted that applying the LSA to the time-averaged flow field (which will be called mean flow for simplicity in the following) instead of the steady base flow, that is an unstable solution to the NS equations, can better predict the shedding frequency of the unsteady flow. Barkley (Reference Barkley2006) conducted a global LSA of the 2-D cylinder flow based on its mean flow and observed that the eigenfrequency is consistent with the nonlinear vortex shedding frequency 
$St$ found by Williamson (Reference Williamson1988) in his famous experiments of the cylinder flow. In addition, Barkley (Reference Barkley2006) also showed that the mean flow is marginally stable. This striking result was explained in the asymptotic analysis by Sipp & Lebedev (Reference Sipp and Lebedev2007) who demonstrated that the 2-D cylinder flow can be well represented by the mean flow value and a single eigenmode, and there is almost no nonlinear interaction between them, validating the usage of the mean flow in the stability analysis. Therefore, the premise of Barkley (Reference Barkley2006) for the global stability analysis of using the time-averaged flow was verified, that is, the Reynolds stress caused by the pulsating wake is not disturbed at the linear order. Sipp & Lebedev (Reference Sipp and Lebedev2007) also provided a theoretical proof for the validity of the global weakly nonlinear analysis of the 2-D cylinder flow near the Hopf bifurcation. More specifically, the theoretical conditions were given on how to keep the mean flow linearly stable and when the eigenfrequency obtained using the mean flow matches the experimental frequency. This condition can be qualitatively described as a situation where the zeroth harmonic (mean flow) is much stronger than the second harmonic. Meanwhile, they also pointed out that the effective results obtained in the case of the 2-D cylinder based on the mean flow are by no means general (in the context of the finite-length cylinder whose axis is perpendicular to the flow direction, the stability analysis based on its mean flow has not been performed and will be conducted here). Besides, because of the weakly nonlinear expansion employed in the theoretical development, the analysis of Sipp & Lebedev (Reference Sipp and Lebedev2007) is strictly valid only very close to the bifurcation point, and cannot fully explain the success of Barkley (Reference Barkley2006)'s global LSA performed on the mean flow across a wide range of Reynolds numbers.
Later, Leontini, Thompson & Hourigan (Reference Leontini, Thompson and Hourigan2010) applied the saddle-point criterion (that the zero group velocity is found at the saddle point in the complex wavenumber plane or the cusp point in the complex frequency plane when local mean flows are analysed in the streamwise direction) to spanwise-averaged 3-D infinite cylinder flow. Their results showed that if the local curvature is not too large and the assumption of the flow changing slowly is valid, even if the flow is three dimensional, the saddle-point criterion can work well. Their global LSA showed that even for the 3-D flow, the spanwise-averaged mean flow remains marginally stable, which supports the following hypothesis: the dynamics of the cylinder wake is mainly dominated by the first linear eigenmode generated on the nonlinear modified mean flow. Turton, Tuckerman & Barkley (Reference Turton, Tuckerman and Barkley2015) proved in a more general manner that if the flow exhibits quasi-monochromatic oscillations, the eigenfrequency of the linearisation operator based on the temporally mean flow is indeed equal to the true nonlinear flow frequency. They further looked into the amplitudes of the waves in the spectral space of the signal and demonstrated that if the amplitude of the first harmonics is sufficiently larger than those of higher-order harmonics, the stability analysis based on the mean flow can represent fairly well the dynamics of the nonlinear flow.
For a weakly non-parallel and strongly convectively unstable flow, whose first singular value largely dominates the others, a strong link between the mean flow and the fully nonlinear dynamics of a turbulent flow has been explained and verified by Beneddine et al. (Reference Beneddine, Sipp, Arnault, Dandois and Lesshafft2016, Reference Beneddine, Yegavian, Sipp and Leclaire2017) from both theoretical and experimental perspectives. According to the RZIF property (real zero imaginary frequency named by Turton et al. Reference Turton, Tuckerman and Barkley2015), Bengana et al. (Reference Bengana, Loiseau, Robinet and Tuckerman2019) identified two successive Hopf bifurcations and classified three situations by using mean flow and base flow LSA in the shear-driven cavity flow and summarizing previous research. Another interesting finding is that although no mathematical proof is given, the cross-eigenvalues would then correspond to the relative stability of the first two eigenmodes, which can qualitatively explain well the hysteresis in the nonlinear flow system.
In summary, it can be concluded that the global LSA based on the mean flow can predict accurately the nonlinear frequency in a wider range of parameters for the quasi-monochromatic oscillation flows. It remains to be seen if this method is applicable to the flow to be studied here.
1.3. The current work
As we can see from the previous sections, there is no regular bifurcation (so-called pitchfork) being reported in the previous research on 3-D finite-length cylinder flows with 
${\small \text{AR}}>2$. It has not been studied systematically how this flow bifurcates when 
${\small \text{AR}}<2$. On the other hand, in the flow around other bluff bodies, many works have confirmed the existence of a regular bifurcation, such as a sphere (Johnson & Patel Reference Johnson and Patel1999; Tomboulides & Orszag Reference Tomboulides and Orszag2000; Thompson, Leweke & Provansal Reference Thompson, Leweke and Provansal2001; Sheard, Thompson & Hourigan Reference Sheard, Thompson and Hourigan2004), axial flow around short cylinders with 
${\small \text{AR}}=1$ (Pierson et al. Reference Pierson, Auguste, Hammouti and Wachs2019) and ellipsoids (Tezuka & Suzuki Reference Tezuka and Suzuki2006; Sheard et al. Reference Sheard, Thompson and Hourigan2008). Thus, one can infer that the aspect ratio may be the decisive factor for educing the regular bifurcation in the finite-length cylinder flows and the regime 
${\small \text{AR}}\le 2$ should be explored to see if this bifurcation exists therein. In this work we will study systematically the influence of 
$Re(\leq 1000$) and 
${\small \text{AR}}\in [0.5,2]$ on the flow stability/instability of the flow past a 3-D finite-length cylinder by nonlinear DNS and linear global LSA (mainly based on temporal mean flow). We will examine whether the LSA results obtained based on the mean flow can be compared with the nonlinear results in the DNS (e.g. in terms of the shedding frequency in the flow). Besides, in order to complete the discussions, we will also compare the LSA results based on the mean flow with the LSA results based on a base flow obtained using the selective frequency damping (SFD) method (Åkervik et al. Reference Åkervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006). The most important result will be summarised in figure 21 on the 
$Re_c$–
${\small \text{AR}}$ relation.
The paper is organised as follows. Section 2 introduces the configuration of 3-D finite cylinder flow, the boundary conditions, the governing equations, i.e. nonlinear NS equations and their corresponding linearised equations about the base state (mean flow or base flow) and the numerical method. Section 3 provides a detailed verification step of nonlinear and linear numerical methods. In § 4 we show the results and discuss the base states, cylinder wake patterns, critical 
$Re$, nonlinear DNS bifurcation scenario, global eigenmodes and effects of 
${\small \text{AR}}$ and 
$Re$ on this flow. Finally, the results are summarised in § 5 and conclusions are provided.
2. Problem formulation and numerical methods
2.1. Governing equations and geometry
We study the 3-D stability of the flow around a finite-length cylinder (of a length 
$L$ and a diameter 
$D$) subjected to a uniform incoming flow in a Cartesian coordinate system. The computational domain size, boundary conditions and geometry of the finite cylinder are shown in figure 1. The origin of the Cartesian coordinate system is located at the centre of the cylinder, the 
$x$-axis points in the flow direction, the 
$z$-axis extends along the centreline of the cylinder and the 
$y$-axis is determined by the right-hand rule. The axis of the cylinder is normal to the incoming flow. The NS equations for the unsteady incompressible flow read as
\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}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{U} = 0, \end{equation}
where 
$\boldsymbol {U}$ is the velocity vector, whose components corresponding to the three directions of 
$x$, 
$y$, 
$z$ are 
$\boldsymbol {U}=(U_x,U_y,U_z)$ and 
$P$ is the pressure. In (2.1) the cylinder diameter 
$D$ is used as the reference length, the velocity 
$U_\infty$ of the uniform incoming flow at infinity as the reference velocity and 
$\rho U^2_\infty$ as the reference pressure. Therefore, the Reynolds number, Strouhal number, drag coefficient 
$C_d$, lift coefficient 
$C_{l}$ and aspect ratio are defined respectively as
\begin{equation} Re = \frac{U_\infty D}{\nu}, \quad St = \frac{fD}{U_\infty}, \quad C_d = \frac{F_d}{\frac{1}{2} \rho U^2_\infty A},\quad C_{l} = \frac{F_l}{\frac{1}{2} \rho U^2_\infty A}, \quad {\small \text{AR}} = \frac{L}{D}, \end{equation}
where 
$f$ is the frequency of vortex shedding (thus, when there is no vortex shedding in the wake, 
$St=0$), 
$\nu$ is the kinematic viscosity of the fluid and 
$\rho$ is the density. Here 
$F_d$ is the drag force on the surface of the cylinder, whose direction is the same as the streamwise direction; 
$F_l$ is the cross-stream lift force acting in either the 
$y$ or 
$z$ direction (see 
$C_{ly}$ and 
$C_{lz}$ below); 
$A$ is the reference area, and for a finite cylinder, 
$A=LD$. For the lift coefficient 
$C_l$, there are two possibilities in the 
$y$ and 
$z$ directions. Among them, 
$C_{ly}$ (in the 
$y$ direction) is mainly discussed in the present work, because its value is directly related to the regular bifurcation; 
$C_{lz}$ (in the 
$z$ direction) or its mean is almost zero in most of the cases due to the symmetry of the flow. Besides, we will also use the letters 
$\bar {C}_d$ and 
$\bar {C}_{l}$ to denote the time-averaged values of 
$C_d$ and 
$C_{l}$.
Figure 1. (
$a$) The computational domain and boundary conditions (not to scale). (
$b$) The projection of the mesh topology on the 
$x$–
$y$ plane.
Next, we define the boundary conditions in the computational domain. As shown in figure 1(
$a$), 
$S_c$ represents the surfaces of the finite cylinder, on which no-slip boundary conditions are applied,
\begin{equation} \boldsymbol{U} =(U_x,U_y,U_z)=0 \quad \textrm{on} \ {S_c}. \end{equation}
Here 
$S_{in}$ and 
$S_{out}$ represent the inlet and outlet surfaces of the rectangular computation domain, whose normal is in the 
$x$ direction. The inlet boundary conditions for the velocity imposed on 
$S_{in}$ are
\begin{equation} U_x=1,\quad U_y=U_z=0 \quad \textrm{on} \ S_{in}. \end{equation}
The pressure and the streamwise derivative of the 
$\boldsymbol {U}$ components are set to zero on the outlet 
$S_{out}$,
\begin{equation} \frac{\partial \boldsymbol{U}}{\partial x} =0,\quad P=0 \quad \textrm{on} \ S_{out}. \end{equation}
Here 
$S_{xy,t}$, 
$S_{xy,b}$, 
$S_{xz,f}$ and 
$S_{xz,b}$ represent the surfaces of the cuboid on the top, bottom, front and back side walls, which are parallel to the 
$xy$, 
$xy$, 
$xz$ and 
$xz$ planes, respectively. Previous research indicates that the computation domain has some influence on the results of the LSA. For example, Giannetti & Luchini (Reference Giannetti and Luchini2007) conducted an eigenvalue sensitivity analysis on the size of the computational domain in the problem of the 2-D cylinder wake flow with symmetric boundaries in the cross-streamwise direction. Their results showed that when the size in the cross-streamwise direction is greater than 
$10D$, the drift of the eigenvalue and eigenmode is small. Following them, we will perform a similar analysis. The symmetry boundary condition is imposed on the surfaces 
$S_{xy,t}$, 
$S_{xy,b}$, 
$S_{xz,f}$ and 
$S_{xz,b}$ at position 
${\pm }10D$, that is, 
$L_y=L_z=10$,
\begin{equation} \left. \begin{array}{ll@{}} \dfrac{\partial U_x }{\partial y}=U_y=\dfrac{\partial U_z }{\partial y}=\dfrac{\partial P}{\partial y}=0 & \textrm{on} \ S_{xz,f},S_{xz,b}, \\ \dfrac{\partial U_x }{\partial z}=\dfrac{\partial U_y }{\partial z}=U_z=\dfrac{\partial P}{\partial z}=0 & \textrm{on} \ S_{xy,t},S_{xy,b}. \end{array} \right\} \end{equation}
The values of 
$L_{x1},L_{x2}$ define the position of the cylinder in the streamwise direction, which, together with the vertical length 
$L_{y}$ and span length 
$L_{z}$, will be discussed in a convergence study of the computational domain size on the linear dynamics.
2.2. Linearisation
In the present work we will study the global linear stability/instability of the flows around the finite-length cylinder. As we have reviewed in the introduction section, we will focus on studying two types of base states (both of which are denoted as 
$(\boldsymbol {U}_b,P_b)$ in the text to follow as long as there is no confusion) to analyse their global instability. The first base state is called a base flow, which is a steady solution to the NS equations (2.1) obtained by the SFD method proposed by Åkervik et al. (Reference Åkervik, Brandt, Henningson, Hœpffner, Marxen and Schlatter2006). This method damps the unsteady temporal oscillations by adding a dissipative relaxation term in the NS equations to let the low-frequency component of the velocity pass. The second base state is the mean flow, which is obtained by time averaging the periodic flow with vortex shedding. For the present global stability analysis, at least ten vortex shedding cycles will be used for the time-average procedure. The linearisation starts with Reynolds’ decomposition 
$\boldsymbol {U} = \boldsymbol {U}_b+\boldsymbol {u}, P=P_b+p$, which will be substituted into the nonlinear governing equations. Then, we discard the nonlinear terms and the terms satisfying the NS solution for the base states and retain terms of the order of the perturbation, yielding the linearised equations for the infinitesimal perturbations residing on these base states
\begin{equation} \frac{\partial \boldsymbol{u} }{\partial t}+(\boldsymbol{U}_b \boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u} + (\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{U}_b ={-}\boldsymbol{\nabla} p+ \frac{1}{Re} \nabla^2\boldsymbol{u}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0, \end{equation}
where 
$\boldsymbol {u}$ is the 3-D perturbation velocity vector 
$\boldsymbol {u}=(u_x,u_y,u_z)$ and 
$p$ is the perturbation pressure. Homogeneous boundary conditions are applied for the perturbed variables
\begin{gather} u_x=u_y=u_z=0 \quad \textrm{on} \ {S}_c, \end{gather}
\begin{gather}u_x=u_y=u_z=0 \quad \textrm{on} \ {S}_{in} , \end{gather}
\begin{gather}\frac{\partial \boldsymbol{u}}{\partial x} =0,\quad p=0, \quad \textrm{on} \ S_{out}, \end{gather}
\begin{gather}\frac{\partial u_x }{\partial y}=u_y=\frac{\partial u_z }{\partial y}=\frac{\partial p}{\partial y}=0 \quad \textrm{on} \ S_{xz,f},S_{xz,b}, \end{gather}
\begin{gather}\frac{\partial u_x }{\partial z}=\frac{\partial u_y }{\partial z}=u_z=\frac{\partial p}{\partial z}=0 \quad \textrm{on} \ S_{xy,t},S_{xy,b}. \end{gather}We found that the numerical results obtained by these boundary conditions in the eigensolver (to be discussed shortly) did not manifest the non-physical oscillations at the outlet boundary encountered by Giannetti & Luchini (Reference Giannetti and Luchini2007). So, we did not set the partial derivative of pressure in the streamwise to zero.
Equations (2.7) can be written in matrix form with 
$\boldsymbol q=(\boldsymbol {u}, p)$ as
\begin{equation} \boldsymbol{M} \frac{\partial \boldsymbol q }{\partial t} = \boldsymbol{A}\boldsymbol q, \end{equation}
where the mass matrix 
$\boldsymbol {M}$ and the Jacobian matrix 
$\boldsymbol {A}$ are given by
\begin{equation} \boldsymbol{M} =\left(\begin{array}{cc} \boldsymbol I & 0\\ 0 & 0 \end{array} \right), \quad \boldsymbol{A} = \left(\begin{array}{cc} -\boldsymbol{U}_b \boldsymbol{\cdot}\boldsymbol{\nabla} - \boldsymbol{\nabla}\boldsymbol{U}_b + Re^{{-}1}\nabla^2 & -\boldsymbol{\nabla}\\ \boldsymbol{\nabla}\boldsymbol{\cdot} & 0 \end{array} \right). \end{equation}
We expand the perturbation 
$\boldsymbol {q}$ in time 
$\boldsymbol {q}(x,y,z,t)=\hat {\boldsymbol {q}}(x,y,z)e^{\lambda t}$, where 
$\lambda =\sigma +\boldsymbol {i}2{\rm \pi} \omega$. Equation (2.9) can be transformed into the following generalized eigenvalue problem by substituting the solution ansatz of 
$\boldsymbol {q}$:
\begin{equation} \lambda \boldsymbol{M} \hat{\boldsymbol{q}} = \boldsymbol{A}\hat{\boldsymbol{q}}. \end{equation}
In this formulation, the stability of the base state 
$\boldsymbol {U}_b$ is dictated by the eigenvalues 
$\lambda$ in the linearised problem with 
$\hat {\boldsymbol {q}}$ being the eigenmode, 
$\sigma$ the temporal growth/decay rate of perturbations and 
$\omega$ the eigenfrequency. Because the flow problem is not spatially periodic or homogeneous in either 
$x,y,z$ directions, 
$\hat {\boldsymbol {q}}$ depends on all the three coordinates, giving rise to a global stability problem (Theofilis Reference Theofilis2011) to be solved by the Arnoldi method. In addition, the eigenfrequency 
$\omega$ of the first eigenvalue also determines whether the base state 
$\boldsymbol {U}_b$ experiences a regular bifurcation (
$\omega =0$) or a Hopf bifurcation (
$\omega >0$) (Bengana et al. Reference Bengana, Loiseau, Robinet and Tuckerman2019).
2.3. Numerical methods: DNS and implicitly restarted Arnoldi method
In order to obtain the accurate wake pattern and base state of the fully 3-D flow past a finite-length cylinder at medium and low Reynolds numbers, we conduct DNS of the flow past the cylinder. The high-order parallelised open-source code Nek5000 (Fischer, Kerkemeier & Peplinski Reference Fischer, Kerkemeier and Peplinski2020) (version 19.0) is used, which is based on the nodal SEM originally proposed by Patera (Reference Patera1984). Hexahedral elements for a polynomial order 
$N = 7$ (the optimal polynomial order suggested by Fischer et al. Reference Fischer, Kerkemeier and Peplinski2020) are used to get a better performance of the code. The time step 
$\Delta t$ is determined by the Courant–Friedrichs–Lewy condition with the target Courant number being in the range of 
$0.5$–
$1.0$. In the following numerical simulations, the boundary layer elements have been refined by the O-type. To achieve the requirement of adequate resolution near the surface of the cylinder, the value of the smallest boundary thickness for each Reynolds number will be discussed in § 3.1.
For a non-parallel 3-D flow past the finite cylinder, the numerical discretisation of linearised NS equation (2.7) will results in a large-scale Jacobian matrix 
$\boldsymbol {A}$ in the generalized eigenvalue problem (2.11). It is often impractical to solve a large-scale eigenvalue problem for its whole eigenspectrum in a 3-D flow. Based on Nek5000 solver and the ARPACK package (Lehoucq, Sorensen & Yang Reference Lehoucq, Sorensen and Yang1998), the matrix-free time-stepper method (Theofilis Reference Theofilis2011; Doedel & Tuckerman Reference Doedel and Tuckerman2012) is used in the present work. In this method one does not need to explicitly construct the matrices, but instead power iterate the linearised NS equations (2.7) in the temporal direction from an initial condition 
$\boldsymbol {u}_0$. In the long-time limit the power iteration will converge to the asymptotic state of the linearised system, that is, the least stable/most unstable eigenmode of (2.11). More eigen-information cannot be extracted from the simple power method. To remedy this, the classical implicitly restarted Arnoldi method (IRAM) (Radke Reference Radke1996; Lehoucq et al. Reference Lehoucq, Sorensen and Yang1998) based on Krylov subspaces will be used to reduce the order of the original matrix, forming the much smaller Hessenberg matrix, to obtain the Ritz eigenmodes approximating the leading eigenmodes of the full system.
3. Validation
3.1. Validation of nonlinear DNS
In this section we will first show the evidence of converged calculations using Nek5000 and explain the choice of the grid resolution. Eight sets of grids from G1 to G8 in table 1 have been used to analyse the influence of spatial resolution on the DNS results at 
$Re=300$ and 
${\small \text{AR}}=1$. The coefficients 
$\bar {C}_d$ and 
$\bar {C}_{ly}$ and the 
$St$ number are used as evaluation indicators. Terms 
$N_{x1}$ and 
$N_{x2}$ are the numbers of the nodes in the 
$x$ direction as shown in the right panel of figure 1(
$b$); 
$N_y$ and 
$N_z$ are the node numbers for the sections related to 
$L_y$ and 
$L_z$, respectively. As previously mentioned, the polynomial order 
$N=7$ for all 
$Re$, which means that the number of the Gauss–Legendre–Lobatto nodes inside each hexahedral element is 
$8\times 8\times 8$ (in three directions). Since the flow in the boundary layer of the cylinder is sheared most, we tested the layer number of the O-type mesh 
$N_b=2, 4, 6, 8$, while keeping the total thickness of the O-type mesh constant. We can see from table 1 that, comparing G5–G7, when 
$N_b\geq 6$, the errors in 
$\bar {C}_d$ and 
$\bar {C}_{ly}$ are generally small. The error in 
$St$ (which appears more sensitive to the grid number) is also less than 1 %. The corresponding smallest boundary layer grid thickness is approximately 0.003 (a dimensionless value with the diameter 
$D$ as the reference length), which is smaller than the estimated value of 0.016 for 
$Re=1000$ in Tomboulides & Orszag (Reference Tomboulides and Orszag2000) who also used the SEM. Next, keeping the parameters of the boundary layer mesh unchanged, we test the effect of the element size in the wake area on the results. Comparing cases G4, G6 and G8, we can see that if 
$(N_{x1}+N_{x2})\times {N_y}\times {N_z}\geq (6+10)\times {6}\times {6}$, the differences in 
$\bar {C}_d$, 
$\bar {C}_{ly}$, 
$St$ among G4, G6 and G8 are also very small. Besides comparing cases G4 and Inoue & Sakuragi (Reference Inoue and Sakuragi2008)'s DNS at 
${\small \text{AR}}=1, Re=300$, the difference of 
$\bar {C}_d$ is approximately 0.53 %; but the difference in 
$St$ is 10.5 % (our 
$St=0.142$ in table 3 using mesh S2 and the DNS results in Inoue & Sakuragi (Reference Inoue and Sakuragi2008) is 
$St = 0.127$), and the reasons for this large difference will be discussed in § 4.3 to follow. Furthermore, as the 
$Re=300$ in table 1 is relatively low, we have also tested the cases with 
$Re=1000$ using the meshes G4 and G8 for a turbulent wake flow. In this case, the time histories of 
$C_d$ and 
$C_{ly}$ are intermittent. We compare the averaged values of theirs within 400 dimensionless time units. As shown in table 2, the difference in 
$\bar {C}_d$ between the two grids is small and the values of 
$\bar {C}_{ly}$ by themselves are also small, close to zero. The latter is because the flow is turbulent and statistically, there is not a preferred direction of the lift force in 
$y$. With these results, we decided to use mesh G4 for 
$Re\leq 500$ and mesh G8 for 
$500\leq Re\leq 1000$ in our numerical simulations.
Table 1. A grid sensitivity test for the case 
$Re=300$, 
${\small \text{AR}}=1$, 
$\Delta t=10^{-3}$. Numbers 
$N_b, N_{x1}$ and 
$N_{x2}$ are shown in figure 1(
$b$). Terms 
$N_y/2$ and 
$N_{z}/2$ are the node numbers on the edges 
$L_y$ and 
$L_z$, respectively. The total number of hexahedral elements in the computational domain is 
$N_{total}$. The errors in the parentheses are computed using as the reference the results of mesh G8.
Table 2. Grid sensitivity test for a high 
$Re=1000$ with 
${\small \text{AR}}=1$ and 
$\Delta t=10^{-3}$.
Next, we deal with the size of the computational domain for G4. As mentioned previously (§ 2.1), the numerical structural stability analysis (Giannetti & Luchini Reference Giannetti and Luchini2007) of the flow around a 2-D cylinder indicates that the results of the stability analysis will change significantly if the computational domain is not sufficiently large and its boundaries are placed close to the flow core area. For a 2-D cylinder, the core area is the recirculation zone near the cylinder. In table 3 we tested four computational sizes for the G4 grid. The blockage ratios (
$=A/S_{in}$) for S0, S1, S2 and S3 are 4 %, 1 %, 0.25 % and 0.098 %, respectively. In order to ensure that the spatial resolution remains relatively the same, the element number outside the O-type mesh zone of mesh S3 increases accordingly with the increase of computational domain size, so its 
$N_{total}$ is slightly larger. We set the results of the S3 case as the reference. It can be seen that when the blockage ratio is less than 1 %, the differences in time-averaged drag and lift coefficients are less than 
$2\,\%$, and 
$St$ is relatively more sensitive. For the blocking ratio less than 0.25 %, the differences in 
$\bar {C}_d$, 
$\bar {C}_{ly}$, 
$St$ are less than 0.6 %. Meanwhile, we also tested the influence of the computational domain size on the eigenvalue solution based on the mean flow. The effect of the computational domain size on the eigenfrequency 
$\omega$ is similar to that on the 
$St$ number. On the other hand, the effect of the computational domain size on the growth rate is stronger than the eigenfrequency. This is because when 
$Re$ is greater than the critical Reynolds number, the linear stability results based on the mean flow are marginally stable (meaning that the growth rate is nearly zero), so any influence of the computational domain size on the growth rate will appear relatively large. These tests are consistent with the numerical tests in Giannetti & Luchini (Reference Giannetti and Luchini2007), whose results indicated that if the boundaries are not placed close to the core zone of instability, one can get reasonable eigenvalues.
Table 3. A sensitivity test for the computational domain size. Here 
$Re=300$, 
${\small \text{AR}}=1$, 
$\Delta t=10^{-3}$. The errors are computed using as a reference the results of the mesh S3.
After the numerical verification step, in the result section, we will use a computational domain for the case of 
${\small \text{AR}}=1$ with a blocking ratio of 0.25 %. For other 
${\small \text{AR}}$ values of the cylinder flow, the computational domain size is increased accordingly to keep the blocking ratio no greater than 0.375 %. For the cases with the same 
${\small \text{AR}}$, when 
$Re$ increases, the computational domain size remains unchanged. In this way, the effect of the computational domain size on the growth/decay rate 
$\sigma$ is fixed to analyse the influence of 
$Re$ on 
$\sigma$.
3.2. Global stability analysis and validation
In this section we will show the global stability analyses of the flow around the finite-length cylinder. Since in the present literature we are not able to find the global LSA results of the fully 3-D flow around a short cylinder with two free ends, we will validate our implementation of the global stability code in the 2-D cylinder flow. The Reynolds number range being explored is 
$35\leq Re\leq 500$ (
$Re$ here is also based on the cylinder's diameter 
$D$). In figure 2(
$a$,
$b$) we compare the real part of the leading eigenmode (in the form of vorticity) by the present LSA based on the time-mean flow with the leading eigenmode calculated by Barkley (Reference Barkley2006). The two results are visually the same. In figure 2(
$c$,
$d$) we further compare the leading eigenfrequency and the growth rate (by both base flow and mean flow) from our IRAM code with the 
$St$ results and growth rates obtained by previous stability analysis (Barkley Reference Barkley2006), DNS (Leontini et al. Reference Leontini, Thompson and Hourigan2010) and experiment (Williamson Reference Williamson1989) as a function of 
$Re$. One can see a good agreement of our results with theirs.
Figure 2. Comparison of eigenpairs between present IRAM code and the results in the literature for a 2-D cylinder wake flow. Panels (
$a$,
$b$) are the vorticity of the leading eigenmode by present code and Barkley (Reference Barkley2006) at 
$Re=100$, respectively. Panel 
$(c)$ is the eigenfrequency in IRAM and the 
$St$ in DNS (Leontini et al. Reference Leontini, Thompson and Hourigan2010) & experiment (Williamson Reference Williamson1989). Panel 
$(d)$ shows the growth rates calculated by the two methods.
In the infinite cylinder wake flow there are two critical values of 
$Re$: one is around 
$46$, across which the flow undergoes a Hopf bifurcation giving rise to the 2-D vortex shedding and the other is approximately 
$180$, across which the 
$St$–
$Re$ curve shows discontinuity (Williamson Reference Williamson1989), which is due to the 3-D nature of the wake flow. Comparing the results in figure 2 regarding the critical Reynolds number, one can see that when 
$Re<46.1$, the eigenvalues obtained using the base flow and the mean flow are the same. This means that for the physical flows without oscillation or with decaying oscillation, stability analysis based on the base flow or the mean flow will yield the same result, which also verifies our implementation of the two methods (one is SFD-based LSA and the other is time-average DNS results). When 
$46.1< Re<180$, we can observe that stability analysis based on the mean flow can reproduce the experimental results but that based on the SFD base flow cannot, which has been observed and discussed in detail by Pier (Reference Pier2002) and Barkley (Reference Barkley2006) among many others. When 
$Re>180$, the result of LSA can no longer reflect the true vortex shedding frequency in the experiment (Williamson Reference Williamson1989), but its eigenfrequency is still close to the vortex shedding frequency obtained by the 2-D DNS (Leontini et al. Reference Leontini, Thompson and Hourigan2010). As previously mentioned, the discontinuous change of the 
$St$ number is caused by the 3-D nature of the wake, so 2-D DNS and LSA are no longer suitable for predicting the vortex shedding frequency in experiments (dot-dashed lines) when 
$Re>180$.
As a further validation step of the 3-D results, we use the power method to solve the linear equation (2.7) (linearised around the mean flow), then compare its growth rate and frequency with the leading eigenvalue obtained by IRAM. The power method can yield the eigenvalue of the largest absolute amplitude of a matrix, provided it is well separated from the second one. As shown in figure 3(
$a$), the leading eigenvalues obtained by the power method and IRAM are very close to each other at 
${\small \text{AR}}=1$, 
$270< Re<350$. In figure 3(
$b$) the growth rates at different 
$Re$ are obtained by the linear interpolation of the DNS results, the power method and the IRAM. The value of the critical 
$Re_{c1}$ obtained by the linear interpolation (of the DNS results) is almost the same as that determined by the IRAM, which are 
$172.2$ and 
$171.4$, respectively (the difference between them is 
$0.46\,\%$). For the cases far away from 
$Re_{c1}$, the growth rates obtained by these two methods appear different (due to the nonlinearity). The critical Reynolds numbers 
$Re_{c2}$ obtained by the power method and the IRAM are also close to each other.
Figure 3. The leading eigenvalues obtained by various methods at 
${\small \text{AR}}=1$, 
$270\leq Re\leq 350$. (
$a$) Growth rates 
$\sigma$ (imaginary part) and eigenfrequencies 
$\omega$ (real part) for different 
$Re$. Each data point is at a certain 
$Re$. (
$b$) Growth rates as function of 
$Re$.
4. Result and discussions
In this section we will present detailed results of global stability analysis and DNS of the flows past a finite-length cylinder. We will show the base states in a nonlinear system in the range of 
$Re\leq 1000$ and 
${\small \text{AR}} \leq 2$, which is characterised by time-averaged aerodynamic coefficients (i.e. 
$C_d, C_l$ etc.). We then study the flow transitions from steady wake pattern P1 to chaotic pattern P4 (see § 4.2) for 
${\small \text{AR}}=1$. This is followed by a discussion of the global modes and how they are connected to the monochromatic flow patterns P3. In the end, we will discuss the effect of 
${\small \text{AR}}$ with the significance of our results summarised in figure 21.
4.1. Base states in nonlinear system
In this section we will show the results of the nonlinear DNS for the flow around the finite-length cylinder. First, we observe that with the increase of 
$Re$, the wake pattern will change, which is reflected in the structure of the mean flow field. The drag and lift coefficients in the mean flow are shown in figure 4. We can fit the relation between 
$\bar {C}_d$ and the Reynolds number for the 
${\small \text{AR}}$ values considered in figure 4(
$a$) using the power function 
$\bar {C}_d=aRe^b+c$. The fitting coefficients 
$(a, b, c)$ are shown in the legend of figure 4(
$a$). The goodness of fit 
$R^2$ for each 
${\small \text{AR}}$ is greater than 0.99. On the other hand, for the lift coefficient, the change in wake pattern has a significant effect on 
$\bar {C}_{ly}$ as shown in figure 4(
$b$). For example, at 
${\small \text{AR}}=1$, each change in the wake pattern by increasing 
$Re$ will cause discontinuity or kinks in 
$\bar {C}_{ly}$, as shown at 
$Re=172, 282$ and 
$550$. These discontinuous or non-smoothness positions correspond to bifurcation points, as we will explain in the following global linear stability analyses.
Figure 4. (
$a$) Time-averaged drag coefficient 
$\bar {C}_d$ as a function of the Reynolds number. The dotted lines are fitted curves by the functions 
$\bar {C}_d =aRe^b+c$, and the fitting coefficients 
$(a, b, c)$ are shown in the parentheses of the legend. (
$b$) Time-averaged lift coefficient 
$\bar {C}_{ly}$ as a function of the Reynolds number. The sign of 
$\bar {C}_{ly}$ is not important because of the symmetry of the flow; here we choose to present it in negative values.
Comparing the values of 
$\bar {C}_{ly}$ of 
${\small \text{AR}}=0.5, 0.75, 1$ and 
$1.5$, we can observe that when 
${\small \text{AR}}$ becomes larger, the value of 
$Re$ at which 
$\bar {C}_{ly}$ becomes non-zero decreases. For these 
${\small \text{AR}}$ values, 
$\bar {C}_{ly}$ decreases to a negative value from zero and then increases and returns to zero. This trend of the 
$\bar {C}_{ly}$ result with increasing 
$Re$ indicates that there is a spatial asymmetry generated in the flow evolving from the laminar solution when 
$Re$ becomes larger and this asymmetry disappears (so 
$C_{ly}\rightarrow 0$) when 
$Re$ becomes even larger. At the onset of non-zero 
$\bar {C}_{ly}$, the gradient of 
$\bar {C}_{ly}$ is quite steep, signalling a sudden transition of the flow states. On the other hand, when 
${\small \text{AR}}$ becomes large, the 
$\bar {C}_{ly}$ results look very different: we can see from panel 
$(b)$ that when 
${\small \text{AR}}$ is large (e.g. 
${\small \text{AR}}=2$), the value of 
$\bar {C}_{ly}$ is always approximately zero for 
$Re$ from 0 to 300 as we investigated here, in which region we can find non-zero values of 
$\bar {C}_{ly}$ for small 
${\small \text{AR}}$ flows. This means that the value of 
${\small \text{AR}}$ is fundamentally important in determining the dynamics of the flow past a finite-length cylinder. We have also superposed the results of sphere (Johnson & Patel Reference Johnson and Patel1999) and interestingly, we find that the 
$\bar {C}_{ly}$ result of 
${\small \text{AR}}=0.75$ in our flow is very close to that of the sphere. Besides, this similarity also exists in the characteristic size of the separation bubble as shown in figure 6(
$a$). Thus, it is reasonable to believe that the finite cylinder with a small 
${\small \text{AR}}$ may behave dynamically similarly to the sphere. As a side note, Inoue & Sakuragi (Reference Inoue and Sakuragi2008) did not report the non-zero value of 
$\bar {C}_{ly}$ in their DNS results, but this phenomenon exists in our results and has also been confirmed in the flow past a sphere (Johnson & Patel Reference Johnson and Patel1999; Tomboulides & Orszag Reference Tomboulides and Orszag2000) and axial flow past a short cylinder (Pierson et al. Reference Pierson, Auguste, Hammouti and Wachs2019).
4.2. Flows past a finite cylinder at 
${\small \text{AR}}=1$
The main parameters in this work are 
$Re$ and 
${\small \text{AR}}$. In order to give a clear presentation of our results, in this section we will mainly use the case 
${\small \text{AR}}=1$ as a prototype. In § 4.4 we will present the effect of changing 
${\small \text{AR}}$.
4.2.1. Steady flows
First, let us discuss the steady flows that one can obtain at different values of 
$Re$. For the cylinder with 
${\small \text{AR}}=1$ and 
$Re$ being larger than approximately 
$10$, the flow begins to separate on the ends of the cylinder, a closed recirculation zone begins to form, and a pair of counter-rotating vortices appear in the wake. At this time, there are two mutually perpendicular symmetric planes (
$x$–
$y$ and 
$x$–
$z$ passing through the origin of the coordinate system) in the wake, and the centroid of the cylinder is on these two planes. We call this flow pattern P1. The corresponding Reynolds number at the onset of pattern P1 is 
$Re_{c0}$. In general, the 
${\small \text{AR}}$ is inversely correlated with 
$Re_{c0}$. The smaller 
${\small \text{AR}}$ is, the closer 
$Re_{c0}$ is to 20 in the case of the sphere wake (Johnson & Patel Reference Johnson and Patel1999), and the larger 
${\small \text{AR}}$ is, the closer 
$Re_{c0}$ is to 4 in the case of infinite-length 2-D cylinder flow. Specifically, the critical Reynolds numbers of the finite cylinder with 
${\small \text{AR}}=0.5, 0.75, 1.0, 1.5$ and 
$2$ are 
$Re_{c0}\approx 18, 15, 10, 8$ and 
$6$, respectively.
A streamline diagram is used to characterise the typical wake structure of the pattern P1 for 
${\small \text{AR}}=1$ and 
$Re=100$ in figure 5. The specific cross-section is the 
$x$–
$y$ and 
$x$–
$z$ planes of the flow. Unless stated otherwise, the streamwise direction of the flow field diagram in this article is from left to right. The topological traits of the pattern P1, such as the separation position on the cylinder's end (
$\theta _s$), the vortex centre point (
$x_b,y_b$) or (
$x_b,z_b$) and the length of the separation bubble (
$x_s$), have also been shown in figure 5. The fluid separates from the cylinder's end at the angle 
$\theta _s$ (the vertex of the angle is at the origin, one edge ends at the forward stagnation point, and the other at the separation point), and converges at a point 
$x_s$ on the 
$ox$-axis, forming a closed separation bubble in the 
$x$–
$z$ plane. The coordinates of the centre of the separation vortex pair in the 
$x$–
$y$ and 
$x$–
$z$ planes are 
$(x_b,y_b)$ and 
$(x_b,z_b)$, respectively, and the velocity at the vortex centre point (
$x_b,y_b,z_b$) is zero. All the lengths are non-dimensionalised by the cylinder diameter. With these parameters, we can quantify the effect of 
$Re$ on the pattern P1, as follows.
Figure 5. A typical steady wake pattern P1 with two symmetrical planes perpendicular to each other, at 
$Re=100$ and 
${\small \text{AR}}=1$. (
$a$) Streamlines on plane 
$x$–
$y$, (
$b$) streamlines on 
$x$–
$z$ plane.
Previous works on flow past a sphere (Fornberg Reference Fornberg1988; Johnson & Patel Reference Johnson and Patel1999; Tomboulides & Orszag Reference Tomboulides and Orszag2000) found and confirmed that when 
$Re>75$, the separation angle and the separation bubble length are logarithmically related to the Reynolds number. Following this idea, we found that a logarithmic relation also holds in the flow past short cylinders, as illustrated in figure 6. For the steady flow, the separation bubble length can be approximately related to 
$Re$ by 
$x_s=a\ln (Re+b)+c$ (see the dotted curves in panel (
$a$) and for the unsteady flow, the separation angle and separation bubble length are measured in the mean flow). The goodness of fit 
$R^2$ for 
${\small \text{AR}}=1$ is 0.9988, and exceeds 
$0.999$ for the other four cylinders. The turning points of curve 
$x_s$–
$Re$ in figure 6(
$a$) correspond to the transition of the wake mode (from steady flows to time-periodic flows), which also can be observed in 
$\bar {C}_{ly}$–
$Re$ (figure 4
$b$). One can see that the 
$x_s$–
$Re$ relation curve for the flow past a sphere is close to the 
${\small \text{AR}}=0.75$ data. As pointed out by one of the referees, this aspect ratio is notable because, interestingly, it has a projected frontal area that is very close to that of a sphere (i.e. 
${\rm \pi} /4\approx 0.785$), which can also explain that the curve 
$\bar {C}_{ly}$–
$Re$ (
${\small \text{AR}}=0.75$) in figure 4(
$b$) almost coincides with that of the sphere.
Figure 6. (
$a$) Separation bubble length and (
$b$) separation angle as a function of the Reynolds number. The dotted lines represent fitting functions 
$x_s=a\ln (Re+b)+c$ and 
$\theta _s =aRe^b+c$ (see figure 5(
$a$) for the schematic definitions). The fitting coefficients 
$(a, b, c)$ are shown in the legends. In panel (
$b$) the black, red and green dashed lines are the averaged value of the separation angles in the 
$x$–
$y$ and 
$x$–
$z$ planes for the short cylinders with 
${\small \text{AR}}=0.5,0.75$ and 
$1$, respectively.
In figure 6
$(b)$ the relationship between the separation angle and the Reynolds number is approximately a power function 
$\theta _s=aRe^b+c$ within the range of 
$Re<600$. It turns out that the separation angle is less sensitive to the transition, and almost no visible kinks are found in these results. At a first glance, the results of the sphere from Johnson & Patel (Reference Johnson and Patel1999) deviate significantly from the result of 
${\small \text{AR}}=0.75$ (note that these results in figure 6(
$b$) are only for 
$\theta _s$ in the 
$x$–
$y$ plane). A more fair comparison entails an additional consideration of the separation angle in the 
$x$–
$z$ plane. With reference to figure 5(
$b$) (plane 
$x$–
$z$), the separation angles at the sharp trailing corners of the short cylinder are fixed at 135 
$^\circ$ (
${\small \text{AR}}=1$), 
$143.13^\circ$ (
${\small \text{AR}}=0.75$) and 
$153.44^\circ$ (
${\small \text{AR}}=0.5$). Hence, we can simply average the two angles in the 
$x$–
$z$ and 
$x$–
$y$ planes and plot the new results as dashed lines in figure 6(
$b$). One can see that in this case, the results of short cylinders (for 
${\small \text{AR}}=0.5,0.75,1$) indeed straddle that of the sphere, consistent with the previous simple argument as suggested by the referee.
When 
$Re$ further increases, larger than a certain critical 
$Re_{c1}$, the current DNS results show that the pattern P1 will be unstable and transition to another steady wake pattern P2 (figure 7) with only one symmetrical plane (
$x$–
$y$), simultaneously inducing non-zero values of 
$\bar {C}_{ly}$. The bifurcation of the short cylinder at 
$Re_{c1}$ is steady, time-independent without vortex shedding (
$St=0$), so the flow experiences a pitchfork bifurcation. It is noted that the pattern P2 may not appear for certain values of 
${\small \text{AR}}$. As we will show in the end of the next section, 
${\small \text{AR}}=1.75$ is such an example. We will present more systematic results of the effect of 
${\small \text{AR}}$ in § 4.4. This is similar to the flow past a sphere, where there is a regular bifurcation before the Hopf bifurcation. The difference is that the wake of the sphere retains a symmetric plane that is randomly oriented, which can be related to the numerical uncertainty or the perturbations in the experimental environment (Tomboulides & Orszag Reference Tomboulides and Orszag2000). In the literature there is no discussion about the P2-type flow pattern (Inoue & Sakuragi Reference Inoue and Sakuragi2008 did not report this) or its associated regular bifurcation for the flow past a short cylinder with its axis being perpendicular to the streamwise direction (in the following, for the ease of discussion, we will call this flow a radial flow; likewise, we call the flow past a short cylinder with its axis parallel to the streamwise direction an axial flow). Recently, Pierson et al. (Reference Pierson, Auguste, Hammouti and Wachs2019) conducted a numerical simulation of an axial flow past a short cylinder at 
${\small \text{AR}}=1, 20\leq Re\leq 420$ and 
${\small \text{AR}}=3, 25\leq Re\leq 250$. They found a regular bifurcation for the axial flow around the cylinder at 
${\small \text{AR}}=1$ and 
$Re\approx 278$ (based on the cylinder diameter) without hysteresis, but this bifurcation was not found for the cylinder with 
${\small \text{AR}}=3$. Sheard et al. (Reference Sheard, Thompson and Hourigan2008) found that there is a regular bifurcation for the radial flow past a short cylinder with hemispherical ends at 
$Re=350\pm 2$. Finally, 3-D transitions of the bluff-body cube have been numerically studied by Saha (Reference Saha2004). The cube wake also undergoes a regular bifurcation (losing one symmetric plane but retaining time independence) at 
$Re=216$ and a Hopf bifurcation at 
$Re=270$. From these works, it can be hypothesised that the aspect ratio 
${\small \text{AR}}$ of the bluff body (including the short cylinder with flat ends as we study here) is related to the existence of a regular bifurcation.
Figure 7. A typical P2 type wake at 
$Re=278$ and 
${\small \text{AR}}=1$, (
$a$) 
$x$–
$y$ plane, (
$b$) 
$x$–
$z$ plane. Vorticity (
$x$-component) isosurface (60 % transparency) is coloured by values 
$-0.1$ and 
$0.1$. The number 22 is the length of the scale ruler. The green object is the cylinder.
Now we present the way we determine 
$Re_{c1}$. We perform nonlinear DNS for 
${\small \text{AR}}=1$, and the critical value 
$Re_{c1}$ is approached by gradually increasing the Reynolds number from a small value. For each 
$Re$, we can calculate the growth rate of disturbance development as in figure 8. In this paragraph we focus on the blue curves for 
$Re=180$ without randomness in the initial condition and will discuss the other curves shortly. The time series of 
$C_{ly}$, when plotted in a logarithmic scale against the time in a linear scale (as in panel 
$b$), will display a period of a straight line, which represents the linear phase. The linear growth or decay rate can then be calculated and the critical 
$Re$ corresponding to the zero growth rate can be interpolated by fitting several data points using a linear function, which is shown in figure 10, to be discussed shortly. For 
${\small \text{AR}}$=1, we can see that 
$Re_{c1}=172.2$, which is consistent with the 
$Re$ of onset of non-zero 
$\bar {C}_{ly}$ shown in figure 4(
$b$). The 
$Re_{c1}$ values of other 
${\small \text{AR}}$ have also been calculated and it is found that when 
${\small \text{AR}}$ is small or large, 
$Re_{c1}$ is relatively large. The 
${\small \text{AR}}$ having the minimum 
$Re_{c1}$ is approximately located around 
$1.25\sim 1.5$. Besides, the 
${\small \text{AR}}$-dependence of the critical Reynolds number will be revealed in § 4.4. In the analysis of Tomboulides & Orszag (Reference Tomboulides and Orszag2000), the authors also performed a logarithmic transformation of the azimuthal velocity at a point located in the wake and near the sphere, and found that the logarithm of the azimuthal velocity is a linear function of time during the initial evolution. The growth rate obtained by this method is consistent with the LSA result of Natarajan & Acrivos (Reference Natarajan and Acrivos1993).
Figure 8. Time histories of lift coefficient 
$C_{ly}$ and 
$\ln (|C_{ly}|)$. 
$(1+0.2\epsilon _{rand})$ in the legend means that the initial condition is the base flow 
$(1,0,0)$ plus a white noise of 0.2 times the base flow in the three velocity components. Likewise for the other initial conditions.
Some numerical experiments have also been conducted to perturb the flows in figure 8. We applied white noise (denoted as 
$\epsilon _{rand}$) to the initial conditions in the DNS at 
${\small \text{AR}}=1, Re=165, 168, 170$ (subcritical, as 
$Re< Re_{c1}=172.2$). The random numbers are shifted to perturb the solution positively and negatively. The amplitude of the white noise is 0.2 times the inlet velocity for the three cases and additionally, for 
$Re=168$, we also examined the amplitudes of the (white-noise) disturbance being 0 and 0.6 times the inlet velocity. The white noise is applied to all the three velocity components through the following initial conditions ((4.1) and (4.2)):
\begin{gather} \left. \begin{array}{ll@{}} U_x = 1.0 +a\epsilon_{rand}, & (a_1=3.0\times10^4, a_2={-}1.5\times10^3, a_3 = 0.5\times10^5),\\ U_y = 0.0 +a\epsilon_{rand}, & (a_1=2.3\times10^4, a_2=2.3\times10^3, a_3 ={-}2.0\times10^5),\\ U_z = 0.0 +a\epsilon_{rand}, & (a_1=2.0\times10^4, a_2=1.0\times10^3, a_3 = 1.0\times10^5), \end{array} \right\} \end{gather}
\begin{gather} \left. \begin{array}{l@{}} \epsilon_1 = a_1 (i_g+x \sin y) + a_2 i_x i_y + a_3 i_x, \\ \epsilon_2 = a_1 (i_g + z \sin\epsilon_1) + a_2 i_z i_x + a_3 i_z, \\ \epsilon_{rand} = \cos(10^3\sin(10^3\sin\epsilon_2)), \end{array} \right\} \end{gather}
Here 
$a$ is the maximum amplitude of the perturbation, 
$\epsilon _{rand}$ denotes white noise with random numbers between 
$-1$ and 1 in (4.2). 
$x,y,z$ are the coordinate values of the element nodes and 
$i_x,i_y,i_z,i_g$ are the numbers of the element nodes. From figure 8 we can see that the lift coefficients 
$C_{ly}$ for all the cases aforementioned eventually converge to zero (because the flows are subcritical) and the linear phase survives for a long time. This demonstrates the robustness of the results. A by-product of these numerical experiments is that they seem to indicate that the first bifurcation of the short cylinder around 
$Re_{c1}$ is probably supercritical for the parameters we consider here, because if the bifurcation was subcritical, the white noise with a large amplitude will probably bring the flow to a stable finite-amplitude nonlinear solution. In order to confirm the supercriticality of the first bifurcation in this flow, we estimate the coefficients of the Landau model from the DNS data (as is well known, the Landau coefficient calculated around the linear critical conditions determines the nature of the bifurcation; see Guckenheimer & Holmes Reference Guckenheimer and Holmes1983). The global variable 
$C_{ly}$ that can directly reflect the transition from P1 to P2 is used to evaluate the Landau coefficient at 
$Re=180$ and 
${\small \text{AR}}=1$. The original time history of 
$C_{ly}$ is shown in figure 8(
$a$), where 
$|C_{ly}|$ grows and finally saturates. Now, based on the Landau equation 
${\textrm {d}A}/{\textrm {d}t}=c_1 A+ c_3 |A|^2$ (where 
$A$ can be viewed as the amplitude of 
$C_{ly}$) or, equivalently, 
$\textrm {d}\log |A|/\textrm {d}t=c_1+ c_3|A|^2$, the Landau coefficients 
$c_1$ and 
$c_3$ can be calculated by plotting 
$\textrm {d}(\log |C_{ly}|)/\textrm {d}t$ vs 
$|C_{ly}|^2$ (Thompson et al. Reference Thompson, Leweke and Provansal2001), as shown in figure 9. Therefore, the vertical intercept point gives an estimation of 
$c_1=0.0201$ and the gradient near this point is an approximation of 
$c_3=-8.78$. The value of growth rate 
$c_1=0.0201$ is close to the linear growth rate 
$0.0196$ (at 
$Re=180$) obtained by calculating the linear slope in figure 8(
$b$). Furthermore, the negative 
$c_3=-8.78$ indicates that the (lowest-order) nonlinearity stabilises the flow. In panel (
$b$) one can also see that the nonlinearity stabilises the flow after the linear regime. Thus, the bifurcation is supercritical transitioning from P1 to P2. We mention in passing that, theoretically, using this method to infer the bifurcation type of a flow is most accurate when the parameter is close to the critical condition. Nevertheless, when this is satisfied, the evolution time of the flow is very long in DNS. As long as the bifurcation type is concerned, the method can be applied slightly away from the critical condition (Henderson & Barkley Reference Henderson and Barkley1996; Carmo, Meneghini & Sherwin Reference Carmo, Meneghini and Sherwin2010; Gao et al. Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013; Feng et al. Reference Feng, Zhang, Vazquez and Shu2021). Besides, we arrive at the above conclusion using only one instance of 
$Re$ and admit that a more rigorous analysis would be to apply the weakly nonlinear stability analysis around the linear critical conditions, which can be pursued in a future work.
Figure 9. Time derivative of log
$|A|$ vs 
$|A|^2$ with the vertical intercept giving the linear growth rate 
$c_1=-0.0201$ and the gradient near this intercept point giving the Landau coefficient 
$c_3=-8.78$, which indicates supercritical bifurcation behaviour. Nonlinear DNS data of case 
${\small \text{AR}}=1$, 
$Re=180$.
To sum up, the P2-type steady flow has been found for short cylinders with small 
${\small \text{AR}}$. Most of the relevant research on the flow past a short bluff body have focused on the sphere or cylinder with a large 
${\small \text{AR}}$, and no previous studies have systematically investigated this P2-type wake mode in the axial flow past a finite cylinder. Through our DNS results in the range of 
$0.5<{\small \text{AR}}<2$, it is found that the 
${\small \text{AR}}$ of the short cylinder has a very significant and important influence on the critical Reynolds number 
$Re_{c1}$. Only when the 
${\small \text{AR}}$ is relatively small (
${\small \text{AR}}< 1.75$) will there be a regular bifurcation. As can be seen in figure 10, when 
${\small \text{AR}}=1.75$, we cannot find the regular transition by either DNS or LSA methods. In fact, when 
${\small \text{AR}}\gtrsim 1.75$, as 
$Re$ increases, the wake of the radial flow around the short cylinder will not experience the P2-type wake, and directly transition from the P1-type to a periodic shedding of hairpin-shaped vortices P3-0 with 
$\bar {C}_{ly}=\bar {C}_{lz}=0$, which has been reported many times in the flow past a bluff body with 
${\small \text{AR}}>2$, for example, in the flow past a cylinder with two free hemispherical ends (Sheard et al. Reference Sheard, Thompson and Hourigan2005, Reference Sheard, Thompson and Hourigan2008), a cylinder with two free flat ends at moderate 
$2 \leq {\small \text{AR}} \leq 10$ (named type III by Inoue & Sakuragi Reference Inoue and Sakuragi2008) and at 
${\small \text{AR}}=3$ and 
$Re>125$ (Pierson et al. Reference Pierson, Auguste, Hammouti and Wachs2019's figure 29), a prolate spheroid at 
$L/D=6$ and 
$Re=100$ (El Khoury, Andersson & Pettersen Reference El Khoury, Andersson and Pettersen2012) and a freely falling cylinder following rectilinear paths at Archimedes number 
$Ar=200$ and 
$L/D=5$ (Toupoint, Ern & Roig Reference Toupoint, Ern and Roig2019's figure 17a). Nevertheless, we will report in the next section that there are two unreported periodic vortex shedding wakes (P3-1 and P3-2) of radial flow around the short cylinder with 
$\bar {C}_{ly}\neq 0$.
Figure 10. The growth rate (
$G_R$) calculated using the DNS results as a function of the Reynolds number. Unstable flows will transition from P1 pattern to P2 pattern via a regular bifurcation and the critical Reynolds number is 
$Re_{c1}$.
4.2.2. Vortex shedding
When 
$Re$ increases further, a second bifurcation appears at 
$Re\approx 282$ (for 
${\small \text{AR}}=1$) with 
$St\neq 0$, which belongs to a Hopf bifurcation. Figure 11 shows the relationship between the linear growth rate and 
$Re$ for short cylinders of different 
${\small \text{AR}}$ values from 
$0.5$ to 2 near the Hopf bifurcation point. The relationship between the growth rate and the Reynolds number is approximately linear when the flow is stable and near the critical condition. The solid lines in figure 11 represent linear fitting results which are shown in the legend. The general trend is that the larger the value of 
${\small \text{AR}}$ is, the smaller the critical value of 
$Re_{c2}$. The denoted eigenmodes A and B will be discussed in § 4.3.
Figure 11. The growth rate 
$\sigma$ of the leading eigenmode in the global LSA (of the time-mean flow) as a function of 
$Re$. Unstable flows will transition from P2 pattern to P3 pattern via a Hopf bifurcation and the critical Reynolds number is 
$Re_{c2}$.
We found that near the Hopf bifurcation point, one of the two periodic wake patterns P3-1 (a transitional state) or P3-2 (a saturated state) may appear and can exist for a long time. The P3-1 wake pattern is a hairpin-like vortex structure shedding from the side surface of the cylinder, as shown in figure 12(
$a$). The P3-1 type wake is symmetric about the 
$x$–
$y$ plane (note that the green ruler is collinear with the cylinder axis, that is, 
$z$), and the lift coefficient 
$C_{lz}$ in the 
$z$ direction is non-oscillating and approximately zero. Fast Fourier transform (FFT) applied to the drag and the lift coefficients gives the same and single-periodic vortex shedding frequency 
$St_{Re=283}=0.125$. Interestingly, we note that this value is comparable to the vortex frequency of the sphere 
$St_{Re=270}=0.1292$, 
$St_{Re=285}=0.1335$ obtained by Tomboulides & Orszag (Reference Tomboulides and Orszag2000) through nonlinear DNS. The cross-section of the P3-1 hairpin vortex on the 
$x$–
$y$ plane is asymmetric with respect to the 
$x$-axis, and the wake is offset in the 
$y$-axis (similar to the P2-type wake). This asymmetric periodic wake pattern is also similar to the asymmetric wakes of bluff bodies, for example, in the axial flow around a short cylinder studied by Pierson et al. (Reference Pierson, Auguste, Hammouti and Wachs2019) (who found that the axial flow becomes unsteady at 
$Re\approx 355\ {\rm and}\ {\small \text{AR}}=1$, but it still maintains a symmetric plane in the wake until 
$Re=420$), in the flow past a sphere by Tomboulides & Orszag (Reference Tomboulides and Orszag2000), Johnson & Patel (Reference Johnson and Patel1999) (who showed that the regular bifurcation at 
$Re=211$, and the asymmetric steady flow becomes unsteady at 
$Re=275$), and in the flow past a short cylinder with hemispherical ends by Sheard et al. (Reference Sheard, Thompson and Hourigan2008) (who found that the regular bifurcation at 
$Re=350\pm 2$). The other 3-D wake pattern P3-2 is shown in figure 12(
$b$). It can be seen that the vortex shedding from the ends begins to be dominant in the wake. The results of FFT applied on the aerodynamic coefficients show that the frequency of 
$C_d$ and 
$C_{ly}$ is the same, and the frequency of 
$C_{lz}$ is 1/2 that of 
$C_d$ and 
$C_{ly}$. Besides, the frequency of vortex shedding of the structure P3-2 is larger than that of the structure P3-1; for example, at 
$Re=290$, we calculated that 
$St_{P3-2}=0.1395$ whereas 
$St_{P3-1}=0.1269$. More results of the 
$St$ data for other values of 
$Re$ are shown in figure 19 to be discussed shortly.
Figure 12. Two typical periodic wake patterns for the case 
$Re=290$, 
${\small \text{AR}}=1$. (
$a$) wake pattern P3-1 stimulated by the initial condition 
$1+0.0\epsilon _{rand}$; (
$b$) wake pattern P3-2 stimulated by the initial condition 
$1+0.2\epsilon _{rand}$. The 
$Q$-criterion isosurfaces 
$Q=0$ are coloured by the 
$x$-component of the vorticity ranging from 
$-10^{-2}$ to 
$10^{-2}$. The green ruler is collinear with the cylinder axis.
In a subcritical flow (with 
$Re<282$ for 
${\small \text{AR}}=1$), both P3-1 and P3-2 patterns are stable and will converge to the P2 pattern. This is demonstrated in figure 13 using DNS where we have used different initial conditions to trigger the appearance of the two P3 patterns. We found that if one uses initial perturbation with a smaller amplitude at 
$Re=278$, the flow will more likely develop to the P3-1 structure; whereas, if a large disturbance is adopted, P3-2 is more likely to appear. To explain, the initial condition is the base flow 
$(U_x, U_y, U_z)$, e.g. (
$1, 0, 0)$, plus white noise with a maximum amplitude of 
$a=0.0, 0.2, 0.5$, denoted as 
$1+a\epsilon _{rand}$. Let us first claim the results that 
$1+a\epsilon _{rand}$ leads to the P3-1 structure, and 
$3+0.2\epsilon _{rand}$ gives rise to the P3-2 structure. The phase trajectory diagram of the lift-drag coefficient is shown in figure 13. For all the initial conditions considered, the 
$C_d$–
$C_{ly}$ curves of P3-1 and P3-2 patterns, as converging spirals, shrink to the same fixed point, that is, the attractor of the subcritical cases is a fixed point, which is the steady flow pattern P2 as shown in figure 7. On the other hand, the oscillation range of 
$C_{lz}$ corresponding to the initial conditions 
$1+a\epsilon _{rand}$ is much smaller than 
$C_{ly}$ and is close to zero, as shown in figure 13(
$b$). Based on this result, we think that the initial conditions 
$1+a\epsilon _{rand}$ leads to the P3-1 pattern. The phase trajectory of 
$C_d$–
$C_{lz}$ corresponding to the initial condition 
$3+0.2\epsilon _{rand}$ is obviously different from that of 
$1+a\epsilon _{rand}$ (the P3-1 pattern) and corresponds to the P3-2 structure.
Figure 13. The diagram of 
$C_d$–
$C_{l}$ at 
$Re=278$ and 
${\small \text{AR}}=1$. The expression 
$1+a \epsilon _{rand}$ in the legends indicates that the initial condition for the velocity is 
$(1,0,0)$ plus white noise with an amplitude of 
$a$. The white noise is three dimensional as seen in (4.2).
In a supercritical flow (with 
$Re>282$ for 
${\small \text{AR}}=1$), the P3-1 pattern is transient and will eventually converge to the saturated P3-2 pattern. However, one may be easily fooled to believe that the P3-1 pattern can be a stable mode because when 
$Re$ is close to the critical condition, the lingering time of the flow around the neighbourhood of P3-1 can be extremely long (before it saturates to the nonlinearly stable P3-2 state). This is shown in figures 14 and 15 for 
$Re=290$ and 
$300$, respectively. As one can see from figure 15 for 
$Re=300$, the two wake patterns and their transition are clear in a relatively short time. When 
$Re=290$, however, the DNS method becomes quite difficult to differentiate the two because one needs to simulate for a very long time to see the transition. In the above two figures, because a high-
$Re$ is more sensitive to numerical error, it makes the 
$Re=300$ flow transition to the saturated state more quickly. Another way to expedite the convergence to the saturated state is to add perturbation in the initial condition, resulting in a high possibility of throwing the flow out of the neighbourhood of P3-1 and converging to the final saturated state more quickly (results not shown).
Figure 14. The time histories of aerodynamic coefficients for cylinder 
${\small \text{AR}}=1$ at 
$Re=290$. The initial condition is the wake P3-1 of a lower-
$Re$ case.
Figure 15. The time histories of aerodynamic coefficients for cylinder 
${\small \text{AR}}=1$ at 
$Re=300$. The initial condition is 
$1+0.0\epsilon _{rand}$.
Some discussions are in order regarding the possible bifurcation type around 
$Re_{c2}$. In the subcritical flows our numerical simulations with different initial conditions cannot converge to a stable flow state other than P2 (which is a linearly stable state). In the supercritical flows, even though the lingering time around P3-1 can be very long, the final saturation of the flow is the P3-2 pattern; different initial conditions will only change the time spent by the flow around P3-1 before converging to the P3-2 wake. With these results, we infer that the nature of the bifurcation around 
$Re_{c2}$ should be supercritical. The Landau model has also been used to study the bifurcation type (following a similar analysis presented in figure 9). Again, 
$c_1=0.0272$ obtained in the Landau model is close to the global stability analysis results based on the SFD base flow (whose linear growth rate is 
$0.02746$). The negative value of 
$c_3$ indicates that the flow bifurcation may be supercritical, which can also be seen in the stabilising effect of the (lowest-order) nonlinearity in figure 16(
$b$). More analyses of the P3-1 and P3-2 patterns and the flow bifurcation will be conducted in § 4.3 on the global modes.
Figure 16. Time derivative of log
$|A|$ vs 
$|A|^2$ with the vertical intercept giving the linear growth rate 
$c_1=0.0272$ and the gradient near this intercept point giving the Landau coefficient 
$c_3=-254.1$, which indicates supercritical bifurcation behaviour. Nonlinear DNS data of case 
${\small \text{AR}}=1$, 
$Re=300$.
4.2.3. Approach to chaos
When 
$Re$ is further increased, the flow becomes more chaotic. Figure 17 shows 
$C_d$–
$C_{lz}$ phase diagrams and power spectral density (PSD) of the 
$C_{lz}$ data for four values of 
$Re$ in 
$338 \leq Re \leq 1000$. When 
$Re$ is larger than 
$332$, a new lower frequency 
$St=0.072$ starts to emerge in the flow field. This can be seen in panels (
$a{,}e$) for 
$Re=338$ where a small peak begins to form around 
$St=0.072$ in the PSD figure. This frequency becomes incommensurable when 
$Re$ further increases and the wake will gradually become chaotic, which is called pattern P4-0 (with 
$\bar {C}_{ly}=\bar {C}_{lz}=0$) and P4-1 (with 
$\bar {C}_{ly}\neq 0$, occurs at 
${\small \text{AR}}<1.75$) in this article. In Tomboulides & Orszag (Reference Tomboulides and Orszag2000)'s DNS study on the sphere, a similar lower frequency value was detected as 
$St_{Re=500}=0.045$ at a relatively large 
$Re$, which was the second highest peak in the PSD of their DNS data. The frequency corresponding to the dominant peak was still 
$St_{Re=500}=0.167$, which was dominant at all positions of the wake. In the work of the axial flow past a short cylinder, Pierson et al. (Reference Pierson, Auguste, Hammouti and Wachs2019) also found a similar low frequency (
$St_{Re=395}=0.03$) when 
$Re \gtrsim 395$ in their DNS study of the axial flow past a short cylinder.
Figure 17. The diagrams of 
$C_{d}$–
$C_{lz}$ in (a–d) and the power spectral density for the time history of 
$C_{lz}$ in (e–f) with 
${\small \text{AR}}=1$, at 
$Re=338$ (
$a$,
$e$), 
$Re=344$ (
$b$,
$f$), 
$Re=410$ (
$c$,
$g$) and 
$Re=1000$ (
$d$,
$h$).
Back to figure 17, at 
$Re=344$, a third new frequency signal 
$St=0.0083$ starts to appear. At this 
$Re$, the attractor of 
$C_d$–
$C_{ly}$ is no longer a limit cycle, but becomes a limit torus. It means that a certain two frequencies in the flow field are incommensurable, and the phase trajectory is no longer closed. The PSD of these two lower frequencies 
$St=0.0706$ and 
$St=0.0083$ is not negligible compared with the dominant frequency 
$St=0.1528$. At 
$Re=410$, the phase trajectory of 
$C_d$–
$C_{ly}$ becomes highly disordered, and the PSD curve shows the characteristics of a chaotic signal. It can be seen from figure 17(
$g$) that the dominant frequency is still around 
$St=0.1611$, and its PSD curve shows that there are many small peaks around 
$St=0.033$, 
$St=0.094$ and 
$St=0.1611$. Finally, when 
$Re=1000$, the phase diagram is highly chaotic and a broad spectrum of frequencies appears signalling the disordered state of the flow. In the literature, Sakamoto & Haniu (Reference Sakamoto and Haniu1995)'s experimental work on the flow past a sphere showed that, when 
$Re>420$, the shedding direction of the hairpin vortex appears intermittent, the oscillation amplitude and waveform caused by the vortex shedding start to become irregular, and the flow field starts to transition from a single-frequency flow to a chaotic state. Pierson et al. (Reference Pierson, Auguste, Hammouti and Wachs2019)'s DNS study (in its supplementary materials) on the axial flow past a finite cylinder (
${\small \text{AR}}=1$) shows that the flow exhibits chaotic properties for the 
$Re$ in 
$420 \leq Re \leq 460$, and the secondary lower frequency (
$St_{Re=460}=0.03$) gradually dominates the entire flow field. These results of the transition from a single frequency to a chaotic state is similar to that of the radial flow past a finite cylinder (
${\small \text{AR}}=1$) as we study here.
4.3. Global modes
The discussions in § 4.2 pertain to the different flow patterns (P1–P4) one can obtain from the nonlinear DNS when increasing the value of 
$Re$ and the critical values of 
$Re$ delimiting the different flow regimes. In this section we will look into the eigenmodes at certain values of 
$Re$ and 
${\small \text{AR}}$ for some base states. We report this because we found that when the parameters change, the type of mode that becomes the most unstable changes. If not stated otherwise, the global LSA to follow is performed with respect to the mean flow.
Based on the symmetrical and asymmetrical stationary base states, the global eigenmodes obtained by IRAM codes can be divided into two categories. The first type is the eigenmode with only one symmetric plane (
$x$–
$y$) as shown in figure 18, and it can be found at the cases 
$Re>Re_{c1}, {\small \text{AR}}<1.75$. The second type is the eigenmode with two mutually perpendicular symmetric planes (
$x$–
$y$ and 
$x$–
$z$) as shown in figure 20 in the parameter ranges 
$Re< Re_{c1}$ and 
${\small \text{AR}}<1.75$ or 
${\small \text{AR}}\geq 1.75$. It can be understood that the global linear eigenmode retains the same structural symmetry as the basic state (either base flow or mean flow), which is consistent with the results of the global LSA of the classical 3-D non-rotating and rotating sphere (Citro et al. Reference Citro, Tchoufag, Fabre, Giannetti and Luchini2016).
Figure 18. The first three eigenmodes at 
$Re=290$ and 
${\small \text{AR}}=1$. (
$a$) Eigenmode A with eigenvalue 
$\lambda =2.959\times 10^{-3}+0.1403\boldsymbol {i}$ (real part is the growth rate and imaginary part the eigenfrequency divided by 
$2{\rm \pi}$), (
$b$) eigenmode B with eigenvalue 
$-3.885\times 10^{-2}+0.1258\boldsymbol {i}$, (
$c$) eigenmode C with eigenvalue 
$-7.491\times 10^{-2}+0.04345\boldsymbol {i}$. The 
$Q$-criterion isosurfaces 
$Q=0$ are coloured by the 
$x$-component of the vorticity ranging from 
$-10^{-3}$ to 
$10^{-3}$.
We first discuss the case 
${\small \text{AR}}=1, Re=290$, which has the first type of eigenmodes. The typical flow structures of the first three eigenmodes in this case are shown in figure 18, which are ordered by the real part (the growth rate) of the eigenvalues from large to small, and are marked as eigenmodes A, B and C, respectively. The imaginary part of the eigenvalues of the first three eigenmodes are all not equal to zero. Modes A and B appear to have smaller flow structures than mode C, whose frequency is smaller. Interestingly, we found that the eigenmodes A and B can be favourably compared with the DNS results and are highly relevant to, respectively, the P3-2 and P3-1 structures that we have identified and discussed earlier. Following Sansica et al. (Reference Sansica, Robinet, Alizard and Goncalves2018), we have also checked that the mode shape and frequency in the dynamic mode decomposition (DMD following Schmid Reference Schmid2010) analysis of the nonlinear flow field snapshots of the P3-1 pattern are almost the same as those of the linear eigenmode B (DMD results are not shown). The same correspondence between the DMD mode of the P3-2 pattern and the global mode A can also be established. Based on these strong links, to further understand the Hopf bifurcation, the relationship between the eigenvalues of these two eigenmodes and the Reynolds number in the stability analyses based on the mean flow and the SFD base flow is shown in figure 19 (with the specific data in table 4). In the two panels, the vertical dash-dot line is 
$Re=Re_{c2}$, and the dotted horizontal line is the zero growth rate 
$\sigma =0$. In panel (
$b$) the dashed cyan and red solid lines are approximate fittings of the non-dimensional vortex frequency 
$St$ measured in the nonlinear DNS (the triangles). The red solid line also indicates stable solutions and the dashed cyan lines transient flow states.
Figure 19. The eigenvalues and Strouhal numbers as a function of 
$Re$ near the Hopf bifurcation at 
${\small \text{AR}}=1$. The red arrow ‘
$\uparrow$’ indicates that the flow will transition to the wake pattern P3-2 at 
$Re>Re_{c2}$ under the effect of disturbance. In panel (
$b$) the red solid line represents stable solutions and the dashed cyan lines represent transient flows. The specific data are shown in table 4.
Table 4. Comparisons between the most unstable eigenvalues and Strouhal numbers near the Hopf bifurcation for the cylinder 
${\small \text{AR}}=1$. The superscripts 
$^{ BF}$ and 
$^{ MF}$ represent SFD base flow and mean flow, respectively.
From panel 
$(a)$ one can see that at a subcritical 
$Re$, the mean flow and the SFD base flow generate the same results for modes A and B. The decay rates of mode B are both smaller than those of mode A (i.e. mode A is more stable), which is consistent with the result that it is sometimes possible to observe wake pattern P3-1 in nonlinear DNS in this range of 
$Re$ (note that mode B corresponds to the P3-1 pattern). In the supercritical range the most unstable mode is now eigenmode A (circles based on the SFD base flow and green filled squares based on the mean flow), whose 
$\sigma$ increases with 
$Re$. The eigenmode B (red crosses) based on the mean flow is stable. Thus, the amplitude of eigenmode A will increase exponentially with time until the nonlinearity becomes important. This can explain that in the evolution of nonlinear DNS, the frequency information of wake pattern P3-2 (corresponding to eigenmode A) is observed in a supercritical condition.
Now we look at panel (
$b$) for 
${\small \text{AR}}=1$. One can immediately observe that there are two clusters of eigenfrequencies and they in fact correspond to the P3-1 and P3-2 structures found in the nonlinear DNS around the Hopf bifurcation. By comparison, one can notice that the vortex shedding frequency of the P3-1 wake is consistent with the eigenfrequency of eigenmode B, while the vortex shedding frequency of the P3-2 wake is consistent with the eigenfrequency of eigenmode A. As we have investigated above, when some perturbation is added on the initial condition (see the red arrow in panel 
$b$), the flow will transition to the P3-2 structure. Besides, the wake P3-2 (eigenmode A) will appear at 
$Re>Re_{c2}$ as the nonlinearly saturated state. All the flow states superposed by the cyan dashed lines are transient. These results are consistent with the DNS investigations in the end of § 4.2.2 that the flow at 
${\small \text{AR}}=1$ bifurcates supercritically around the Hopf bifurcation. Moreover, in the validation § 3.1, we mentioned that the difference of vortex shedding frequency between the present result (
$St_{P3\text {-}2}=0.142$) and the DNS result in Inoue & Sakuragi (Reference Inoue and Sakuragi2008) (
$St=0.127$) at 
${\small \text{AR}}=1, Re=300$ is 10.5 %. We think that the frequency obtained by Inoue & Sakuragi (Reference Inoue and Sakuragi2008) may be the frequency of P3-1 (our 
$St_{P3\text {-}1}=0.1295$, which is closer to their value with 
$Err=1.9\,\%$) in present work. According to the results of our nonlinear DNS and LSA, P3-1 is a transitional state. These results pertain to 
${\small \text{AR}}=1$. In the case of other 
${\small \text{AR}}$ values, the flow bifurcation may be different. In a different but related context, Sheard et al. (Reference Sheard, Thompson and Hourigan2004) showed that the geometric feature size 
${\small \text{AR}}$ can change the properties of bifurcation in the flow past spheres and circular cylinders (rings). Thus, the bifurcation types of the radial flow at the other values of 
${\small \text{AR}}$ deserve to be investigated further in the future.
When 
${\small \text{AR}}\geq 1.75$, no matter whether the flow is steady or unsteady, the basic state obtained by the SFD method and the time-average method has two mutually perpendicular symmetric planes. The eigenmodes obtained based on such symmetrical basic states also retain the symmetry of their basic state, and the typical structures are shown in figure 20. The figure shows the first five eigenmodes for the case 
${\small \text{AR}}=1.75, Re=225$, ordered by the real part of the eigenvalues from large to small from panels 
$(a)$–
$(e)$. The eigenfrequencies (the imaginary parts of the eigenvalues) of the eigenmodes II and III (panels 
$b$ and 
$c$) are both zero and the corresponding eigenmodes have flow structures which appear like long strips extending streamwise. On the other hand, the other eigenmodes with non-zero eigenfrequencies have distinct wave structures in the streamwise direction.
Figure 20. The first five eigenmodes I–V based on the mean flow with two symmetric planes at 
$Re=225$ and 
${\small \text{AR}}=1.75$. The eigenvalues of (
$a$–
$e$) are 
$1.432\times 10^{-6}+ 0.1201\boldsymbol {i}$, 
$-2.662\times 10^{-3}+0\boldsymbol {i}$, 
$-9.895\times 10^{-3}+0\boldsymbol {i}$, 
$-2.794\times 10^{-2}+0.06716\boldsymbol {i}$ and 
$-3.526\times 10^{-2}+0.08008\boldsymbol {i}$, respectively. The 
$Q$-criterion isosurfaces (
$Q=0$) are coloured by the 
$x$-component of the vorticity ranging from 
$-10^{-3}$ to 
$10^{-3}$.
4.4. $Re_c$–${\small \text{AR}}$ diagram
To summarise this work, we present a 
$Re_c$–
${\small \text{AR}}$ diagram in figure 21. The preceding §§ 4.3 and 4.2 mainly discussed the various flow patterns and global modes with reference to 
${\small \text{AR}}=1$. A key focus of this work is to study the effect of 
${\small \text{AR}}$ on the flow transition in the flow past a finite-length cylinder (experiencing different flow patterns), which has not been investigated systematically in previous works for this flow. Moreover, combining present works and the previous research on the flow past a finite cylinder with 
${\small \text{AR}}>2$, the transition scenarios of this radial flow past a cylinder with two free ends can be better understood.
Figure 21. Different flow regimes in the 
$Re$–
${\small \text{AR}}$ plane. The dashed lines represent the neutral stability. The triangle ‘
$\triangle$’ represents the experimental data of the flow past a finite cylinder with hemispherical ends for the transition from a stationary state to a time-dependent flow from Schouveiler & Provansal (Reference Schouveiler and Provansal2001). The circle ‘
$\circ$’ represents the DNS of Inoue & Sakuragi (Reference Inoue and Sakuragi2008) (using a finite-difference method) results of the finite cylinder for the transition from steady flow to periodic vortex shedding. Critical Reynolds numbers in the flow past a sphere (Johnson & Patel Reference Johnson and Patel1999) and the axial flow past a cylinder with 
${\small \text{AR}}=1$ (Pierson et al. Reference Pierson, Auguste, Hammouti and Wachs2019) have also been superposed. The data of sphere (Johnson & Patel Reference Johnson and Patel1999) is placed along the vertical line of 
${\small \text{AR}}={\rm \pi} /4\approx 0.785$.
Before embarking on the discussion, we mention that in addition to 
$Re_{c1}$ (red) and 
$Re_{c2}$ (green) the definitions of which we have discussed in detail in previous sections, we also define 
$Re_{c3}$ which marks the transition from periodic flow to chaotic flow. More specifically, Sipp & Lebedev (Reference Sipp and Lebedev2007) and Turton et al. (Reference Turton, Tuckerman and Barkley2015) have proved theoretically and numerically that if a flow is dominated by a single frequency (i.e. monochromatic wave oscillation), the eigenfrequency of the linearisation operator based on the mean flow is equal to the nonlinear frequency. Based on this, 
$Re_{c3}$ is defined as follows. We increase 
$Re$ by 5 each simulation from a lower value of 
$Re$ and keep track of the dominant eigenfrequency as the Reynolds number increases. When the difference between the nonlinear vortex frequency and the dominant eigenfrequency obtained by the linear IRAM code is consecutively greater than 1 %, we define the corresponding 
$Re$ as the third critical Reynolds number 
$Re_{c3}$. For example, when 
${\small \text{AR}}=1$, as we increase 
$Re$, the first time when we have a pair of consecutive 
$Re$ whose nonlinear and linear frequencies differ by greater than 1 % is 
$Re=322$ and 
$327$ and we set 
$Re_{c3}=322$ for this parameter. The physical meaning of the bifurcation point 
$Re_{c3}$ is the transition from a single frequency dominated flow to the coexistence of multi-frequency oscillations. As for 
$Re_{c4}$, it is reminded that this value delimits the two chaos states P4-0 (with 
$\bar {C}_{ly}=\bar {C}_{lz}=0$) and P4-1 (with 
$\bar {C}_{ly}\neq 0$, occurs at 
${\small \text{AR}}<1.75$), as we have discussed in the previous section.
In figure 21 we can see that with the increase of 
${\small \text{AR}}$, the critical Reynolds number 
$Re_{c1}$ first decreases and then increases, while the critical Reynolds numbers 
$Re_{c2}$ and 
$Re_{c3}$ decrease monotonically in the range of 
${\small \text{AR}}$ from 0 to 2. The relationships between 
$Re_{c2}$–
${\small \text{AR}}$ and 
$Re_{c3}$–
${\small \text{AR}}$ approximately conform to power functions, as shown by the green and blue dashed lines in the figure with the corresponding fitting formula being 
$Re_ {c2}=158.9{\small \text{AR}}^{-1.241}+117$ and 
$Re_{c3}=136{\small \text{AR}}^{-1.572}+186.8$, respectively. Besides, it can be seen from figure 4
$(b)$ that as 
$Re$ increases, 
$\bar {C}_{ly}$ will turn to zero again, that is, the wake pattern P4-1 of 
$\bar {C}_{ly}\neq 0$ will disappear, and the black-dashed neutral curve is used to fit the critical Reynolds number of this change of wake. Then, the areas where the flows with 
$\bar {C}_{ly}\neq 0$ exist are combined into a thumb-shaped area surrounded by the black and red dashed lines in figure 21.
We have also superposed the results of the critical 
$Re$ (transitioning from a steady flow to a periodic flow, so should be compared with our 
$Re_{c2}$) in the radial flow around a short cylinder with double hemispherical free ends in Schouveiler & Provansal (Reference Schouveiler and Provansal2001) (experiments) and the double plane free ends in Inoue & Sakuragi (Reference Inoue and Sakuragi2008) (numerical simulations using the finite-difference method). Their results are supposed to be compared with our green curves and squares. It can be seen that our results are closer to the experimental results of Schouveiler & Provansal (Reference Schouveiler and Provansal2001), and are different from the DNS results of Inoue & Sakuragi (Reference Inoue and Sakuragi2008). The Hopf bifurcation result of Inoue & Sakuragi (Reference Inoue and Sakuragi2008) appears to be close to our 
$Re_{c1}$–
${\small \text{AR}}$ curve obtained. But the regular bifurcation of a radial flow past the cylinder (
${\small \text{AR}}<1.75$), generating the steady flows with one symmetric plane, found in this paper was not reported in Inoue & Sakuragi (Reference Inoue and Sakuragi2008). Nevertheless, this regular bifurcation phenomenon has been widely recognised in the flow past small 
${\small \text{AR}}$ bluff bodies such as a sphere (Johnson & Patel Reference Johnson and Patel1999; Tomboulides & Orszag Reference Tomboulides and Orszag2000; Thompson et al. Reference Thompson, Leweke and Provansal2001; Sheard et al. Reference Sheard, Thompson and Hourigan2004), axial flow around short cylinders with 
${\small \text{AR}}=1$ (Pierson et al. Reference Pierson, Auguste, Hammouti and Wachs2019) and ellipsoids (Tezuka & Suzuki Reference Tezuka and Suzuki2006; Sheard et al. Reference Sheard, Thompson and Hourigan2008). Indeed, we have also superposed the 
$Re_{c1}$ and 
$Re_{c2}$ results of the flows past a sphere (Johnson & Patel Reference Johnson and Patel1999) in figure 21. Following the discussion in the previous section, we place these data along the vertical line for 
${\small \text{AR}}={\rm \pi} /4\approx 0.785$. The 
$Re_{c1}$ value for the cylinder with 
${\small \text{AR}}=0.75$ is close to that of a sphere in the range of 
$210< Re_{c1}<212$ (reversed triangle), while the value of 
$Re_{c2}$ at 
${\small \text{AR}}=0.75$ is larger than that of the sphere's (
$270< Re_{c2}<280$), see the black triangle pointing right. Note that this is related to the fact that we place the results of Johnson & Patel (Reference Johnson and Patel1999) at 
${\small \text{AR}}\approx 0.785$. Besides, the results of Pierson et al. (Reference Pierson, Auguste, Hammouti and Wachs2019) on axial flows past an 
${\small \text{AR}}=1$ cylinder have also been superimposed in the figure. The current three critical Reynolds numbers of the radial flow past a cylinder with 
${\small \text{AR}}=1$ are all smaller than those in Pierson et al. (Reference Pierson, Auguste, Hammouti and Wachs2019)'s DNS results (
$Re_{c1}\approx 278, Re_{c2}\approx 355$ and 
$Re_{c3}\approx 395$, data found in their supplementary material) for the axial flow. (The definition of 
$Re_{c3}$ in Pierson et al. (Reference Pierson, Auguste, Hammouti and Wachs2019) is above which the attractor in the drag coefficient-lift coefficient diagram becomes chaotic. This is in principle different than our definition, but both definitions to some extent measure the onset of chaos in the flow.) This to some extent indicates that, for the cylinder with 
${\small \text{AR}}=1$, the flow separation due to the curved surface (in our radial flow) is a more efficient destabilising mechanism than the sharp-edge separation due to the front plane end in an axial flow. This may be helpful to the discussion of how to place the finite-length cylinder in a flow to suppress or stimulate flow instability. In this vein, a more general work will be needed to study the flow transition and bifurcation past a cylinder for a series of yaw angles.
5. Conclusions
In this paper we studied systematically the effect of aspect ratio (
$0.5\leq {\small \text{AR}}\leq 2$) and Reynolds number (
$Re\leq 1000$) on the wake patterns and bifurcations (of regular and Hopf types) in the flow past a finite-length short cylinder by performing nonlinear DNS and linear stability analyses (based on both time-mean flow and base flow). The detailed transition paths of the radial flow around the short cylinder (shown in figure 21) are obtained, and combined with the previous research on 
${\small \text{AR}}>2$, the present work will help us to better understand the wake transitions of the cylinder with two free ends.
With the 
$Re$ increasing from small to large, the first flow pattern is a steady wake P1 with two mutually perpendicular symmetric planes. When 
$Re>Re_{c1}$, the second flow pattern P2 emerges and it is also a steady wake but with only one symmetric plane (which is perpendicular to the cylinder axis). This P2 pattern was not discovered in the DNS study by Inoue & Sakuragi (Reference Inoue and Sakuragi2008) on the radial flow past a short cylinder. In the studies of the flows past a sphere, ellipsoid and the axial flow past a short cylinder, a similar flow structure with 
$\bar {C}_l\neq 0$ has also been found in the experimental and numerical results. According to our results, this wake pattern P2 only exists in short cylinders with aspect ratio 
${\small \text{AR}}<1.75$; when 
${\small \text{AR}}\gtrsim 1.75$, as the 
$Re$ increases, the wake will directly transition from the pattern P1 to an unsteady wake pattern P3-0 with zero 
$\bar {C}_l$.
When 
$Re>Re_{c2}$, vortex shedding happens and the flow becomes periodic and unsteady. Two wake patterns P3-1 and P3-2 have been observed. They have different frequencies and shapes and respectively represent the vortex shedding from the side surface and free ends of the cylinder. When 
${\small \text{AR}}=1$, the P3-1 structure is transient and will eventually transition to the nonlinearly stable P3-2 pattern. The transition can take a very long time if 
$Re$ is close to the critical value and there is no disturbance in the flow field. According to our numerical analyses, the bifurcation around 
$Re_{c2}$ at 
${\small \text{AR}}=1$ is supercritical as we cannot find a nonlinearly stable solution in the subcritical regime when we change the initial conditions. We have not investigated in detail the type of Hopf bifurcation in this flow at other values of 
${\small \text{AR}}$. When 
$Re$ is further increased and is larger than 
$Re_{c3}$, the flow becomes chaotic. The main message of this work is summarised in the diagram on the relation between the aspect ratio 
${\small \text{AR}}$ and the critical Reynolds numbers. The values of 
$Re_{c2}$ and 
$Re_{c3}$ decreases as 
${\small \text{AR}}$ increases. The critical 
$Re_{c1}$ first decreases and then increases with the increase of 
${\small \text{AR}}$, and finally intersects the 
${\small \text{AR}}$–
$Re_{c2}$ curve at 
${\small \text{AR}}\approx 1.75$.
Based on mean flow and base flow, global linear stability analyses have been conducted to analyse the instability and bifurcation in the flow past the short cylinder. By linearly interpolating the real part of the dominant eigenvalue, the aforementioned two critical Reynolds numbers 
$Re_{c1}$ (regular bifurcation) and 
$Re_{c2}$ (Hopf bifurcation) obtained from DNS results can also be determined in the LSA. When 
$Re< Re_{c2}$, the eigenvalues obtained based on mean flow and base flow are the same; when 
$Re>Re_{c2}$, similar to the 2-D cylinder wake flow, the nonlinear vortex frequency is consistent with the most unstable eigenfrequency based on mean flow, but the difference between the nonlinear vortex shedding frequency and the eigenfrequency based on the mean flow gradually becomes larger when 
$Re$ increases to a non-monochromatic flow state. Finally, around the Hopf bifurcation, the previously discussed wake patterns P3-1 and P3-2 can be connected to the first two global eigenmodes A and B in the global LSA based on the mean flow. In our numerical analysis, we can find that the frequencies and structures of the first two eigenmodes A, B are very close to those of nonlinear wakes P3-2 and P3-1, respectively. Therefore, based on this strong link, by analysing the linear growth rates of these eigenmodes, the evolution of the flow can be better understood and a more accurate reduced-order model of the flow can be proposed in the future.
For future directions, an interesting topic is to investigate the bifurcation type for other values of 
${\small \text{AR}}$ using a global weakly nonlinear stability analysis or similar stability analysis following Bengana 
$et ~al.$ (Reference Bengana, Loiseau, Robinet and Tuckerman2019) based on the mean flow past a finite cylinder in LSA to determine the area of hysteresis and the transition conditions of the flow.
Acknowledgements
The simulations were performed at National Supercomputing Centre, Singapore (NSCC).
Funding
We acknowledge the financial support from the Ministry of Education, Singapore (a Tier 2 grant with the WBS no. R-265-000-661- 112).
Declaration of interests
The authors report no conflict of interest.
