2. Numerical methodology

2.1. Governing equations

The flow is modelled using the three-dimensional, incompressible Navier–Stokes equations, expressed in non-dimensional form. These equations describe the evolution of the velocity field $\boldsymbol{u} = (u, v, w)$ and the pressure field $p$ in the spatial coordinates $\boldsymbol{x} = (x, y, z)$ , corresponding to the streamwise, vertical and spanwise directions, respectively:

(2.1a) \begin{align}{\qquad\qquad\quad \boldsymbol{\unicode{x2207} \cdot u} = 0\quad \text{in }\varOmega , } \end{align}

(2.1b) \begin{align} { \dfrac {\partial \boldsymbol{u}}{\partial t} +(\boldsymbol{u \cdot \unicode{x2207}})\boldsymbol{u} = -{\boldsymbol{\unicode{x2207} }} p + \frac {1}{\textit{Re}}{\boldsymbol{\unicode{x2207} }}^2 \boldsymbol{u} \quad \text{in }\varOmega . }\end{align}

Here, $\varOmega$ represents the computational domain and $\textit{Re} = U_{\infty }D/u$ is the Reynolds number, our sole control parameter. It is defined using the characteristic length $D$ of the bluff body (that is, the side length for the cube and the diameter for the sphere), the velocity of the free flow $U_{\infty }$ and the fluid kinematic viscosity $u$ . The equations are non-dimensionalised using $D$ as the length scale, $U_{\infty }$ as the velocity scale and $D/U_{\infty }$ as the time scale.

2.2. Computational domain and boundary conditions

In our simulations, the bluff bodies under consideration, namely the cube and the sphere, are centrally positioned within the computational domain $\varOmega$ and have a characteristic length $D = 1$ . The computational domain is defined by dimensions of $(\text{In}, \text{Out}, \text{Sides}) = (10D, 50D, 10D)$ , which means that the bluff body is located $10D$ downstream from the inlet and $50D$ upstream from the outlet, with side boundaries $10D$ away from the centre of the body in both vertical $(y)$ and spanwise ( $z$ ) directions. These domain dimensions are chosen to balance accuracy and computational cost, while ensuring minimal influence of domain boundaries on the wake dynamics. This set-up aligns with or improves on previous reference studies, as discussed in Appendix B.1. By using the same domain size for both the cube and the sphere, we maintain a consistent numerical environment, facilitating direct comparisons of the flow dynamics between these two canonical bluff-body shapes and also isolating the influence of body geometry on the wake dynamics and transition phenomena.

We enforce a Dirichlet boundary condition with a uniform velocity field $\boldsymbol{u} = (U_{\infty }, 0, 0) = (1, 0, 0)$ , assuming steady freestream flow at the inlet. At the outlet, we prescribe an open boundary condition with zero normal pressure gradient (Neumann condition), allowing the flow to exit the domain without artificial reflections. A no-slip condition is imposed on the body surface to model fluid–solid interactions. Free-slip boundary conditions are imposed on all side boundaries (upper, lower and lateral), allowing the use of a coarser mesh in regions far from the body, where high resolution is not required. Figure 1 provides a detailed illustration of the computational domain, including spatial configurations and boundary conditions.

Figure 1. Schematic representation of the computational domain for flow past a cube or sphere. The cube or sphere is centrally positioned within the domain. Dashed lines indicate symmetry boundaries, while dotted lines represent open boundary conditions at the inlet and outlet.

2.4. Linear stability analysis

We characterise the stability of flows around a bluff body at various $\textit{Re}$ using global linear stability analysis, which assesses the evolution of small perturbations to a base flow, and predicts the onset of instabilities and bifurcations. The base flow, computed from the Navier–Stokes equations (2.1), is either a steady state (fixed point), $\boldsymbol{Q} = (\boldsymbol{U}, P)^{\textrm{T}}$ , or a time-periodic solution (periodic orbit), $\boldsymbol{Q}(t) = (\boldsymbol{U}(t), P(t))^{\textrm{T}}$ . We employ a matrix-free, time-stepping approach based on Krylov methods – specifically the Krylov–Schur algorithm (Stewart Reference Stewart2001) – to solve large-scale eigenvalue problems. This method supports both types of base flows, enabling a comprehensive stability analysis of fully three-dimensional configurations.

We superimpose infinitesimal perturbations $\boldsymbol{q}(\boldsymbol{x}, t)=(\boldsymbol{u},p)^{\textrm{T}}$ onto $\boldsymbol{Q}$ :

(2.2) \begin{equation} \boldsymbol{Q} + \varepsilon\boldsymbol{q} = (\boldsymbol{U} + \varepsilon\boldsymbol{u}, P + \varepsilon p)^T ,\quad \text{with } \varepsilon \ll 1. \end{equation}

This applies to steady states and extends to periodic cases, $\boldsymbol{Q}(t) + \varepsilon , \boldsymbol{q}(t)$ , as is covered in §2.4.2. Substituting (2.2) into the evolution operator and linearising with respect to $\varepsilon$ yields the linearised Navier–Stokes equations (LNSEs):

(2.3a) $$\begin{align}\boldsymbol{\unicode{x2207} \cdot u} = 0\quad \text{in}\ \varOmega,\end{align}$$

(2.3b) $$\begin{align}\dfrac {\partial \boldsymbol{u}}{\partial t} + (\boldsymbol{u} \boldsymbol{\cdot }{\boldsymbol{\unicode{x2207} }}) \boldsymbol{U} + (\boldsymbol{U} \boldsymbol{\cdot }{\boldsymbol{\unicode{x2207} }}) \boldsymbol{u} = -{\boldsymbol{\unicode{x2207} }} p + \frac {1}{\textit{Re}} \boldsymbol{\unicode{x2207}} ^2 \boldsymbol{u}, \quad \text{in }\varOmega,\end{align}$$

where $\boldsymbol{u}$ is the velocity perturbation and $p$ the pressure perturbation. The LNSE can be further compactly represented in matrix form as

(2.4) \begin{equation} \mathscr{B} \dfrac {\partial \boldsymbol{q}}{\partial t} = \mathscr{L} \boldsymbol{q}, \end{equation}

where $\mathscr{B}=((\textbf{I}, 0), (\textbf{0},0))^{\textrm{T}}$ is the (singular) mass matrix due to incompressibility and $\mathscr{L}$ is the Jacobian matrix. Assuming a normal mode decomposition via a Laplace transform in time, the infinitesimal perturbation has the form

(2.5) \begin{equation} \boldsymbol{q}(\boldsymbol{x}, t) = \hat {\boldsymbol{q}}(\boldsymbol{x}) e^{\lambda t} + \text{c.c.}, \end{equation}

where $\hat {\boldsymbol{q}}(\boldsymbol{x})$ is the spatial eigenfunction of the perturbation, dependent on all three spatial coordinates $\boldsymbol{x} = (x, y, z)$ , and $\lambda$ is the complex eigenvalue of the considered normal mode. Substituting (2.5) into (2.4) and neglecting higher-order terms yields the generalised eigenvalue problem:

(2.6) \begin{equation} \lambda \mathscr{B}\hat {\boldsymbol{q}} = \mathscr{L} \hat {\boldsymbol{q}}, \end{equation}

where the complex eigenvalue $\lambda = \sigma + i \omega$ governs the temporal behaviour of the perturbation. The real part, $\sigma$ , represents the growth rate, dictating whether perturbations amplify ( $\sigma \gt 0$ ) or decay ( $\sigma \lt 0$ ), while the imaginary part, $\omega = 2\pi f$ , corresponds to the frequency $f$ . A base flow $\boldsymbol{U}$ is deemed globally unstable if the eigenvalue problem (2.6) admits at least one solution with a positive growth rate. The corresponding eigenmode, $\hat {\boldsymbol{q}}(\boldsymbol{x})$ , termed a global eigenmode, represents the spatial structure of the instability throughout the domain. The emergence of global eigenmodes requires that at least one area of the flow exhibits a synchronised, self-sustained oscillation described by a global eigenfunction (Chomaz, Huerre & Redekopp Reference Chomaz, Huerre and Redekopp1988). In the opposite case, where all eigenvalues have negative growth rates, the base flow is globally stable (not temporally self-sustained).

The next two subsections address solving eigenvalue problem (2.6) for steady-state (fixed point) and time-periodic (periodic orbit) base flows.

2.4.1. Fixed points

To solve the problem (2.6), building and storing the Jacobian matrix $\mathscr{L}$ is computationally impractical due to the large number of degrees of freedom in three-dimensional simulations. Instead, we employ a time-stepping approach, based on Krylov methods, that leverages the linearity of the problem and computes the action of the Jacobian over time. We approximate the action of the Jacobian $\mathscr{L}$ via the matrix exponential of the operator $\textbf{M}_{\tau } = e^{\mathscr{L} \tau }$ over a finite sampling time $\tau$ . This leads to the large-scale eigenvalue problem:

(2.7) \begin{equation} \textbf{M}_{\tau } \hat {\boldsymbol{q}} = { \mu} \hat {\boldsymbol{q}}, \end{equation}

where ${ \mu}$ are the eigenvalues of $\textbf{M}_{\tau }$ , and $\hat {\boldsymbol{q}}$ are the corresponding eigenmodes. The eigenvalues of $\mathscr{L}$ are then recovered with the natural logarithm transformation $\lambda = \ln ({ \mu} ) / \tau$ .

Figure 2 illustrates common eigenspectrum structures and how the leading eigenvalue governs the bifurcation type. If the leading eigenvalue crosses the imaginary axis $(\sigma \gt 0)$ with zero imaginary part ( $\omega = 0$ ) and symmetry breaking occurs, the system undergoes a pitchfork bifurcation (figure 2 a), resulting in a steady-state solution. If the leading eigenvalue crosses into the unstable half-plane with $\omega \neq 0$ , an Hopf bifurcation occurs (figure 2 b), leading to periodic dynamics (periodic orbits). To efficiently solve (2.7), we use the NekStab toolbox, which integrates the highly parallel time-stepping capabilities of Nek5000 with Krylov methods. The effectiveness of this methodology in fixed point cases has been demonstrated in several studies (Loiseau et al. Reference Loiseau, Robinet, Cherubini and Leriche2014; Bucci et al. Reference Bucci, Puckert, Andriano, Loiseau, Cherubini, Robinet and Rist2018; Picella et al. Reference Picella, Loiseau, Lusseyran, Robinet, Cherubini and Pastur2018; Loiseau et al. Reference Loiseau, Bucci, Cherubini and Robinet2019; Bucci et al. Reference Bucci, Cherubini, Loiseau and Robinet2021).

Figure 2. Schematic of bifurcations experienced by a steady-state solution (i.e. fixed point). This diagram illustrates the two types of bifurcations that occur in steady-state flows as a parameter (e.g. Reynolds number) is varied. The shaded region indicates the stable part of the spectrum, i.e. the lower complex half-plane $\textit{Re}\{\lambda \} = \sigma \lt 0$ . Adapted from Frantz et al. (Reference Frantz, Loiseau and Robinet2023).

2.4.2. Periodic orbits

Now, we assess the stability of a periodic orbit $\boldsymbol{Q}(t) = (\boldsymbol{U}(t), P(t))^{\textrm{T}}$ by performing a Floquet analysis. This method, designed for time-periodic systems, evaluates the net effect of small perturbations over one full period $\tau$ , determining whether they grow or decay. This evolution is captured by the monodromy matrix $\textbf{M}_\tau$ , which encapsulates the cumulative effect of the linearised dynamics over the period $\tau$ . In practice, it is obtained by integrating the linearised system:

(2.8) \begin{equation} \textbf{M}_\tau = \int _0^{\tau } \mathscr{L}(t) \,\textrm {d}t, \end{equation}

where $\mathscr{L}$ is the linearised Jacobian operator around the periodic base flow $\boldsymbol{Q}(t)$ from which we can assess how perturbations behave over successive periods. The eigenvalue problem to solve therefore writes

(2.9) \begin{equation} \textbf{M}_\tau \hat {\boldsymbol{q}} = { \mu} \hat {\boldsymbol{q}}, \end{equation}

where $\hat {\boldsymbol{q}}$ denote the Floquet modes (i.e. the eigenvectors of the monodromy matrix $\textbf{M}_\tau$ ) and ${ \mu}$ are the corresponding Floquet multipliers (i.e. eigenvalues). The magnitude of the Floquet multipliers indicates the stability of the periodic orbit: if $|{ \mu} | \lt 1$ for all ${ \mu}$ in the spectrum of $\textbf{M}_\tau$ , perturbations decay over one period and the periodic orbit is stable to small disturbances; if there exists at least one eigenvalue ${ \mu}$ such that $|{ \mu} | \gt 1$ , perturbations grow over one period, indicating that the periodic orbit is unstable. The Floquet multiplier ${ \mu}$ is related to the Floquet exponent $\lambda$ by ${ \mu} = e^{\lambda \tau }$ , where $\lambda = \sigma + i \omega$ is a complex number. Here, $\sigma$ represents the average growth rate per orbital period (real part) and $\omega$ corresponds to an additional frequency component introduced by the perturbation (imaginary part).

Figure 3. Schematic of bifurcations experienced by a time-periodic base flow (or periodic orbit). This figure depicts three typical bifurcation types in time-periodic systems. The shaded region indicates the stable part of the spectrum, i.e. the unit disk $|{ \mu} | = 1$ . Adapted from Frantz et al. (Reference Frantz, Loiseau and Robinet2023).

Figure 3 illustrates common scenarios of how the leading Floquet multiplier governs the type of bifurcation that the system may undergo.

1. (i) Pitchfork bifurcation (figure 3 a) ( $|{ \mu} | \gt 1$ and ${ \mu} \in \mathbb{R}^+$ ): new periodic solutions with different amplitudes or phases emerge from the original orbit. This type of bifurcation often occurs in systems with certain symmetries and can lead to multiple stable periodic solutions coexisting.

2. (ii) Period-doubling bifurcation (figure 3 b) ( $|{ \mu} | \gt 1$ and ${ \mu} \in \mathbb{R}^-$ ): a period-doubling or flip bifurcation occurs. Here, the Floquet multiplier moves through ${ \mu} = -1$ , indicating that perturbations change sign every period. A new periodic solution with twice the period of the original orbit emerges.

3. (iii) Neimark–Sacker bifurcation (figure 3 c) ( $|{ \mu} | \gt 1$ and ${ \mu} \in \mathbb{C}$ with a non-zero imaginary part): a pair of complex conjugate Floquet multipliers move across the unit circle, introducing a new frequency into the system. This bifurcation is thus also known as a secondary Hopf bifurcation. The resulting motion is quasiperiodic, characterised by oscillations with two incommensurate frequencies.

As in the fixed point methodology, we employ the matrix-free approach to circumvent the computational impracticability of constructing the monodromy matrix $\textbf{M}_\tau$ . Instead of forming $\textbf{M}_\tau$ explicitly, we compute its action on perturbation vectors by integrating the linearised system over one period $\tau$ . Again, we employ the NekStab toolbox, which uses Krylov methods to compute dominant Floquet multipliers and corresponding modes, enabling the stability assessment of periodic orbits in large-scale three-dimensional flows; for a comprehensive overview of Floquet analysis and its applications, we refer the reader to the review by Frantz et al. (Reference Frantz, Loiseau and Robinet2023).

2.4.3. Numerical parameters

For the stability analysis of fixed points, the time-stepping approach is based on an integration sampling period equal to $\tau = 1.25$ . This value is determined by the product of the number of iterations ( $n$ ) required by the linear solver (i.e. (2.3)), and the time step ( ${\textrm d}t$ ) used in both cube and sphere simulations ( $\tau =n\times {\textrm d}t$ ). This choice follows the recommendations of Frantz et al. (Reference Frantz, Loiseau and Robinet2023), which enables the efficient execution of the eigenproblem within a single supercomputer run. For the Floquet analysis involving periodic orbits, $\tau$ is set to precisely match the period of the orbit, whose evaluation is described in Appendix D.2. Despite the large domain size, sponge regions were used at the inlet (one element) and outlet (five elements) to minimise potential boundary reflections, acting to gradually zero perturbations, in line with best practices for stability analysis. For the eigenvalue computations, a Krylov basis of size $m = 100$ was employed, with convergence achieved when at least 20 eigenvalues reached an absolute residual of $10^{-6}$ , while absolute tolerances of $10^{-10}$ were employed for the velocity and pressure solvers.

3. Primary bifurcation

Figure 4. Comparative visualisation of base flows at the onset of a primary pitchfork bifurcation. (a,c,e) Orthogonallysymmetric steady base flow around a cube at a supercritical $\textit{Re}=207\gt \textit{Re}_{c,1}$ . (b,d,f) Axisymmetric steady base flow around a sphere at a supercritical $\textit{Re}=214\gt \textit{Re}_{c,1}$ . For each geometry, the figure is extracted at perpendicular cross-sections $z=0$ , $y=0$ and $x=0$ . Streamwise velocity maps with $u\in [-0.1,1.1]$ are overlaid by white isolines for $u=0$ , which delineate the recirculation regions indicative of reverse flow, and black isolines for $u=1$ , highlighting the freestream flow.

Figure 5. Primary bifurcation for (a,c) the cube and (b,d) the sphere. (a,b) Eigenspectra in the $\sigma$ – $St^L$ plane, based on the steady symmetric base flow, evidencing a pitchfork bifurcation. (c,d) Trajectory of the leading eigenvalue’s amplification rate $\sigma$ in relation to the Reynolds number $\textit{Re}$ . The vertical dashed line defines the critical Reynolds number $\textit{Re}^L_{c,1}=206.39$ for the cube case, interpolated from the fitting linear equation $\sigma (\textit{Re}) = 0.0019 \boldsymbol{\cdot }(\textit{Re} - \textit{Re}^L_{c,1})$ in panel (c) and $\textit{Re}^L_{c,1}=213.11$ for the sphere case, interpolated from the fitting linear equation $\sigma (\textit{Re}) = 0.0028 \boldsymbol{\cdot }(\textit{Re} - \textit{Re}^L_{c,1})$ in panel (d).

Cube. From nonlinear simulations, Meng et al. (Reference Meng, An, Cheng and Kimiaei2021) identified a primary bifurcation in the flow past a static cube at a critical Reynolds number $\textit{Re}_{c,1}\approx 207$ . Our investigation, which started with a DNS at $\textit{Re}=200$ , confirmed the steady state of the flow and revealed, as reported by Meng et al. (Reference Meng, An, Cheng and Kimiaei2021), an orthogonally symmetric recirculation zone directly behind the cube.

Using the Newton–Krylov solver described in Appendix D.1, we meticulously traced the evolution of the orthogonally symmetric flow solution from $\textit{Re}=200$ to $\textit{Re}=212$ . This process unveiled a transition to an unstable base flow at $\textit{Re}=207$ , perfectly orthogonally symmetric, as depicted in perpendicular cross-sections ( $z=0, y=0, x=0$ ) in figure 4(a,c,e). It is crucial to note that our analysis did not impose additional symmetry constraints (such as half- or quarter-domain simulations), ensuring that the flow symmetry emerges naturally from the governing equations.

Using a linear solver (discretising the LNSE around the present base flow solution), the large-scale eigenvalue problem (2.7) is then solved, as described in §2.4.1. The results are reported in figures 5(a) and 5(c). Figure 5(a) represents the eigenspectra in the $\sigma$ – $St^L$ (where the superscript $L$ denotes results obtained from linear analysis) plane corresponding to three distinct Reynolds numbers and based on the orthogonal-symmetric steady base flow shown in figure 4 ( $St^L$ and $\omega$ are linked by the relation $St^L = \omega / (2 \pi )$ ). These eigenspectra indicate the destabilisation of the base flow as a purely real eigenvalue ( $St^L=0$ ) migrates into the positive $\sigma$ half-plane, uncovering a pitchfork bifurcation.

Figure 5(c) illustrates the relationship between the amplification rate $\sigma$ and the Reynolds number. We systematically sampled Reynolds numbers from $\textit{Re}=200$ to $\textit{Re}=212$ with increments $\varDelta _{\textit{Re}}=\{3,2,1,1,2,3\}$ , ensuring a monotonically increasing sequence. This sampling strategy reveals a clear linear trend in the amplification rates as the Reynolds number increases. The monotonic nature of this relationship allows for reliable linear interpolation, enabling us to estimate the primary bifurcation threshold with high accuracy. Using this method, we determine that the critical Reynolds number is $\textit{Re}^L_{c,1}=206.39$ . This value is closely aligned with those reported in the literature (cf. §1.2.4 and tables 3, 4).

Sphere. Adopting a similar analytical trajectory for the sphere, over Reynolds numbers $\textit{Re}=207$ to $\textit{Re}=219$ with an identical $\varDelta _{\textit{Re}}$ step sequence, we identified an axisymmetric yet unstable base flow at $\textit{Re}=214$ , shown in figure 4(b,d,f). From this base flow, the eigenspectra analysis reported in figures 5(b) and 5(d) determines the primary (pitchfork) bifurcation threshold for the sphere at $\textit{Re}^L_{c,1}=213.11$ . This threshold is in agreement with the literature-reported value of $\textit{Re}_{c,1}\approx 212$ (cf. §1.2.1 and tables 3, 4).

Discussion and comparison. By comparing cube and sphere base flows (figure 4), it can be observed that for an identical distance to $\textit{Re}^L_{c,1}$ , the recirculation length behind the cube is longer than that of the sphere. This may be explained by the presence of the secondary recirculation bubble formed along the lateral and top/bottom faces of the cube (see the $u=0$ isolines in figure 4 a,c,e), which participates in the elongation of the main recirculation zone behind the cube (Klotz et al. Reference Klotz, Goujon-Durand, Rokicki and Wesfreid2014). In their study, Meng et al. (Reference Meng, An, Cheng and Kimiaei2021) also evidenced the presence of four additional pairs of vortices generated on the lateral faces of the cube, directly ‘feeding’ the main vortex core. This phenomenon does not exist in the sphere case due to the curvature of the body, which prevents the formation of any secondary recirculation bubble.

Figure 6. Primary bifurcation. Three-dimensional view of isocontours of (a) streamwise $(u')$ , (b) vertical $(v')$ and (c) spanwise $(w')$ perturbation velocities plotted for the levels $\pm 0.01$ of the real part of the leading unstable mode at $\textit{Re} = 207$ for the cube case and $\textit{Re} = 214$ for the sphere case.

Concerning eigenanalysis, we observed that at the respective critical values $\textit{Re}^L_{c,1}$ , the emergence of a positive real eigenvalue signals a pitchfork bifurcation for both geometries. This bifurcation marks the destabilisation of spatial symmetry, heralding a transition from a maximally symmetric wake to a planar symmetric state as the Reynolds number surpasses $\textit{Re}^L_{c,1}$ . Furthermore, figure 6 shows the unstable mode of the primary (pitchfork) instability, capturing its development in proximity to the steady symmetric fixed point (i.e. at $\textit{Re} \gtrapprox \textit{Re}_{c,1}$ ) for both cube and sphere configurations. The isosurfaces of the components of the perturbation velocity field show identical spatial symmetries in the $(x, y)$ and $(x, z)$ planes between the cube and the sphere case. These results therefore lead us to the conclusion that the corner-induced vortices which induce secondary flows around the cube surface do not impact the global dynamics, but act as an extra source making the main recirculation bubble larger and therefore destabilising the flow faster. This is captured by the lower critical Reynolds number in the cube case.

4. Secondary bifurcation

Figure 7. Comparative visualisation of base flows at the onset of a secondary Hopf bifurcation. (a,c,e) Planar-symmetric steady base flow around a cube at a supercritical $\textit{Re}=252\gt \textit{Re}_{c,2}$ . (b,d,f) Planar-symmetric steady base flow around a sphere at a supercritical $\textit{Re}=275\gt \textit{Re}_{c,2}$ . For each geometry, the figure is extracted at perpendicular cross-sections $z=0$ , $y=0$ and $x=0$ . Streamwise velocity maps with $u\in [-0.1,1.1]$ are overlaid by white isolines for $u=0$ and black isolines for $u=1$ .

Figure 8. Secondary bifurcation for (a,c) the cube and (b,d) the sphere. (a,b) Eigenspectra in the $\sigma$ – $St^L$ plane, based on the bifurcated steady planar-symmetric base flow, evidencing a Hopf bifurcation. (c,d) Evolution of the amplification rate with respect to the Reynolds number. The vertical dashed line defines the critical Reynolds number $\textit{Re}^L_{c,2}=250.61$ for the cube case, interpolated from the fitting linear equation $\sigma (\textit{Re}) = 0.0039 \times (\textit{Re} -\textit{Re}^L_{c,2})$ in panel (c) and $\textit{Re}^L_{c,2}=273.67$ for the sphere case, interpolated from the fitting linear equation $\sigma (\textit{Re}) = 0.0037 \times (\textit{Re}-\textit{Re}^L_{c,2})$ in panel (d).

Figure 9. Secondary bifurcation. Three-dimensional view of isocontours of (a) streamwise $(u')$ , (b) vertical $(v')$ and (c) spanwise $(w')$ perturbation velocities plotted for the levels $\pm 0.01$ of the real part of the periodic leading unstable mode associated to the shedding frequency $St_1$ , at $\textit{Re} = 252$ for the cube case and $\textit{Re} = 275$ for the sphere case.

Our investigation of secondary bifurcations looks into the post-primary instability dynamics in both cube and sphere configurations.

Cube. At $\textit{Re}=245$ , surpassing the primary instability threshold, a direct simulation in the cube scenario revealed a transition to a non-orthogonal asymmetric steady wake. We then performed a rigorous stability analysis, employing the Newton–Krylov solver to compute new fixed points, uncovering a planar-symmetric base flow which is temporally unstable at $\textit{Re}=252$ , delineated in figure 7(a,c,e). The corresponding eigenspectra, depicted in figure 8(a), show a pair of complex conjugate eigenvalues (i.e. with $St^L \neq 0$ ) that become unstable as the $\textit{Re}$ increases from 245 to 257, exhibiting the emergence of an Hopf bifurcation and a fundamental frequency $St_1^L$ , responsible for periodic vortex shedding in the wake region. Then, figure 8(c) gives the linear relation between $\sigma$ and $\textit{Re}$ , for Reynolds numbers ranging from $\textit{Re}=245$ to $\textit{Re}=257$ with steps $\varDelta _{\textit{Re}}=\{3,2,1,1,2,3\}$ , allowing us to estimate the secondary bifurcation threshold at $\textit{Re}^L_{c,2}=250.61$ . This critical value aligns with the ranges reported in the literature, notably by Khan et al. (Reference Khan, Sharma and Agrawal2020b ) ( $250 \lt \textit{Re}_{c,2} \lt 276$ ) and Meng et al. (Reference Meng, An, Cheng and Kimiaei2021) ( $\textit{Re}_{c,2}\approx 252$ ). Using this interpolated critical value, we can interpolate again and obtain the associated critical frequency $St_{c,1}^L=0.0927$ , highlighted by the dash-dotted vertical line in figure 8(a). We note that the estimated critical frequency is also in very good agreement with the frequency $St_1 = 0.0975$ found by Meng et al. (Reference Meng, An, Cheng and Kimiaei2021) in their nonlinear simulations at $\textit{Re}=270$ .

Sphere. Similarly, the Newton–Krylov solver computes an unstable planar-symmetric base flow for the sphere at $\textit{Re}=275$ (figure 7 b,d,f) and the stability analysis, showcased in figures 8(b) and 8(d) and spanning from $\textit{Re}=268$ to $\textit{Re}=280$ , led to the interpolation of the secondary bifurcation (Hopf) threshold at $\textit{Re}^L_{c,2}=273.67$ and to an estimated critical fundamental frequency equal to $St_{c,1}^L=0.1282$ . These findings are in close agreement with the thresholds found in the literature ( $\textit{Re}_{c,2} = 275 \pm 6$ ) (cf. §1.2.1 and tables 3, 4) and specifically with the stability analysis conducted by Citro et al. (Reference Citro, Siconolfi, Fabre, Giannetti and Luchini2017) ( $\textit{Re}_{c,2} = 271.8$ and $St_{c,1} = 0.129$ ), emphasising the precision of the present methodology.

Discussion and comparison. This stability analysis first highlights striking similarities between the cube and sphere base flows, as well as in the unstable global eigenmode associated with the secondary bifurcation (Hopf), as shown in figure 9. In terms of discrepancies between the two bluff bodies, one can notice from the base flow computation (cf. figure 7), and for the same reason as that explained at the end of §3, the recirculation length $L_r$ for the cube is longer than the sphere’s counterpart $L_r^{\textit{cube}} \gt L_r^{\textit{sphere}}$ . The stability analysis then emphasises that the critical shedding frequency in the cube case is lower than that of the sphere, $St^{\textit{cube}} \lt St^{\textit{sphere}}$ (figure 8) and that the opposite relation stands regarding the streamwise wave length $\varLambda$ of the unstable mode, i.e. $\varLambda _{\textit{cube}} \gt \varLambda _{\textit{sphere}}$ (qualitatively observable in figure 9). The existence of a proportional dispersion relation between $\varLambda$ and $St$ in the range of Reynolds numbers under consideration in this study implies that , and having the following relation between the recirculation length and wave length therefore makes the quantitative discrepancies between the cube and the sphere cases in terms of $L_r$ , $St$ and $\varLambda$ values naturally correlated.

Figure 10. Secondary bifurcation at $\textit{Re}=275$ (sphere case only). $\varOmega _r$ isocontours of: (a) ‘BF’, planar-symmetric steady base flow; (b) ‘M’, secondary unstable global eigenmode; (c) ‘BF + M’, linear combination of the base flow and global eigenmode at an arbitrarily fixed amplitude; (d) ‘DNS’, nonlinear solution. The top-left box shows a side view. [Movie 1] (cube), [Movie 2] (sphere).

We now turn to examining the spatial mechanisms that drive the dynamical evolution of the planar-symmetric base flow following the Hopf bifurcation. Figure 10 presents a synchronised comparison of $\varOmega _r$ -criterion isocontours for the sphere case, illustrating: (a) the planar-symmetric base flow (BF); (b) the unstable global eigenmode (M); (c) the linear superposition of base flow and eigenmode at a fixed amplitude (20 % of the maximum streamwise velocity of BF) that clearly reveals the eigenmode’s effect (BF + M), and (d) the direct numerical simulation at $\textit{Re}=275 \gt \textit{Re}^L_{c,2}$ (DNS). The $\varOmega _r$ -criterion employed in figure 10 follows the Omega-Liutex/Rortex methodology (Dong, Gao & Liu Reference Dong, Gao and Liu2019), which identifies vortices by quantifying the ratio of local rigid-rotation vorticity component to total vorticity (rotation + deformation). A distinctive advantage of this criterion is its robust threshold value of $\varOmega _r = 0.52$ , which serves as a case-independent demarcation for vortex boundaries. When $\varOmega _r \gt 0.52$ , the rotational component dominates the deformation component, effectively capturing both weak and strong vortical structures within the same flow domain – an improvement over previous criteria that relied on adhoc, flow-dependent thresholds.

We note that the fundamental mechanisms at play in the cube case are equivalent, so for clarity and conciseness, the following explanation focuses on the sphere’s spatial flow structures. The ‘DNS’ part of figure 10 clearly shows the well-known periodic vortex-shedding process, featuring alternating hairpin-like structures. In our view, this phenomenon aligns with the minimal flow-elements model proposed by Cohen, Karp & Mehta (Reference Cohen, Karp and Mehta2014), which provides a ‘recipe’ for generating packets of hairpin vortices in shear flows. Their model – consisting of a shear flow, a counter-rotating vortex pair (CVP) and a streamwise oscillation – explains the mechanics of the instability leading to hairpin vortices and, more precisely, how the CVP modifies the base flow to create an inflection point, thereby amplifying the wavy vortex sheet. This amplification, in turn, destabilises the CVP, leading to its segmentation and the formation of hairpin vortex packets. Analysing more precisely the vortex shedding process in the present case, and as already observed in many numerical and experimental studies (Sakamoto & Haniu Reference Sakamoto and Haniu1990; Johnson & Patel Reference Johnson and Patel1999; Kim & Choi Reference Kim and Choi2002), the vortices shed from the sphere are alternating in a non-symmetric fashion in the $(x,y)$ plane, with the creation of hairpin rings occurring in an apparently fixed orientation, namely, in the upper branch of the wake. The present figure 10 provides more detail on this phenomenon by spatially representing the eigenmode’s action on the base flow. In fact, in the part ‘(BF + M)’ of the figure, a delayed spatial amplification of the structures formed in the lower branch of the wake can be observed compared with their counterparts formed in the upper branch. This delay in amplification is directly related to the asymmetry of the base flow in the $(x,y)$ plane, with a recirculation bubble shifted in the upper half $(x,y)$ plane, implying different local velocity gradients in the two half-planes. Due to nonlinear effects, the amplification delay observed in the spatial structures of the eigenmode (M) and its linear superimposition with the base flow (BF + M) manifests in the ‘DNS’ solution with the formation of non-alternating hairpin rings (horseshoe-shaped). These rings are aligned with the orientation of the asymmetric base flow recirculation zone.

Base flow asymmetry and hairpin formation. In other words, the formation of hairpin vortices in the nonlinear flow can be traced back to the interplay between the asymmetric base flow and the unstable eigenmode. Although the linear mode alone does not manifest as a clear hairpin, it provides the initial ‘seed’ of vorticity perturbation. Because the base flow is vertically asymmetric (due to the pitchfork bifurcation), amplification differs between the upper and lower parts of the wake (i.e. local gradients), causing a noticeable lag or phase shift in vortex roll-up. By visualising the superposition of the base flow and the linear mode, one sees that the spatial evolution of the instability directly depends on local velocity gradients; the upper region of the wake experiences earlier, stronger amplification, while the lower region lags behind. In fully nonlinear simulations, this asymmetry in amplification eventually leads to hairpin structures forming predominantly on one side first, corroborating the linear predictions that spatial gradients in the base flow drive an uneven – yet ultimately coherent – shedding pattern.

5. Tertiary bifurcation

As previously shown, a limit cycle associated with periodic dynamics is created in the vicinity of the steady planar-symmetric fixed point when the secondary bifurcation (first Hopf) occurs. The investigation of tertiary bifurcation through linear stability analysis thus implies the computation and the stabilisation of the periodic base flow (a periodic orbit). In this section, we present results of the stability characterisation of stabilised periodic orbits based on Floquet theory. We present what we believe to be the first computation of three-dimensional Floquet modes, which reveal instabilities within the limit cycle and indicate the onset of tertiary bifurcations. To the best of our knowledge, the stability analysis of tertiary bifurcations in the flow past a cube and a sphere has not been previously reported. Therefore, in this section, we present the results not in terms of the ‘cube’ and ‘sphere’ cases sequentially, but by following the two main steps of the Floquet analysis: (i) computation of the periodic base flows for the cube and sphere; and (ii) computation of the Floquet modes for both obstacles.

Computation and analysis of the UPO. Figure 11 illustrates, for both cube and sphere cases, the residual reduction during the Newton–Krylov algorithm used for the computation of periodic orbits (cf. Appendix D.2). Figure 12 delineates the variation of the limit cycle frequency $St$ and related period $\tau$ against the $\textit{Re}$ , using the Newton–Krylov solver to stabilise periodic orbits with naturally selected periods for both cases. We used initial guesses for the orbit period based on the frequency peak from the lift coefficient timeseries from preliminary simulations. The data highlight a discernible trend, fitted with a second-order polynomial, along with the critical $\textit{Re}$ interpolated from linear stability analysis results (presented hereafter). In particular, the critical frequency of the periodic base flow is interpolated and depicted for comparative analysis. Preliminary DNS at subcritical $\textit{Re}$ provided a robust initial estimate for the orbit period. This allowed a gradual increment of the $\textit{Re}$ to supercritical values, ensuring proximity to the desired solution branch and thus ensuring convergence of the Newton method.

Figure 11. Residual deflation as a function of Newton–Krylov iterations ( $ k$ ) for the stabilisation of unstable periodic orbits (UPOs) in the cube and sphere cases. The initial solutions correspond to solutions at $ \textit{Re} = 280$ for the cube and $ \textit{Re} = 330$ for the sphere. Due to the proximity of these initial solutions to the target solutions, the initial residuals are of the order of $ 10^{-4}$ . The dotted line represents the specified tolerance level, while the curves demonstrate the convergence of the residual norm ( $ \|r\|$ ).

Figure 12. Strouhal numbers ( $St$ ) and periods ( $\tau$ ) of the limit cycle computed using the Newton–GMRES solver as functions of the Reynolds number ( $\textit{Re}$ ) in (a) the cube case and (b) the sphere case. The black dashed lines represent second-order polynomial fits to the data points.

A visualisation of the stabilised unstable limit cycle around the cube at supercritical Reynolds number $\textit{Re}=281 \gt \textit{Re}_{c,3} = 280.31$ and around the sphere at supercritical Reynolds number $\textit{Re}=332 \gt \textit{Re}_{c,3}=330.75$ , is shown in figure 13. The stabilised orbit period $\tau$ is here denoted as $\tau _1$ . The $\varOmega _r = 0.52$ isocontours (figure 13 c) are plotted at every quarter of the shedding period $\tau _1$ , showing the periodic shedding of the hairpin vortical structures and their transport in the wake. Notably, there is a perfect match between the snapshots of the UPO at $t = 0$ and $t =\tau _1$ . The time-periodic base flow of flow past the cube and the sphere are also represented by the $(u(t),v(t),\partial _t v(t))$ phase-portraits (figure 13 a,b), depicting for both obstacles a clear limit cycle of period $\tau _1$ . Figure 14 represents the stabilised UPO at $t=0$ in different perpendicular cross-sections, highlighting the planar-symmetry of the periodic base flow for both the cube and the sphere.

Figure 13. Comparative visualisation of the stabilised unstable limit cycle (i.e. periodic base flow) at the onset of a tertiary Neimark–Sacker bifurcation, at $\textit{Re}=281 \gt \textit{Re}_{c,3}$ for the cube case and $\textit{Re}=332\gt \textit{Re}_{c,3}$ for the sphere case. (a and b) $(u,v,\partial v / \partial t)$ phase-portrait of periodic base flow in cube and sphere case, respectively. (c) Perspective view of $\varOmega _r = 0.52$ isocontours of the periodic base flow represented for every quarter period of a shedding cycle. Green and purple arrows depict, in the cube and sphere case, respectively, the spatial location of a given hairpin through a complete shedding period $\tau _1$ .

Figure 14. (Complement to figure 13). Comparative visualisation of a single snapshot on the stabilised unstable limit cycle in perpendicular cross-sections at $t = 0$ , at the onset of a tertiary Neimark–Sacker bifurcation. The figure highlights the planar-symmetry of the cube and sphere’s stabilised periodic base flow.

Figure 15. Tertiary bifurcation for (a,c,e) the cube and (b,d,f) the sphere. (a,b) Eigenspectra of the Floquet multipliers developing in the vicinity of the stabilised limit cycles at various Reynolds numbers, evidencing a Neimark–Sacker bifurcation. (c,d) Eigenspectra in the $\sigma$ – $St^L$ plane, based on the stabilised limit cycles at various $\textit{Re}$ , evidencing the new frequency $St_{c,2}^L$ . Note that zero eigenvalues corresponding to the base flow derivative are omitted from the plot. (e,f) Evolution of the amplification rate with respect to the $\textit{Re}$ . The vertical dashed line defines the critical Reynolds number $\textit{Re}^L_{c,3}=280.31$ for the cube case, interpolated from the fitting linear equation $\sigma (\textit{Re}) = 0.0020\boldsymbol{\cdot }(\textit{Re}-\textit{Re}^L_{c,3})$ in panel (e) and $\textit{Re}^L_{c,3}=330.75$ for the sphere case, interpolated from the fitting linear equation $\sigma (\textit{Re}) = 0.0011\boldsymbol{\cdot }(\textit{Re}-\textit{Re}^L_{c,3})$ in panel (f).

Computation and analysis of the Floquet modes. Figure 15 represents the results of the Floquet stability analysis (cf. §2.4.2) for the cube (panels a, c, e) and the sphere case (panels b,d,f). Figures 15(a) and 15(b) display for each obstacle the eigenspectrum of the Floquet multipliers ${ \mu} _i$ (corresponding to the eigenvalues of the monodromy matrix ${\textbf{M}}_{\tau }$ ) relative to the unit disk $\vert { \mu} \vert =1$ , where the inner grey area represents the stable portion of the spectrum. For both cube and sphere, this eigenspectrum reveals the evolution of a complex conjugate pair of multipliers, progressively leaving the unit circle out of the $\mathcal{R}\textit{e} ({ \mu} )$ axis as $\textit{Re}$ grows, indicating a Neimark–Sacker bifurcation (also called a secondary Hopf bifurcation). As a result, the mode associated to a new frequency $St_2$ (incommensurate with respect to the shedding frequency $St_1$ that arises after the onset of the first Hopf bifurcation) becomes unstable, leading to quasiperiodic dynamics (the quasiperiodic regime will be temporally and spatially analysed in §6). Figures 15(c) and 15(d) give the mapping of the spectrum of the monodromy matrix (figure 15 a,b) towards the corresponding spectrum in the $\sigma {-}St^L$ plane of the Jacobian matrix, following the relation $\lambda _i = \ln ({ \mu} _i)/ \tau$ . Eventually, figures 15(e) and 15(f) portray the linear progression of the growth rates $\sigma$ with respect to $\textit{Re}$ .

Concerning the cube case, the tertiary bifurcation is identified at a critical Reynolds number of $\textit{Re}^L_{c,3}=280.31$ (figure 15 e), aligning with DNS predictions of $\textit{Re}_{c,3}=282$ given by Meng et al. (Reference Meng, An, Cheng and Kimiaei2021). The new frequency of the unstable leading modes is interpolated from the critical value of $\textit{Re}^L_{c,3}$ and is estimated to $St_{c,2}^L=0.0276$ (see the vertical dash-dotted line in figure 15 c), which matches the $St_2 = 0.0281$ value reported by Meng et al. (Reference Meng, An, Cheng and Kimiaei2021) in their nonlinear simulations at $\textit{Re}=282$ .

Regarding the sphere case, the Neimark–Sacker bifurcation emerges at a critical Reynolds number of $\textit{Re}^L_{c,3}=330.75$ (figure 15 f), which is slightly below the broad range of $350 \lt \textit{Re}_{c,3} \lt 480$ reported in the DNS results of Kalro & Tezduyar (Reference Kalro and Tezduyar1998), Mittal (Reference Mittal1999), Mittal & Najjar (Reference Mittal and Najjar1999), Lee (Reference Lee2000), Tomboulides & Orszag (Reference Tomboulides and Orszag2000). The new lower frequency is estimated to $St_{c,2}^L=0.0385$ (see figure 15 d), which lies in a reasonable range regarding the only available reference values reported at $\textit{Re}=500$ by Mittal & Najjar (Reference Mittal and Najjar1999) ( $St_2=0.05$ ), Lee (Reference Lee2000) ( $St_2=0.043$ ) and Tomboulides & Orszag (Reference Tomboulides and Orszag2000) ( $St_2=0.045$ ). It is important to note that no hysteresis was observed. This was verified by performing nonlinear simulations at a supercritical $\textit{Re}$ , followed by a gradual reduction of the control parameter. The system smoothly returned to the periodic state, confirming the supercritical nature of the third bifurcation for both the cube and sphere flows.

Finally, figure 16 presents the 3-D isocontours of the leading unstable eigenvector at $\textit{Re}=281 \gt \textit{Re}^L_{c,3}$ for the cube and at $\textit{Re}=332\gt \textit{Re}^L_{c,3}$ for the sphere. This mode is associated to the new low-frequency $St_2$ and is therefore a periodic eigenfunction, of period $\tau _2$ , represented in figure 16 by a snapshot extracted at the beginning of this period. This eigenvector, emerging in the vicinity of the UPO, leads to the Neimark–Sacker bifurcation. It exhibits complex structures that are very comparable between the cube and the sphere case.

Figure 16. Tertiary bifurcation. Three-dimensional view of isocontours of (a) streamwise $(u')$ , (b) vertical $(v')$ and (c) spanwise $(w')$ perturbation velocities plotted for the levels $\pm 0.01$ of the real part of the periodic leading unstable mode associated to the secondary low-frequency $St_2$ , at $\textit{Re} = 281$ for the cube case and $\textit{Re} = 332$ for the sphere case.

As an overview of the linear stability analyses performed in §§3, 4 and 5, table 1 presents our linear predictions for the sequence of bifurcations, detailing the type, nature and critical thresholds in the flow past a cube and a sphere (a consistency verification of our critical frequencies for Hopf (II) and Neimark–Sacker (III) bifurcations is proposed in Appendix C by comparing them with frequencies extracted from nonlinear DNS at very similar values of $\textit{Re}$ ). These results strongly support the conjectured sequence of: (I) pitchfork; (II) Hopf; and (III) Neimark–Sacker bifurcations in flows around CSC-BBs.

Table 1. Summary of the present linear stability analysis with critical Reynolds numbers ( $\textit{Re}^L_c$ ) and critical Strouhal numbers ( $St^L_c$ ). All bifurcations are of supercritical nature.

It should be noted that the critical values reported in this study are all based on the body diameter $D$ as the characteristic length scale. Normalising $\textit{Re}^L_c$ and $St^L_c$ using an alternative reference length scale could potentially collapse these critical thresholds for the cube and sphere cases, providing perhaps an even closer metric for flow behaviour across different body shapes, as suggested by attempts such as those of Gennaro, Colaciti & Medeiros (Reference Gennaro, Colaciti and Medeiros2013). However, upon closer inspection, it can easily be shown from table 1 that this normalisation process cannot rely on the same length scale for $\textit{Re}^L_c$ and $St^L_c$ values; and despite the fact that a length based on the size of the recirculation bubble appears to be a good candidate for such a normalisation procedure, our investigations found that it does not allow such unification across all bifurcations.

6. Dynamics of the quasiperiodic regime

This final section provides a detailed temporal and spatial characterisation of the flow dynamics in the quasiperiodic regime that emerges following the onset of the Neimark–Sacker bifurcation in both the cube and sphere cases. All analyses are based on DNS.

6.1. Temporal description

As previously mentioned, the quasiperiodic motion is characterised by two incommensurate fundamental frequencies $St_1$ and $St_2$ (i.e. ), and the trajectory of the quasiperiodic system is represented by a curve on a so-called $T^2$ -torus that wraps around the torus without ever exactly repeating itself. To investigate in detail the dynamic nature of a quasiperiodic flow past a cube and a sphere, we use DNS. Our analysis focuses on the vertical lift coefficient, $C_y$ , examining its time signal, power spectral density (PSD) and phase-space representation. This phase-space representation is augmented by the incorporation of the time derivative, $\partial _t C_y$ , and the drag coefficient, $C_x$ . The drag coefficient and the vertical lift coefficient are defined as follows: $C_x = {F_x}/({q\mathcal{A}}), C_y = {F_y}/({q\mathcal{A}})$ , where with being the fluid density, $U_\infty$ and $D$ the freestream velocity and the characteristic body length, and where $F_x$ and $F_y$ are the streamwise and vertical forces on the surface of the body, respectively, and $\mathcal{A}$ the body reference area equal to $\mathcal{A}=D^2$ for the cube and to $\mathcal{A}=\pi D^2/4$ for the sphere. Since the spanwise component $C_z$ is effectively zero within machine precision in the range of $\textit{Re}$ studied in this paper, we will therefore omit it from the visual representations.

It should be noted that our analysis relies on interpreting the temporal evolution of surface-integrated quantities, which, on the one hand, are relevant with respect to the global analysis performed in this work, and, on the other hand, correlate strongly with spatial probe data (e.g. velocity probes located in the wake) whose spectra show identical frequencies (see Appendix E for supporting results), thus indicating a robust synchronisation between the energy injected by the body and its redistribution by the wake for the considered $\textit{Re}$ .

Figure 17. $T^2$ -torus in the quasiperiodic (QP) regime. The $C_y$ time signal, its PSD and the $(C_x, C_y, \partial _t C_y)$ phase-space representation are shown for (a) the cube at $\textit{Re}=282$ and (b) the sphere at $\textit{Re}=332$ . In each PSD, the primary frequency $St_1$ , the secondary frequency $St_2$ , their harmonics and the frequencies resulting from nonlinear triadic interactions are labelled. (c) Schematic of the frequency decomposition (in linear scale), showing how nonlinear triadic interactions generate frequency peaks of the form $St_1 + k St_2$ ( $k \in \mathbb{Z}$ ), separated by the constant offset $St_2$ . The fundamental frequencies $St_1$ and $St_2$ are incommensurate (i.e. ).

Figures 17(a) and 17(b) depict the $C_y$ time signal, its PSD and the $(C_x, C_y, \partial _t C_y)$ phase-space representation of the flow past a cube and past a sphere at a Reynolds number slightly above the tertiary bifurcation point, that is, $\textit{Re}=282$ for the cube and $\textit{Re}=332$ for the sphere. Both signals span approximately $90$ periods of the secondary frequency $St_2$ , ensuring fine precision in the identification of the PSD peak, and equivalent spectral resolution between the cube and sphere cases. From the right-hand sides of figures 17(a) and 17(b), one can first observe that the curves drawn by the phase space representations clearly characterise a $T^2$ -torus. Concerning the Fourier spectra, one can identify two high-amplitude peaks corresponding to the fundamental frequencies $St_1$ and $St_2$ ( $St_1=0.0975$ and $St_2=0.0281$ for the cube at $\textit{Re}=282$ ; $St_1=0.1338$ and $St_2=0.0392$ for the sphere at $\textit{Re}=332$ ), as well as their respective harmonics.

The other well-defined peaks in the PSD constitute a frequency comb centred around $St_1$ , with the constant offset $St_2$ and with peak intensities decreasing in the comb as $\vert k \vert$ increases (one refers to figure 17 c for a sketch of a typical $T^2$ -torus Fourier spectrum). These frequencies arise from nonlinear triadic interactions between the two fundamental frequencies $St_1$ and $St_2$ , which satisfy the relation of the form $St_1 + k St_2$ (with integer $k \in \mathbb{Z}$ ). This process is accompanied by a corresponding nonlinear energy transfer from the primary shedding mode to these newly formed frequencies.

This spectrum analysis aligns well with observations by Bohorquez et al. (Reference Bohorquez, Sanmiguel-Rojas, Sevilla, Jiménez-González and Martínez-Bazán2011), who provided, in the context of flow past blunt-based bluff bodies, a rare analysis of the secondary low-frequency arising in the quasiperiodic regime, stating that ‘the excitation of new frequencies in the wake (i.e. the low-frequency modulation and the harmonics of the main frequency) leads to a decrease of the energy associated with the leading velocity fluctuations due to the nonlinear energy transfer from the natural frequency ( $St_1$ ) to the newly excited ones.’ This redistribution of energy suggests that the Neimark–Sacker bifurcation plays a crucial role in initiating the transition towards a more complex, eventually turbulent, state by broadening the range of excited frequencies and slowly breaking the energy balance in the wake.

6.2. Spatial description

An intuitive way to grasp the tertiary bifurcation and the ensuing quasiperiodic regime is by examining the spatial flow structures associated with the unstable Floquet mode and contrasting them with the fully nonlinear flow fields obtained from DNS. This subsection discusses how the newly emerging low-frequency motion ( $St_2$ ) combines with the primary shedding frequency ( $St_1$ ) to produce modulated vortex dynamics in the wake.

Figure 18. Superimposition of the streamwise velocity component of the periodic base flow, shedding at frequency $St_1$ (grey) with the real part of the leading unstable Floquet mode (streamwise component of velocity perturbation coloured with the levels $-0.01$ (blue) and $0.01$ (red)) represented for every quarter period of the Floquet mode’s period $\tau _2$ . One notes that while the Floquet mode covered one period $\tau _2$ , the periodic base flow covered $ St_1 /St_2 \simeq 3.5$ shedding cycles of period $\tau _1$ . (a) Cube case at $\textit{Re}=281$ [Movie 3]. (b) Sphere case at $\textit{Re}=332$ [Movie 4].

Floquet mode structure. Figure 18 illustrates, for both cube and sphere cases, the three-dimensional Floquet mode (coloured in blue and red), that becomes unstable at the tertiary bifurcation, for every quarter of its own period $\tau _2$ . The periodic base flow superimposed here in grey – already shedding at frequency $St_1$ – serves as the ‘time-periodic equilibrium’. Spatially, one sees from the 3-D mode representation that velocity perturbations appear downstream of the body and extend through the wake, concentrating mostly around the hairpin heads. Moreover, the videos given as Supplementary material in figure 18’s caption allow one to notice that, for both obstacles, the Floquet mode does not propagate at constant speed along the base flow; it accelerates between two hairpins. The unstable mode therefore manifests as a gentle modulation or ‘waviness’ applied to the hairpin-like structures characteristic of the primary shedding. Since the velocity perturbations preserve the planar symmetry of the base flow, the Floquet mode does not disrupt the underlying reflectional symmetry. Instead, it introduces a slow, global oscillation in phase, which, as supported hereafter by the nonlinear results, manifests itself as a streamwise lag modulation between the heads of successive hairpin vortices.

Quasiperiodic dynamics. Figure 19 confirms the picture in figure 18 given by Floquet theory through DNS snapshots of the wake in the supercritical regime, where both frequencies $St_1$ and $St_2$ are excited. Unlike the single-frequency shedding observed in the previous periodic regime (cf. figure 13), the flow now develops an additional low-frequency modulation under the action of $St_2$ .

Figure 19. Quasiperiodic wake states for (a,c) cube at $\textit{Re}=282$ and (b,d) sphere at $\textit{Re}=332$ , illustrating the interaction of two incommensurate frequencies. Panels (a) and (b) show timeseries of verticalvelocity $v$ , where the red markers indicate multiples of the primary period $\tau _1$ , and blue markers the secondary period $\tau _2$ . Panels (c) and (d) display isosurfaces of $\varOmega _r = 0.52$ at selected multiples of $\tau _1$ . The secondary frequency induces a convective modulation of the hairpin structures shedding (materialised by the non-equispaced grey bands), producing a quasiperiodic flow pattern that never repeats exactly. [Movie 5] (cube), [Movie 6] (sphere).

Exactly as observed by Meng et al. (Reference Meng, An, Cheng and Kimiaei2021) in the cube case, figures 19(c) and 19(d) clearly show that over approximately 7 periods of vortex shedding (i.e. $7 \tau _1$ ), the cube’s and sphere’s wakes undergo approximately 2 cycles (i.e. $2 \tau _2$ ) of shedding modulation. This observation is perfectly corroborated by figures 19(a)and 19(b) depicting for cube and sphere the time series of vertical velocity $v$ along one quasiperiod. Additionally, the time series reveal that the vortex shedding is modulated both in phase and amplitude. This phase and amplitude modulation can be seen as a slow spatial drift which is due to the superposition of two incommensurate frequencies ( $St_1$ from the original limit cycle, plus $St_2$ from the Floquet mode). This interplay produces a quasiperiodic wake, whose trajectory in phase space (e.g. lift $C_y$ versus drag $C_x$ ) wraps around a two-dimensional torus rather than closing into a single cycle (see the phase portraits in figure 17). The fact that the quasiperiodic pattern is completed within approximately 7 periods of vortex shedding or 2 cycles of wake oscillation is directly linked to the $St_1/St_2 \approx 7/2$ ratio. In the present simulations, performed at supercritical $\textit{Re}$ close to the Neimark–Sacker bifurcation, the phase oscillation of the hairpin shedding manifests solely in the $(x, y)$ plane (the flow remains planar-symmetric). However, one may hypothesise that as the $\textit{Re}$ increases to values not explored in this work, the flow could potentially develop additional three-dimensional shifts, perhaps showing a spanwise wobble in $z$ , like that observed by Sakamoto & Haniu (Reference Sakamoto and Haniu1990) and Chrust et al. (Reference Chrust, Goujon-Durand and Wesfreid2013).