2.1. Numerical simulation method and computational grids

The flow around a circular cylinder at $\textit{Re} = 200$ was computed using a numerical simulation of the incompressible Navier–Stokes equations. The governing equations are as follow:

(2.1) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0, \end{equation}

(2.2) \begin{equation} \frac {\partial {\boldsymbol{u}}}{\partial {t}}=-(\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{u}-\frac {1}{\rho }\boldsymbol{\nabla }\!{p}+\frac {1}{\textit{Re}}{{\nabla} ^2}\boldsymbol{u}, \end{equation}

where $\boldsymbol{u}$ denotes the velocity vector (with bold symbols indicating vector quantities), $p$ is the pressure, and $\rho$ is the fluid density. The Reynolds number $\textit{Re}$ is defined as

(2.3) \begin{equation} \textit{Re}=\frac {U_{\infty }D}{\nu }, \end{equation}

where $U_{\infty }$ is the free-stream velocity, $\nu$ is the kinematic viscosity, and $D$ is the diameter of the cylinder.

The governing equations were discretised using the fractional step method proposed by Le & Moin (Reference Le and Moin1991). Time advancement was performed with a third-order, three-stage Runge–Kutta scheme for the advection term and a second-order implicit Crank–Nicholson scheme for the viscous term. The time step was chosen such that the maximum Courant–Friedrichs–Lewy (CFL) number (Zang, Street & Koseff Reference Zang, Street and Koseff1994) across all grid cells remained below $0.5$ . Validation of the time-step size is provided in Appendix A.1. Spatial derivatives were computed using a second-order central difference scheme (Kajishima & Taira Reference Kajishima and Taira2017) and the quadratic upwinding interpolation for convective kinematics method (Leonard Reference Leonard1979). The pressure Poisson equation was solved using the bi-conjugate gradient stabilised method (van der Vorst Reference van der Vorst1992). Further details of the numerical implementation can be found in previous studies (Nakamura et al. Reference Nakamura, Sato and Ohnishi2024a , Reference Nakamura, Sato and Ohnishib ).

The three-dimensional computational mesh is shown in figure 1. The far-field boundary of the domain was extended to a distance of $60D$ , where $D$ is the diameter of the circular cylinder. The mesh consisted of $240$ cells in the wall-normal direction and $440$ cells in the wall-parallel direction. The height of the first layer adjacent to the cylinder surface was $1.0 \times 10^{-3}$ . A periodic boundary condition was imposed in the spanwise direction. As noted by Jiang et al. (Reference Jiang, Cheng and An2017), the characteristics of Mode A are influenced by the spanwise boundary conditions. While symmetry conditions are suitable for predicting the critical Reynolds number associated with the transition from two-dimensional to three-dimensional flow (including Mode A), periodic conditions are more appropriate for examining the three-dimensional structure of a fully developed wake. The spanwise domain length $L_z$ was varied by changing the number of grid points, while maintaining a constant spanwise grid spacing of ${\rm d}z = 0.1D$ . These grid settings and boundary conditions follow previous DNS by Jiang et al. (Reference Jiang, Cheng, Draper, An and Tong2016a , Reference Jiang, Cheng, Tong, Draper and Anb ), and Jiang & Cheng (Reference Jiang and Cheng2017).

Figure 1. Computational grid around a circular cylinder: (a) overall grids (b) close-up view of near the cylinder.

2.2. The DMD algorithm

The DMD provides a means of decomposing a flow field into a sequence of spatial modes, each associated with a characteristic frequency and growth rate. While various extensions of the DMD algorithm have been proposed (Schmid Reference Schmid2022), we introduce the standard DMD formulation, which is suitable for large-scale time-series data obtained from numerical simulations (Hemati, Williams & Rowley Reference Hemati, Williams and Rowley2014; Ohmichi, Ishida & Hashimoto Reference Ohmichi, Ishida and Hashimoto2017).

Let $\boldsymbol{u}(\boldsymbol{x},t_{\!j})$ denote the velocity field at time $t_{\!j}$ ( $j = 1, 2, \ldots , M$ ), expressed as a vector containing the three velocity components ( $u$ , $v$ , $w$ ) at all grid points, with $N$ total elements. In DMD, we construct two data matrices $X$ and $X'$ from the following snapshots:

(2.4) \begin{align} \begin{split} X=[\boldsymbol{u}(\boldsymbol{x},t_1), \boldsymbol{u}(\boldsymbol{x},t_2), {\cdots}\,, \boldsymbol{u}(\boldsymbol{x},t_{M-1})] \in \mathbb{R}^{N \times M-1},\\ X'=[\boldsymbol{u}(\boldsymbol{x},t_2), \boldsymbol{u}(\boldsymbol{x},t_3), {\cdots}\,, \boldsymbol{u}(\boldsymbol{x},t_{M})] \,\,\,\,\,\, \in \mathbb{R}^{N \times M-1}. \end{split} \end{align}

These matrices are related by a linear operator $A \in \mathbb{R}^{N \times N}$ that approximates the time evolution of the system, such that

(2.5) \begin{align} \begin{split} X' = AX. \end{split} \end{align}

The goal of DMD is to compute the eigenvalues and eigenvectors of $A$ . However, for numerically simulated flows, the size of $A$ is typically quite large, making direct computation impractical. To address this, singular value decomposition (SVD), also known as proper orthogonal decomposition (POD) (Lumley Reference Lumley1967), is applied to reduce the dimensionality of the problem:

(2.6) \begin{equation} X \approx U_rS_rV_r^{\rm T}, \end{equation}

where $U_r \in \mathbb{R}^{N \times r}$ and $V_r \in \mathbb{R}^{r \times (M-1)}$ contain the leading $r$ left and right singular vectors, respectively, and $S_r \in \mathbb{R}^{r \times r}$ is a diagonal matrix of non-negative singular values. Here, $U_r$ corresponds to the POD modes, computed without subtracting the mean flow $\boldsymbol{\bar {u}}$ defined as

(2.7) \begin{align} \begin{split} \boldsymbol{\bar {u}}(\boldsymbol{x}) = \sum ^{M-1}_{l = 1} \frac {1}{M-1}\boldsymbol{u}(\boldsymbol{x},t_l). \end{split} \end{align}

In the context of DMD, the implications of subtracting or retaining the mean field have been discussed by (Tu et al. Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2013). In this study, we apply DMD without subtracting the temporal mean from the dataset; however, the rank is determined based on the energy of the fluctuations. Specifically, the rank is chosen using the cumulative contribution ratio $R_{\textit{CCR}}$ , defined as the fraction of the total fluctuation energy captured by the retained left singular vectors. Under this criterion, the contribution of the mean-flow energy, embedded in the leading singular values, must be excluded, and the rank $r$ is determined such that

(2.8) \begin{align} \begin{split} \frac {\sum _{k=1}^{r} s^2_k}{\sum _{k=1}^{M} s^2_k} \geqslant \frac {\rho +R_{\textit{CCR}}}{\rho +1}, \end{split} \end{align}

where $s_k$ is $k$ th singular value. Here, $\rho$ represents the energy ratio between the mean field and the fluctuations:

(2.9) \begin{align} \begin{split} \rho &= \frac {(M-1)\|\boldsymbol{\bar {u}}(\boldsymbol{x})\|^2}{\sum _{l=1}^{M-1}\|\boldsymbol{u}(\boldsymbol{x},t_l) - \boldsymbol{\bar {u}}(\boldsymbol{x})\|^2}, \end{split} \end{align}

where $\| \boldsymbol{\cdot }\|$ is vector norm. Details of this criterion are provided in Nakamura et al. (Reference Nakamura, Sato and Ohnishi2024c ). To efficiently handle large datasets, we employ the incremental POD method (Ross et al. Reference Ross, Lim, Lin and Yang2008; Ohmichi Reference Ohmichi2017; Ohmichi et al. Reference Ohmichi, Ishida and Hashimoto2017) to construct $U_r$ . Furthermore, the POD is parallelised using the method of Asada & Kawai (Reference Asada and Kawai2025).

Using the reduced basis $U_r$ , we define the low-rank matrix $\tilde {X}$ as

(2.10) \begin{equation} \tilde {X}={U_r}^{\rm T}X. \end{equation}

The reduced-order approximation $\tilde {A}$ of the operator $A$ is then computed by

(2.11) \begin{equation} \tilde {A}={U_r}^{\rm T}X'\tilde {X}^{{\rm T}}(\tilde {X}\tilde {X}^{\rm T})^{\dagger }, \end{equation}

where the superscript $\dagger$ denotes the Moore–Penrose pseudoinverse. This expression is equivalent to $U_r^{\rm T} X' V_r S_r^{-1}$ , as shown in Hemati et al. (Reference Hemati, Williams and Rowley2014). Given the $k$ th eigenvector $\boldsymbol{\tilde {\varphi }}^{\textit{DMD}}_k$ of $\tilde {A}$ , the corresponding DMD mode in the full space is obtained as

(2.12) \begin{equation} \boldsymbol{\varphi }^{\textit{DMD}}_k = U_r\boldsymbol{\tilde {\varphi }}^{\textit{DMD}}_k. \end{equation}

The non-dimensional frequency $f^{\textit{DMD}}_k$ and growth rate $\sigma ^{\textit{DMD}}_k$ of this mode are computed from the associated eigenvalue $\tilde {\lambda }^{\textit{DMD}}_k$ as

(2.13) \begin{equation} f^{\textit{DMD}}_k = \frac {\text{Im}\{\log ({\tilde {\lambda }^{\textit{DMD}}_k})\}}{2\pi \Delta t}, \end{equation}

(2.14) \begin{equation} \sigma ^{\textit{DMD}}_k = \frac {\textrm{Re}\{\log ({\tilde {\lambda }^{\textit{DMD}}_k})\}}{\Delta t}, \end{equation}

where $\textrm{Re}(\boldsymbol{\cdot })$ and $\text{Im}(\boldsymbol{\cdot })$ denote the real and imaginary parts, respectively, and $\Delta t$ is the non-dimensional time interval. The function $\log (\boldsymbol{\cdot })$ denotes the complex logarithm, with the argument (phase angle) defined in $-\pi$ to $\pi$ . To enhance the spectral resolution and robustness of the DMD, we employ Hankel matrices for $X$ and $X'$ , following the approaches in Brunton et al. (Reference Brunton, Brunton, Proctor, Kaiser and Kutz2017) and Asada & Kawai (Reference Asada and Kawai2025).

2.3. The BMD algorithm

The BMD (Schmidt Reference Schmidt2020) identifies spatial structures associated with triadic interactions in statistically stationary flow fields. In such interactions, the frequencies $\{f_{\!p}, f_q, f_r\}$ satisfy a resonance condition of the form

(2.15) \begin{align} \begin{split} f_r \pm f_{\!p} \pm f_q=0. \end{split} \end{align}

Within the BMD framework, both the interaction relationships, which are characterised via the bispectrum, and their corresponding spatial structures are extracted based on the quadratic nonlinearity in the Navier–Stokes equations.

Let $\boldsymbol{u}(\boldsymbol{x},t_m)$ ( $m = 1, 2, \ldots , M$ ) denote the time series of the velocity field. To estimate an asymptotically consistent power spectrum and bispectrum, Welch’s method is employed. The time series is divided into $N_{\textit{blk}}$ segments, each containing $N_{\textit{FFT}}$ snapshots, with an overlap of $N_{{ovlp}}$ snapshots between consecutive segments.

The discrete Fourier transform of the $l$ th segment, computed numerically via the fast Fourier transform (FFT), is given by

(2.16) \begin{align} \begin{split} \hat {\boldsymbol{u}}^{l}(\boldsymbol{x},f_{\!p}) = \sum ^{N_{\textit{FFT}}-1}_{m=0}\boldsymbol{u}(\boldsymbol{x},t_m)e^{-\frac {2\pi {\rm i}pm}{N_{\textit{FFT}}}}\,\,\,\,\,\,\,(p=0,1,{\cdots} N_{\textit{FFT}}-1). \end{split} \end{align}

The sampling frequency of $\hat {\boldsymbol{u}}$ is given by $1/\Delta t$ . The bispectrum at frequencies $f_{\!p}$ and $f_q$ is defined as

(2.17) \begin{align} \begin{split} \mathcal{B}(f_{\!p},f_q) = E[\langle \hat {\boldsymbol{u}}^*(\boldsymbol{x},f_{\!p}) \circ \hat {\boldsymbol{u}}^*(\boldsymbol{x},f_q) , \hat {\boldsymbol{u}}(\boldsymbol{x},f_{\!p}+f_q)\rangle ], \end{split} \end{align}

where $\circ$ denotes the element-wise (Hadamard) product, $E[\boldsymbol{\cdot }]$ is the expectation operator, and the superscript $*$ indicates complex conjugation. In the context of discrete frequency analysis, the index $r$ corresponds to $p+q$ when the frequency triad satisfies the condition in (2.15). From the perspective of the Navier–Stokes equations, the component $\hat {\boldsymbol{u}}(\boldsymbol{x},f_{\!p}+f_q)$ arises from the quadratic nonlinear interaction between $\hat {\boldsymbol{u}}(\boldsymbol{x},f_{\!p})$ and $\hat {\boldsymbol{u}}(\boldsymbol{x},f_q)$ (Schmidt Reference Schmidt2020; Freeman et al. Reference Freeman, Martinuzzi and Hemmati2024). Thus, the bispectrum characterises the interaction relationship between the product $\hat {\boldsymbol{u}}(\boldsymbol{x},f_{\!p}) \circ \hat {\boldsymbol{u}}(\boldsymbol{x},f_q)$ (cause) and the resulting response $\hat {\boldsymbol{u}}(\boldsymbol{x},f_{\!p}+f_q)$ (effect).

The bispectrum is computed from the spectra of all $N_{\textit{blk}}$ segments. Specifically, the cause and effect within each segment are related through the coefficient $a_k(f_{\!p}, f_q)$ for $k = 1, 2, \ldots , N_{\textit{blk}}$ , as follows:

(2.18) \begin{align} \begin{split} \boldsymbol{\phi }_{f_{\!p} \circ f_q}({\boldsymbol{x},f_{\!p},f_q})&=\sum _{k=1}^{N_{\textit{blk}}}a_k(f_{\!p},f_q)\{\hat {\boldsymbol{u}}^{k}(\boldsymbol{x},f_{\!p}) \circ \hat {\boldsymbol{u}}^{k}(\boldsymbol{x},f_q)\},\\ \boldsymbol{\phi }_{f_{\!p} + f_q}({\boldsymbol{x},f_{\!p}+f_q})&=\sum _{k=1}^{N_{\textit{blk}}}a_k(f_{\!p},f_q)\{\hat {\boldsymbol{u}}(\boldsymbol{x},f_{\!p} + f_q)\}, \end{split} \end{align}

where $ \boldsymbol{\phi }_{f_{\!p} \circ f_q}({\boldsymbol{x},f_{\!p},f_q})$ and $\boldsymbol{\phi }_{f_{\!p} + f_q}({\boldsymbol{x},f_{\!p}+f_q})$ represent the cause and effect, respectively, in terms of bispectrum summarised on all segments. The core objective of BMD is to identify the coupling coefficients $a_k$ that maximise the bispectral interaction. Accordingly, $a_k$ is determined by solving the following maximisation problem:

(2.19) \begin{align} \begin{split} \boldsymbol{a} = \underset {\|{\boldsymbol{a}}\|=1} {\operatorname {argmax}} \big | E \big [\langle \boldsymbol{\phi }_{f_{\!p} \circ f_q} , \boldsymbol{\phi }_{f_{\!p} + f_q}\rangle \big ] \big |. \end{split} \end{align}

This maximisation problem is solved by computing the eigenvector corresponding to the largest eigenvalue of the matrix $B_{p,q}$ , defined as follows:

(2.20) \begin{align} \begin{split} B_{p,q}&=U_{p \circ q}^{H}\textit{WU}_{p + q} \in \mathbb{R}^{N_{\textit{blk}} \times N_{\textit{blk}}},\\ U_{p \circ q}&=[\hat {\boldsymbol{u}}^1_{p \circ q}, \hat {\boldsymbol{u}}^2_{p \circ q}, {\cdots}\, \hat {\boldsymbol{u}}^{N_{\textit{blk}}}_{p \circ q}] \in \mathbb{R}^{N \times N_{\textit{blk}}},\\ U_{p + q}&=[\hat {\boldsymbol{u}}^1_{p + q}, \hat {\boldsymbol{u}}^2_{p + q}, {\cdots}\, \hat {\boldsymbol{u}}^{N_{\textit{blk}}}_{p + q}] \in \mathbb{R}^{N \times N_{\textit{blk}}}, \end{split} \end{align}

where superscript $H$ denotes the Hermitian transpose. Here, $\hat {\boldsymbol{u}}^k_{p \circ q}$ and $\hat {\boldsymbol{u}}^k_{p+q}$ are presented as follows:

(2.21) \begin{align} \begin{split} \hat {\boldsymbol{u}}^k_{p \circ q} &= \hat {\boldsymbol{u}}^{k}(\boldsymbol{x},f_{\!p}) \circ \hat {\boldsymbol{u}}^{k}(\boldsymbol{x},f_q),\\ \hat {\boldsymbol{u}}^k_{p+q}&=\hat {\boldsymbol{u}}(\boldsymbol{x},f_{\!p}+f_q). \end{split} \end{align}

The largest eigenvalue, denoted by $\lambda ^{\textit{BMD}}_{f_{\!p},f_q}$ , represents the maximised bispectrum associated with the frequency pair $(f_{\!p}, f_q)$ . The eigenvalue and its corresponding eigenvector of $B{p,q}$ are computed using the ZGEEV routine from the LAPACK library.

The bispectrum can be evaluated at arbitrary pairs of frequencies $(f^{\textit{BMD}}_{\!p}, f^{\textit{BMD}}_q)$ . However, due to the constraint imposed by the Nyquist frequency $f_c$ , valid computations are restricted to the region defined by $-f_c\leqslant f^{\textit{BMD}}_{\!p} \leqslant f_c$ , $-f_c\leqslant f^{\textit{BMD}}_q \leqslant f_c$ and $-f_c\leqslant f^{\textit{BMD}}_{\!p} + f^{\textit{BMD}}_q \leqslant f_c$ . Moreover, the bispectral magnitude $\lambda ^{\textit{BMD}}_{f_{\!p},f_q}$ remains invariant under the exchange of $f^{\textit{BMD}}_{\!p}$ and $f^{\textit{BMD}}_q$ , as well as under complex conjugation of the bispectrum. As a result, it suffices to compute $\lambda ^{\textit{BMD}}_{f_{\!p},f_q}$ within the principal region indicated in grey in figure 2.

Figure 2. Symmetry of the bispectrum. Owing to the Nyquist frequency limit, conjugate symmetry and the symmetry between $f^{\textit{BMD}}_{\!p}$ and $f^{\textit{BMD}}_q$ , the bispectrum is evaluated only within the shaded region.

In the context of BMD, the cause distribution $\boldsymbol{\phi }_{f_{\!p} \circ f_q}$ and its corresponding effect $\boldsymbol{\phi }_{f_{\!p} + f_q}$ are referred to as the cross-frequency field and the bispectral mode, respectively. The interaction between these two fields is characterised by the following spatial mode:

(2.22) \begin{equation} \boldsymbol{\tau }_{f_{\!p} \circ f_q}(\boldsymbol{x}) = \big |\boldsymbol{\phi }_{f_{\!p} \circ f_q}(\boldsymbol{x}) \circ \boldsymbol{\phi }_{f_{\!p} + f_q}(\boldsymbol{x})\big |, \end{equation}

where $\left | \boldsymbol{\cdot }\right |$ denotes the magnitude of a complex quantity. The field $\boldsymbol{\tau }_{f_{\!p} \circ f_q}(\boldsymbol{x})$ represents the spatial distribution of interaction strength, as its integral over space equals $ |\lambda ^{\textit{BMD}}_{f_{\!p},f_q}|$ . Accordingly, $\boldsymbol{\tau }_{f_{\!p} \circ f_q}$ is referred to as the interaction map.

3.1. Numerical simulation results

Numerical simulations were conducted with a spanwise domain size of $L_z = 12D$ . According to Jiang et al. (Reference Jiang, Cheng and An2017), for $L_z \geqslant 12D$ , the time-averaged drag coefficient, the root-mean-square of the lift coefficient, and the Strouhal number ( $St$ ) associated with Kármán vortex shedding become independent of $L_z$ . Therefore, the simulation at $L_z = 12D$ can be considered effectively free from spanwise domain-size effects. Figure 3 shows an isosurface of the Q-value, coloured by the streamwise velocity component, from the simulation at $L_z = 12D$ . The wake behind the cylinder exhibits a periodic vortex structure in the spanwise direction. The observed spanwise wavelength $\lambda$ was approximately $4D$ , as three periodic structures were present across the domain. This wavelength is consistent with previous findings reported by Jiang et al. (Reference Jiang, Cheng and An2017, Reference Jiang, Cheng, Tong, Draper and An2016b).

Figure 3. Isosurfaces of the $Q$ -value at $0.1$ , coloured by the streamwise ( $x$ -direction) velocity, obtained in the computational domain with $L_z = 12D$ : (a) overall view of the flow field and (b) view in the $y$ -normal direction. The flow exhibits a three-dimensional vortex structure with a spanwise wavelength of approximately $4D$ .

As reported by Jiang et al. (Reference Jiang, Cheng and An2017), both the spanwise wavelength $\lambda$ and the temporal behaviour of the force coefficients vary with the spanwise domain length $L_z$ . To examine how temporal behaviour depends on $L_z$ , numerical simulations were performed for several values of $L_z$ . Following previous studies (Jiang et al. Reference Jiang, Cheng and An2017; Jiang & Cheng Reference Jiang and Cheng2017; Rolandi et al. Reference Rolandi, Fontane, Jardin, Gressier and Joly2023; Nakamura, Sato & Ohnishi Reference Nakamura, Sato and Ohnishi2025), the time variation of the drag and lift coefficients acting on the cylinder is used as a key indicator of the flow’s temporal characteristics. The drag coefficient $C_D$ and lift coefficient $C_L$ were computed from the simulation results using the following expression (Akhtar Reference Akhtar2008):

(3.1) \begin{align} C_D = \frac {1}{\frac {1}{2}\rho {U}_{\infty }^2DL_z}\int \int _\varOmega \left(p\cos \theta -\frac {1}{\textit{Re}}{\omega }_z\sin \theta \right)\,{\rm d}\theta {\rm d}z, \end{align}

(3.2) \begin{align} C_L = \frac {1}{\frac {1}{2}\rho {U}_{\infty }^2DL_z}\int \int _\varOmega \left(\!-p\sin \theta -\frac {1}{\textit{Re}}{\omega }_z\cos \theta \right)\,{\rm d}\theta {\rm d}z, \end{align}

where $\varOmega$ denotes the surface of the cylinder, $\theta$ is the angle measured from the stagnation point, and $\omega _z$ is the spanwise vorticity component, defined as

(3.3) \begin{equation} \omega _z = \frac {\partial v}{\partial x} - \frac {\partial u}{\partial y}. \end{equation}

The derivation of the fluid force expressions is provided in Appendix B.

Figure 4 presents the time histories of the drag coefficient for spanwise domain sizes of $L_z = 3.4D$ , $3.7D$ , $3.9D$ and $12D$ , highlighting characteristic differences in temporal behaviour. At $L_z = 3.4D$ , the drag coefficient exhibited a constant amplitude, resembling the well-known two-dimensional periodic flow around a circular cylinder (Nakamura et al. Reference Nakamura, Sato and Ohnishi2024b ), and differing markedly from the more complex dynamics observed at $L_z = 12D$ . At $L_z = 3.7D$ , small-amplitude, low-frequency fluctuations appeared in the drag signal, although their magnitude remained smaller than those observed at $L_z = 12D$ . At $L_z = 3.9D$ , more pronounced low-frequency fluctuations were evident, comparable to those at $L_z = 12D$ . As noted by Henderson (Reference Henderson1997), the presence of multiple frequency components contributes to the emergence of low-frequency fluctuations. The observed fluctuations in the present results suggest that multiple frequencies were present in the drag coefficient signals.

Figure 4. Time variation of the drag coefficient for different spanwise domain sizes: (a) $L_z = 3.4D$ , (b) $L_z = 3.7D$ , (c) $L_z = 3.9D$ and (d) $L_z = 12D$ . As $L_z$ increases, low-frequency fluctuations become more pronounced.

The time series of the drag and lift coefficients was analysed using the FFT to compute their power spectral densities. Figure 5 shows the power spectral densities for $L_z = 3.4D$ , $3.7D$ , $3.9D$ and $12D$ . The FFT was performed over a length of time window $T = 0.1 \times 4096$ , yielding a frequency resolution of $\Delta f = 1/T = 1/409.6 \approx 2.44 \times 10^{-3}$ . A Hamming window was applied to mitigate spectral leakage. For $L_z = 3.4D$ , the dominant non-zero frequency peak occurs at $f/\Delta f = 77$ , whereas for $L_z = 3.7D$ , $3.9D$ and $12D$ , the peak appears consistently at $f/\Delta f = 75$ .

Figure 5. Power spectral density of drag and lift coefficients: (a) $L_z = 3.4D$ , (b) $L_z = 3.7D$ , (c) $L_z = 3.9D$ and (d) $L_z = 12D$ . As $L_z$ increases, additional low-frequency components become apparent. Subharmonic peaks appear at $f_t$ and $f_s$ , and a low-frequency component $f_b$ emerges for $L_z = 3.9D$ , indicating the presence of triadic interactions.

It is well established that the oscillation frequency of the lift coefficient corresponds to the vortex-shedding frequency in the wake. The non-dimensional vortex-shedding frequency, known as the Strouhal number $St$ , is defined as

(3.4) \begin{equation} St = \frac {\check {f}_K D}{U_{\infty }}, \end{equation}

where $\check {f}_K$ denotes the dimensional frequency of wake vortex shedding. Although the shedding frequency can fluctuate in unsteady flows, such as those exhibiting Mode A, many previous studies (Henderson Reference Henderson1997; Posdziech & Grundmann Reference Posdziech and Grundmann2001) have estimated it using the time-averaged oscillation frequency of the lift coefficient or transverse velocity at a single point. Jiang & Cheng (Reference Jiang and Cheng2017) showed that the peak frequency obtained from the FFT of the lift-coefficient time series closely matches this time-averaged value. Based on their report, we regard the peak frequency in the FFT spectrum of the lift coefficient as the representative value of $St$ . In the present study, the observed peak frequencies, $f = 77 \Delta f \approx 0.188$ for $L_z = 3.4D$ and $f = 75 \Delta f \approx 0.183$ for $L_z = 3.7D$ , $3.9D$ and $12D$ , are in good agreement with the vortex-shedding frequencies reported by Jiang et al. (Reference Jiang, Cheng and An2017). Hereafter, for each $L_z$ , the peak non-dimensional frequency of the lift-coefficient spectrum is denoted as $f_K$ . As discussed in a later section, this frequency corresponds to the most energetic flow structure in the entire field.

For the case of $L_z = 3.4D$ , a distinct peak in the drag coefficient spectrum is observed at approximately twice the frequency of $f = f_K$ . This behaviour is characteristic of periodic flows, such as the two-dimensional flow past a circular cylinder, where the spectrum contains the fundamental frequency and its harmonics. Within the frequency range shown in the figure, only the peaks at $f_K$ and $2f_K$ are clearly visible; however, additional harmonics at $3f_K$ , $4f_K$ and higher were also identified in the higher-frequency range.

For the case of $L_z = 3.7D$ , several additional spectral peaks are observed besides the peaks at $f/\Delta f = 75$ and $f/\Delta f = 150$ , which correspond to $f_K$ and its harmonic $2f_K$ . Notably, peaks also appear at $f/\Delta f = 25$ and $f/\Delta f = 50$ . These frequencies correspond to one-third and two-thirds of the shedding frequency, respectively, indicating the presence of subharmonic components in the flow.

For the case of $L_z = 3.9D$ , a large number of additional peaks appear, although the spectral peaks at $f/\Delta f = 25$ , $50$ and $75$ are also observed. The lowest-frequency peak is identified at $f/\Delta f = 6$ in the drag coefficient spectrum. Notably, neighbouring peaks in the lift-coefficient spectrum at $f/\Delta f = 69$ and $81$ satisfy $75 - 69 = 81 - 75 = 6$ , suggesting a triadic interaction involving $f/\Delta f = 6$ and the vortex-shedding frequency $f_K$ . Similar relationships are observed for other lift-coefficient peaks: $f/\Delta f = 19$ and $31$ satisfy $25 - 19 = 31 - 25 = 6$ , and $f/\Delta f = 44$ and $56$ satisfy $50 - 44 = 56 - 50 = 6$ , further supporting the presence of triadic interactions. Additionally, a peak at $f/\Delta f = 12$ is observed in the drag coefficient spectrum, which likely corresponds to the second harmonic of $f/\Delta f = 6$ . A peak at $f/\Delta f = 18$ is also found in the drag coefficient spectrum and may be related to $f/\Delta f = 6$ as well. However, based on the current spectral data, it remains inconclusive whether this component is synchronised with the lift-coefficient peak at $f/\Delta f = 19$ .

In the following discussion, the frequency $25\Delta f$ observed for $L_z = 3.7D$ and $3.9D$ is referred to as $f_t$ , and the frequency $50\Delta f$ as $f_s$ . Additionally, the lowest frequency observed for the $L_z = 3.9D$ case, corresponding to $f = 6\Delta f \approx 0.0146$ , is denoted as $f_b$ . As discussed in a later section, frequencies, including the vortex-shedding frequency $f_K$ , may vary depending on the spanwise domain size. Therefore, the frequency labels defined here correspond specifically to the values observed for $L_z = 3.7D$ and $3.9D$ . However, similar spatial structures may appear at different frequencies in cases with other spanwise extents.

For the case of $L_z = 12D$ , clear spectral peaks are observed at $f/\Delta f = 75$ and $f/\Delta f = 150$ , corresponding to $f_K$ and its harmonic. In contrast to the case of $L_z = 3.9D$ , no distinct peaks appear at $f \lt = f_K$ . In particular, the low-frequency fluctuation evident in the time history of the drag coefficient (figure 4 d) does not manifest as a pronounced peak in the spectrum. The influence of spanwise domain size on low-frequency fluctuations is evident in the temporal variation of the drag coefficient.

Furthermore, the power spectrum density of the drag and lift coefficient at $L_z=3.7D$ and $3.9D$ captures primary peaks in $f\lt f_K$ . However, the low-frequency fluctuations of $L_z = 12D$ , as shown by the temporal variation of the drag coefficient, do not appear as evident peak frequencies in the power spectral density. For small values of $L_z$ , the flow is constrained in the spanwise direction, which limits the allowable wavelengths and, in turn, alters the temporal characteristics of the cylinder wake. A more detailed analysis is needed to quantify the frequencies associated with each $L_z$ and to clarify the relationships among the different spectral components.

3.2. The DMD-based analysis

For the case of $L_z = 12D$ , the FFT analysis does not reveal any distinct peaks in the low-frequency range ( $f \lt f_K$ ) of the power spectral density. However, low-frequency fluctuations are clearly visible in the time history of the drag coefficient. In this section, we employ DMD to investigate the existence of low-frequency components. One of the core features of DMD is that it can extract structures associated with the top- $r$ left singular vectors, that is, energetic modes that oscillate at representative growth rates and frequencies. Consequently, even when many closely spaced frequency components exist, such that FFT does not exhibit a clear peak, DMD can still identify the corresponding frequency elements, which appear explicitly as eigenvalues. This enables confirmation of their presence even when spectral leakage obscures the peaks.

To confirm the presence of the low-frequency components in cases where no clear spectral peak exists, DMD was performed for the time-series data of the flow fields in the $L_z=12D$ case. In this case, the number of snapshots was $3000$ , and the number of modes  $r$ in the low-rank approximation of SVD was $194$ , determined based on the cumulative contribution ratio exceeding $99.5 \,\%$ . Figure 6 shows the eigenvalue distribution of the DMD mode. Here, the $f_K$ value indicated in figure 6 corresponds to the eigenvalue of the DMD mode whose frequency is closest to the $f_K$ obtained from the power spectral density of the lift coefficient. Some of the DMD eigenvalues lie inside the unit circle, indicating that the associated frequencies are decaying over time. These decaying modes are found at frequencies higher than $f_K$ , whereas the low-frequency modes lie on the unit circle, indicating neutrally stable oscillations. These low-frequency components correspond to the low-frequency fluctuations in the drag coefficient shown in figure 4. The absence of distinct peaks in the power spectral density shown in figure 5(d) is therefore likely due to the presence of multiple low-frequency components with comparable energy levels, which may spread the spectral content.

The time variations of the drag coefficient, shown in figure 4, indicate that computations in a smaller domain suppress the emergence of low-frequency fluctuations. This suppression is also expected to simplify the flow dynamics due to the constraints imposed by the small spanwise size. Applying DMD to the $L_z = 3.4D$ , $3.7D$ and $3.9D$ cases may help clarify this effect. In the DMD-mode computation, the rank $r$ in the SVD approximation is truncated based on a cumulative contribution ratio exceeding $99.8\,\%$ . The truncated ranks for the $L_z = 3.4D$ , $3.7D$ and $3.9D$ cases corresponding to figure 4 are summarised in table 1. The rank $r$ increases with increasing $L_z$ , indicating that high-order modes (associated with smaller singular values) contribute non-negligibly to the flow representation. Consequently, a broader range of frequency components is captured by DMD in the larger- $ L_z$ cases. This result also implies that, as $L_z$ increases, the flow contains a greater number of distinct frequency components, making it more difficult to identify clear peaks in the power spectral densities of lift and drag coefficients shown in figure 5.

Table 1. The SVD rank used in the DMD computation, determined based on the cumulative contribution ratio exceeding 99.8 $\,\%$ .

Figure 6. Eigenvalues of DMD modes obtained from the flow field at $L_z = 12D$ . The SVD rank $r = 194$ was determined based on the cumulative contribution ratio of the singular values associated with the fluctuating component.

Figure 7 shows the eigenvalue distribution of the DMD modes for the same $L_z$ cases presented in table 1 and figure 4. For the $L_z = 3.4D$ case, no frequency component below $f^{\textit{DMD}}_k =f_K$ is observed, and the identified frequencies are uniformly spaced. This indicates that only harmonics of $f^{\textit{DMD}}_k =f_K$ are captured, which is consistent with the FFT results for both drag and lift coefficients shown in figure 5(a). In the $L_z = 3.7D$ and $3.9D$ cases, harmonics of $f_K$ are also identified. However, frequency components lower than $f^{\textit{DMD}}_k =f_K$ become significant. At $L_z = 3.7D$ , the frequencies remain uniformly spaced, similar to the $3.4D$ case. This suggests that only harmonics of the lowest frequency are present in the $L_z = 3.7D$ case. In this case, the lowest frequency is $f_t$ .

Figure 7. Eigenvalues of DMD modes obtained from flow fields at (a) $L_z = 3.4D$ , (b) $L_z = 3.7D$ and (c) $L_z = 3.9D$ . The SVD rank used in each case is listed in table 1. At $L_z = 3.4D$ , only harmonics of the shedding frequency $f_K$ are observed, whereas for $L_z = 3.7D$ and $3.9D$ , additional eigenvalues corresponding to frequencies other than the harmonics of $f_K$ appear.

For the case of $L_z = 3.9D$ , the emergence of frequency components significantly lower than $f_t$ (whose frequency is $f_K/3$ at $L_z=3.9D$ ), as revealed in the power spectral density of the lift and drag coefficients, is further supported by the distribution of DMD eigenvalues. The lowest frequency identified by the DMD analysis is $f^{\textit{DMD}}_k = 0.0149$ , which closely corresponds to the lowest peak observed in the drag coefficient spectrum. Based on this agreement, we refer to $f^{\textit{DMD}}_k = 0.0149$ as $f_b$ for the $L_z = 3.9D$ case.

Based on the eigenvalue distribution obtained from the DMD analysis, we selected the DMD mode corresponding to the shedding frequency $f_K$ for the $L_z = 3.4D$ case, and DMD modes with frequencies lower than $f_K$ for the $L_z = 3.7D$ and $3.9D$ cases. Specifically, we focus on two DMD modes, $f^{\textit{DMD}}_k = f_t$ and $f_s$ , for the $L_z = 3.7D$ case, and the DMD mode, $f^{\textit{DMD}}_k = f_b$ , for the $L_z = 3.9D$ case. Figure 8 presents the spatial distributions of the DMD modes associated with these four selected frequencies. The DMD mode at $f^{\textit{DMD}}_k = f_K$ corresponds to the classical Kármán vortex fluctuations in the wake of the cylinder. This structure is consistently observed across all $L_z$ cases, though it is omitted from the figure for brevity.

Figure 8. Real part of the streamwise ( $x$ -direction) velocity component of selected DMD modes: (a) $f^{\textit{DMD}}_k = f_K$ for the $L_z = 3.4D$ case; (b) $f^{\textit{DMD}}_k = f_s$ and (c) $f^{\textit{DMD}}_k = f_t$ for the $L_z = 3.7D$ case; and (d) $f^{\textit{DMD}}_k = f_b$ for the $L_z = 3.9D$ case. The modes in (b) and (c) are chosen for their lower frequencies relative to $f_K$ at $L_z = 3.7D$ , while the mode in (d) represents the lowest-frequency DMD mode identified for $L_z = 3.9D$ .

The DMD mode at frequency $f^{\textit{DMD}}_k = f_t$ exhibits an asymmetric structure concerning $y=0$ in the cylinder wake, in contrast to the symmetric structure observed at $f^{\textit{DMD}}_k = f_s$ . Such asymmetry relative to $y=0$ is considered to contribute to lift fluctuations. This interpretation is supported by the power spectral density shown in figure 5(b), where the asymmetric mode at $f^{\textit{DMD}}_k = f_t$ aligns with a spectral peak in the lift coefficient. In contrast, the symmetric mode at $f^{\textit{DMD}}_k = f_s$ corresponds to a peak in the drag coefficient. This tendency, where modes asymmetric concerning $y=0$ are associated with lift fluctuations, and symmetric modes with drag fluctuations, is also observed for the DMD modes at $f^{\textit{DMD}}_k = f_K$ and $f^{\textit{DMD}}_k = f_b$ shown in figure 5.

The low-frequency components (DMD modes) shown in figure 7, extracted via DMD for $L_z = 3.4D$ , $3.7D$ and $3.9D$ , may be influenced by the limited spanwise domain size. Before further analysing these low-frequency components, we first verify whether similar frequency components appear in flow fields computed within a sufficiently large domain. From the distribution of DMD eigenvalues for the $L_z = 12D$ case shown in figure 6, DMD eigenvalues corresponding to $f^{\textit{DMD}}_k \approx f_K$ , $f_s$ , $f_t$ and $f_b$ were identified. Figure 9 presents the spatial distributions of the corresponding DMD modes. These four modes exhibit the same symmetries concerning $y=0$ as their counterparts observed in the $L_z = 3.4D$ , $3.7D$ and $3.9D$ cases. This comparison suggests that the presence of these low-frequency components is not a numerical artefact resulting from the limited spanwise extent, but rather a distinct feature of the flow dynamics.

Figure 9. The DMD modes corresponding to the same frequencies as in figure 8, but for the $L_z = 12D$ case. These distinct frequency components appear in the fully developed Mode A regime and occur at frequencies lower than $f_K$ .

To investigate the temporal evolution of the flow in detail, numerical simulations were performed using finely adjusted spanwise domain sizes: $L_z = 3.2D$ , $3.5D$ , $3.6D$ , $3.8D$ , $4.0D$ , $4.2D$ , $4.7D$ and $5.0D$ . For $L_z = 3.2D$ and $5.0D$ , the flow remained two-dimensional, as indicated by the negative growth rate of the spanwise-velocity component. The onset of three-dimensional flow occurred at values of $L_z$ consistent with the stable spanwise-wavelength range of Mode A, as predicted by Floquet analysis (Barkley & Henderson Reference Barkley and Henderson1996; Rolandi et al. Reference Rolandi, Fontane, Jardin, Gressier and Joly2023). This range corresponds to a regime in which spanwise perturbations remain sufficiently small that nonlinear effects in the spanwise velocity can be neglected. This suppression arises because, under the constraint of a limited spanwise extent, only wavelengths with negative growth rates are permitted during the linear development from the two-dimensional flow. As a result, the onset of three-dimensionality is inhibited, and low-frequency components do not emerge. This relationship between spanwise domain size and the onset of three-dimensionality in numerical simulations has also been reported by Jiang et al. (Reference Jiang, Cheng and An2017).

The DMD was applied to the time-series flow field data for spanwise domain sizes of $L_z = 3.4D$ , $3.5D$ , $3.6D$ , $3.7D$ , $3.8D$ , $4.0D$ , $4.2D$ , $4.3D$ and $4.7D$ . The cases with $L_z = 3.7D$ and $3.8D$ represent transitional regimes, where the periodic characteristics of the flow begin to weaken, and low-frequency DMD modes gradually emerge, making relatively small contributions. To capture these weak components in the $L_z = 3.7D$ and $3.8D$ cases, we adopted a fixed number of DMD modes, $r = 194$ , based on the rank determined for the $L_z = 12D$ case. For the other cases, the rank was selected such that the cumulative energy contribution exceeded 99.8 %. From the resulting DMD eigenvalue distributions, we extracted neutrally stable modes (i.e. zero growth rate) within the frequency range $0 \lt f^{\textit{DMD}}_k \leqslant 0.2$ to isolate the low-frequency dynamics. Figure 10 presents the DMD-mode frequencies $f^{\textit{DMD}}_k$ obtained for each spanwise domain size. Frequencies marked with red triangles indicate the $f_K$ associated with vortex-shedding frequency, which closely match the values reported by Jiang et al. (Reference Jiang, Cheng and An2017), shown as blue triangles.

Figure 10. The DMD-mode frequencies $f^{\textit{DMD}}_k$ in the range of $0$ to $0.2$ for various spanwise domain sizes: (a) $3.2D \leqslant L_z \leqslant 5.0D$ ; (b) $L_z = 12D$ . Two notable transitions are observed, the second of which marks the onset of multiple low-frequency DMD modes.

The first notable change in the emergence of low-frequency DMD modes occurs between $L_z = 3.5D$ and $3.6D$ , marked by the appearance of frequencies at $f^{\textit{DMD}}_k = f_t$ and $f_s$ . Between this transition and the next, which occurs between $L_z = 3.7D$ and $3.8D$ , only three dominant frequency components, $f^{\textit{DMD}}_k = f_t$ , $f_s$ and $f_K$ , are observed. Following the transition at $L_z = 3.8D$ , frequency components lower than $f_t$ begin to appear, for example, in the case of $L_z = 3.9D$ , this corresponds to $f^{\textit{DMD}}_k = f_b$ , along with multiple components clustered around $f^{\textit{DMD}}_k = f_K$ , $f_t$ and $f_s$ . As $L_z$ increases further, the number of distinct frequencies near $f^{\textit{DMD}}_k = f_K$ , $f_s$ and $f_t$ gradually increases. These observations suggest that the emergence of DMD modes with $f^{\textit{DMD}}_k \lt f_t$ marks the onset of multi-frequency behaviour in the flow field.

Figure 11 presents the frequencies identified in figure 10, normalised by the vortex-shedding frequency $f_K$ corresponding to each value of $L_z$ . As discussed, for $3.6D \leqslant L_z \leqslant 3.7D$ , DMD modes corresponding to $f^{\textit{DMD}}_k = f_t$ are clearly observed. In the range $3.8D \leqslant L_z \leqslant 4.0D$ , components with $f^{\textit{DMD}}_k/f_K \lt 1/3$ appear alongside those at $f^{\textit{DMD}}_k = f_t$ . However, as $L_z$ increases further, the frequencies increasingly deviate from the $f^{\textit{DMD}}_k = f_t$ line. This is attributed to the emergence of additional low-frequency components in the flow, which interact with dominant frequencies such as $f^{\textit{DMD}}_k = f_t$ , $f_s$ and $f_K$ , leading to frequency deviation. Despite this deviation, the DMD modes with frequencies near $f_t$ and $f_s$ (shown in figures 9 b and 9 c) retain symmetry concerning $y=0$ , similar to that observed in the $L_z = 3.9D$ case. This similarity suggests that the $L_z = 3.7D$ and $3.9D$ cases provide representative examples for analysing the behaviour of these frequency components, offering insight into the underlying dynamics even when the spanwise domain is sufficiently large.

Figure 11. The DMD-mode frequencies $f^{\textit{DMD}}_k$ , normalised by the shedding frequency $f_K$ , in the range $0$ to $0.2$ for various spanwise domain sizes: (a) $3.2D \leqslant L_z \leqslant 5.0D$ ; (b) $L_z = 12D$ .

3.3. The BMD-based analysis

Previous studies have suggested that the low-frequency fluctuations in Mode A arise from the presence of multiple frequencies that are close to, but distinct from, the shedding frequency. In the case of $L_z = 3.9D$ , the sum and difference of $f_K$ and $f_b$ yield frequencies near $f_K$ but offset from it. Furthermore, the peak frequencies observed in the power spectral densities of the lift and drag coefficients at $L_z = 3.9D$ indicate the presence of components at $f = f_K \pm f_b$ , $f_t \pm f_b$ and $f_s \pm f_b$ (see figure 5). These observations suggest that combinations of frequency components coincide with other frequencies, implying the presence of nonlinear interactions (Phillips Reference Phillips1960). To explore these interaction mechanisms in Mode A, we apply BMD to the flow fields for the $L_z = 3.7D$ and $3.9D$ cases.

The BMD was applied to the time-series data for the $L_z = 3.7D$ case. Details regarding the validation of FFT parameters used in the BMD algorithm are provided in Appendix C. Figure 12 shows the mode bispectrum obtained from the BMD analysis. Spectral peaks appear at intervals of $f_t$ , consistent with the frequency patterns identified in the power spectral densities (figure 5) and DMD results (figure 10). The appearance of a bispectral peak confirms triadic interactions among the three frequency components $f = f_K$ , $f_s$ and $f_t$ . The presence of only integer multiples of $f_t$ across the entire spectrum indicates that the harmonics of the lowest-frequency component coincide with the shedding frequency $f_K$ . In other words, in the case of $L_z=3.7D$ , the lowest-frequency component is the third subharmonic, exactly equal to $f_t$ . The frequency pattern shown in figure 10 suggests that the number `3’ is a characteristic feature of the system within the range $3.5D \lt L_z \lt 3.8D$ . The appearance of only doublet waves based on the third subharmonic indicates that the flow field at $L_z = 3.7D$ exhibits overall periodic behaviour. This is further supported by the time variation of the drag coefficient in figure 4(b), which displays a repeating pattern with three distinct peaks.

Figure 12. Mode bispectrum corresponding to the two frequencies $f^{\textit{BMD}}_{\!p}$ and $f^{\textit{BMD}}_q$ , obtained via BMD for the case $L_z = 3.7D$ : (a) wide frequency range; (b) close-up of the range $0$ to $f_K$ . The absolute values of the bispectrum are displayed using a logarithmic colour scale. Nonlinear interactions are observed only at harmonics of $f_t$ . The spatial structures of modes (A) and (B) are shown in figure 13.

High peaks are observed along the $0$ -frequency line $f^{\textit{BMD}}_q = 0$ , as well as at points associated with the $f_K$ -frequency component, such as ( $f^{\textit{BMD}}_{\!p}$ , $f^{\textit{BMD}}_q$ ) = ( $f_K$ , $f_K$ ), ( $2f_K$ , $f_K$ ), ( $f_K$ , $-f_K$ ) and ( $2f_K$ , $-f_K$ ). These features indicate an energy cascade into harmonics of the $f_K$ -frequency component. Along the $f^{\textit{BMD}}_{\!p} = f_K$ line, the interactions between $f_K$ - and $f_t$ -frequency components are more pronounced than those between $f_K$ - and $f_s$ -frequency components, suggesting a stronger nonlinear interaction between the Kármán vortex fluctuation and the $f_t$ -frequency component. This implies a possible relationship between the Kármán vortex fluctuation and the emergence of the $f_t$ -frequency component.

Figure 13 presents the spatial distributions and interaction relationships of the $(f^{\textit{BMD}}_{\!p}, f^{\textit{BMD}}_q) = (f_s, f_t)$ and $(f_K, -f_t)$ modes, which exhibit strong triadic interactions among the frequencies $f_t$ , $f_s$ and $f_K$ . In the figure, each arrow points in the direction of the resulting frequency component, given by $f^{\textit{BMD}}_{\!p} + f^{\textit{BMD}}_q$ . In the $(f_s, f_t)$ interaction, the bispectral mode $\phi _{f_s + f_t}$ corresponds to a Kármán vortex fluctuation with a sum of frequency $f^{\textit{BMD}}_{p+q} = f_K$ . For the $(f_K, -f_t)$ interaction, the bispectral mode $\phi _{f_K - f_t}$ closely matches the spatial structure of the DMD mode shown in figure 8(b). Because bispectral modes of periodic flows coincide with Fourier modes, and because DMD and Fourier modes at the same frequency exhibit similar spatial structures (Tu Reference Tu2013), the bispectral and DMD modes are expected to share the same spatial distribution. Accordingly, the interactions captured by BMD are consistent with those identified by DMD. The interaction map shows a coherent spatial structure across the three frequency components, particularly in the downstream region of the cylinder, where aligned features in the streamwise ( $x$ ) direction emphasise the strength of their triadic interaction.

Figure 13. Isosurfaces of bispectral mode and interaction map obtained for the $L_z = 3.7D$ case. All isosurfaces represent the streamwise ( $x$ -direction) velocity component. The BMD results indicate a triadic interaction among the three frequency components $f = f_t$ , $f_s$ and $f_K$ . The central mode in the figure corresponds symbolically to the DMD mode at frequency $f_t$ shown in figure 8(c).

The three frequency components at which bispectral peaks were identified using BMD are, from a fluid dynamics perspective, energetically coupled through the nonlinear terms of the governing equations. According to the formulation by Freeman et al. (Reference Freeman, Martinuzzi and Hemmati2024), the magnitude and direction of energy transfer resulting from nonlinear interactions identified by BMD are given by

(3.5) \begin{equation} E_{f_{\textit{dn}} \stackrel {f_{\textit{ct}}}{ \rightarrow } f_{\textit{rc}}} = -\textrm{Re}\!\left [ \int \boldsymbol{\hat {u}}^H(\boldsymbol{x},f_{\textit{rc}})\left \{\boldsymbol{\hat {u}}(\boldsymbol{x},f_{\textit{ct}})\boldsymbol{\cdot }\boldsymbol{\nabla }\right \}\boldsymbol{\hat {u}}(\boldsymbol{x},f_{\textit{dn}}) {\rm d}\boldsymbol x\right ]\!, \end{equation}

where the frequencies $f_{\textit{rc}}$ , $f_{\textit{ct}}$ and $f_{\textit{dn}}$ satisfy the triadic relation $f_{\textit{rc}} = f_{\textit{ct}} + f_{\textit{dn}}$ . The term $E_{f_{\textit{dn}} \stackrel {f_{\textit{ct}}}{ \rightarrow } f_{\textit{rc}}}$ quantifies the energy transferred from the donor frequency $f_{\textit{dn}}$ to the recipient frequency $f_{\textit{rc}}$ via the catalyst frequency $f_{\textit{ct}}$ . This formulation, known as the recipient–donor framework, implies that the donor frequency $f_{\textit{dn}}$ transfers energy to the recipient frequency $f_{\textit{rc}}$ through the catalyst frequency $f_{\textit{ct}}$ , which facilitates the nonlinear interaction without directly receiving or donating energy itself.

For the frequencies at which bispectral peaks were observed, and triadic interactions were inferred, the energy transfer terms defined in (3.5) were computed. Figure 14 presents the energy transfer terms for selected recipient frequencies, $f_{\textit{rc}} = f_K$ , $2f_K$ , $f_t$ and $f_s$ , corresponding to prominent bispectral peaks. The horizontal axis indicates the donor frequency $f_{\textit{dn}}$ , while the catalyst frequency is determined as $f_{\textit{ct}} = f_{\textit{rc}} - f_{\textit{dn}}$ . The component at the vortex-shedding frequency $f_K$ gains energy from the $0$ -frequency mode via nonlinear transfer. In contrast, it loses energy to the second harmonic ( $2f_K$ ), suggesting a transfer of energy from $f_K$ to $2f_K$ . This interpretation is supported by figure 14(b), where the $f_{\textit{rc}} = 2f_K$ component shows a positive energy gain from $f_K$ . Additionally, as shown in figure 14(a), the $f_K$ component exhibits negative energy transfer to itself. This corresponds to the case where $f_{\textit{ct}} = f_K - f_K = 0$ , representing convection by the $0$ -frequency component. This result implies that energy at $f_K$ is being convected out of the computational domain by the $0$ -frequency component. Similar negative self-convection has also been reported by Yeung et al. (Reference Yeung, Chu and Schmidt2024).

Focusing on frequencies lower than the vortex-shedding frequency $f_K$ , the $f_t$ -frequency component is found to receive the most significant amount of energy from the $f_K$ -frequency component. It also gains energy through a difference interaction between the $0$ - and $f_K$ -frequency components. Additionally, $f_t$ -frequency component transfers energy to the $-f_s$ -frequency component, with the $f_K$ -frequency component serving as a catalyst for this transfer term. The $f_s$ -frequency component receives most of its energy from the $0$ -frequency component, followed by the $-f_t$ -frequency component. In the latter case, the energy transfer is again facilitated by the $f_K$ -frequency component. This suggests that the $f_s$ -frequency component does not arise from a self-interaction of the $f_t$ -frequency component (i.e. $2 \times f_t$ ), but rather from a difference interaction of the form $f_K - f_t$ . In other words, unless the $f_t$ -frequency component is an exact third subharmonic of the $f_K$ -frequency component, the $f_s$ -frequency component cannot be regarded as the second harmonic of the $f_t$ -frequency component. For instance, if the interacting frequency is $f_K/4$ , the difference $f_K - (f_K/4) = 3f_K/4$ does not correspond to a second harmonic of $f_K/4$ . This interpretation is further supported by the weak bispectral peak at $(f^{\textit{BMD}}_{\!p}, f^{\textit{BMD}}_q) = (f_t, f_t)$ , indicating that the self-interaction of the $f_t$ -frequency component is weak. Furthermore, in the $L_z = 4.0D$ case shown in figure 11, a DMD-derived frequency slightly below $f_K/3$ and another slightly above $2f_K/3$ are observed. This can be interpreted as the result of a difference interaction between $f_K$ and a component slightly lower than $f_K/3$ , denoted as $f_K/3 - \delta$ . The resulting frequency, $f_K - (f_K/3 - \delta ) = 2f_K/3 + \delta$ , lies slightly above $2f_K/3$ , consistent with the DMD eigenvalue distribution.

Figure 14. Energy transfer terms associated with various $f_{\textit{dn}}$ at $L_z = 3.7D$ : (a) $f_{\textit{rc}} = f_K$ ; (b) $f_{\textit{rc}} = 2f_K$ ; (c) $f_{\textit{rc}} = f_t$ ; (d) $f_{\textit{rc}} = f_s$ .

We performed BMD for the $L_z = 3.9D$ case to investigate interactions associated with the $f_b$ -frequency component. Figure 15 shows the distribution of the mode bispectrum obtained from the BMD analysis. Compared with the $L_z = 3.7D$ case, where no lowest-frequency component other than $f_t$ was identified, a greater number of bispectral peaks were observed. Primary peaks appear at $f_t$ , $f_s$ and $f_K$ , similar to the $L_z = 3.7D$ case. However, for $L_z = 3.9D$ , additional satellite peaks are observed around these primary frequencies, a feature not present in the earlier case. Along the line $f^{\textit{BMD}}_q = f_b$ , distinct peaks appear at $(f^{\textit{BMD}}_{\!p}, f^{\textit{BMD}}_q) = (f_b, f_t)$ , $(f_b, f_s)$ and $(f_b, f_K)$ , as highlighted by the white circles in figure 15(b). These correspond to interactions between the $f_b$ -frequency component and the $f_t$ -, $f_s$ -, and $f_K$ -frequency components. These satellite peaks thus indicate nonlinear interactions between the $f_b$ -frequency component and the other primary frequency components $f^{\textit{BMD}}_{\!p} = f_t$ , $f_s$ and $f_K$ .

Figure 15. Mode bispectrum corresponding to the two frequencies $f^{\textit{BMD}}_{\!p}$ and $f^{\textit{BMD}}_q$ , obtained via BMD for the $L_z = 3.9D$ case: (a) wide frequency range; (b) close-up of the range $0$ to $f_K$ . The absolute values of the mode bispectrum are presented using a logarithmic colour scale. With the emergence of the $f_b$ -frequency component, a satellite peak appears near the $f_K$ -, $f_t$ - and $f_s$ -frequency peaks. The interaction points labelled (C), (D), (E) and (F), representing interactions between the $f_K$ component and the $f_b$ -frequency component, are shown in figures 16 and 17.

In the $L_z = 3.9D$ case, the spectral magnitudes of satellite peaks such as ${f_{\!p}^{\textit{BMD}}} = f_t + 2f_b$ , $f_t - 2f_b$ and $f_K - 2f_b$ , which result from interactions with the $2f_b$ -frequency component, were smaller than those of peaks at $f_t + f_b$ , $f_t - f_b$ and $f_K - f_b$ , which arise from interactions with the $f_b$ -frequency component. This suggests that the nonlinear interaction involving the $2f_b$ -frequency component is weaker than that involving the $f_b$ -frequency component. However, as shown in figure 10, the number of frequency components identified by the DMD increases with the spanwise domain size $L_z$ . This indicates that frequency components initially associated with weak peaks become detectable by the DMD. Therefore, in a flow field with a sufficiently large spanwise extent, nonlinear interactions involving the lowest-frequency components enhance the emergence of additional frequency components throughout the flow.

Figure 16. Isosurfaces of bispectral modes and interaction maps at $(f^{\textit{BMD}}_{\!p}, f^{\textit{BMD}}_q) = (f_K, f_b)$ and $(f_K^+, -f_b)$ , obtained for the $L_z = 3.9D$ case. All isosurfaces represent the streamwise ( $x$ -direction) velocity component. The central mode in the figure corresponds symbolically to the DMD mode at frequency $f_b$ shown in figure 8(d).

We examined the mode distributions at the interaction points where large peaks were observed in figure 15. For simplicity, the frequencies $f_K + f_b$ and $f_K - f_b$ are hereafter denoted as $f_K^+$ and $f_K^-$ , respectively. Figure 16 presents the interaction relationships and spatial distributions of the modes obtained from BMD at $(f^{\textit{BMD}}_{\!p}, f^{\textit{BMD}}_q) = (f_K, f_b)$ and $(f_K^+, -f_b)$ . These correspond to interactions mediated by the $f_b$ -frequency component between the primary peak and its satellite peaks. The bispectral modes $\phi _{f_K^+ - f_b}$ and $\phi _{f_K+f_b}$ both represent Kármán vortex fluctuations in the wake at different frequencies. Hence, the presence of several shedding frequencies is associated with low-frequency fluctuations, which supports the assertions of Henderson (Reference Henderson1997) and Jiang et al. (Reference Jiang, Cheng and An2017). The interaction map indicates a strong interaction near the centreline of the wake ( $y = 0$ ), particularly close to the cylinder. This localised interaction implies that the Kármán vortex-shedding frequency in the wake is determined primarily by the shedding position or timing at the cylinder.

Figure 17. Isosurfaces of cross-frequency fields, bispectral modes and interaction maps at $(f^{\textit{BMD}}_{\!p}, f^{\textit{BMD}}_q) = (f_b, -f_b)$ and $(0, f_b)$ , obtained for the $L_z = 3.9D$ case. All isosurfaces represent the streamwise ( $x$ -direction) velocity component.

Figure 17 shows the spatial distributions of the modes obtained from BMD at $(f^{\textit{BMD}}_{\!p}, f^{\textit{BMD}}_q) = (f_b, 0)$ and $(f_b, -f_b)$ , along with their interaction relationships. The interactions involving the $0$ -frequency (mean flow) and $f_b$ -frequency component, namely $(f^{\textit{BMD}}_{\!p}, f^{\textit{BMD}}_q) = (f_b, 0)$ and $(f_b, -f_b)$ , exhibit a symmetric structure about $y = 0$ . The cross-frequency fields $\phi _{f_b \circ 0}$ and $\phi _{f_b \circ (-f_b)}$ are localised near the wake recirculation region associated with the $0$ -frequency component, represented by the bispectral mode $\phi _{f_b - f_b}$ . The interaction map also indicates strong coupling near the recirculation region. These distributions suggest a correlation between the vortical structures formed near the recirculation region and the $f_b$ -frequency component.

Figure 18 presents the energy transfer terms computed for various recipient frequencies at $L_z = 3.9D$ . The recipient frequency labels $f_{\textit{rc}}$ shown in the figure are based on the frequencies identified from the FFT results. In figure 18(a), the energy transfer characteristics at $f_{\textit{rc}} = f_t$ display similar directions and magnitudes to those observed at the same frequency in the $L_z = 3.7D$ case. This indicates that the energy transfer relationships among the $f_t$ -, $f_s$ - and $f_K$ -frequency components identified for $L_z = 3.7D$ also hold for $L_z = 3.9D$ .

From the energy transfer terms of $f_{\textit{rc}} = f_b$ shown in figure 18(b), it is evident that the $f_b$ -frequency component receives the largest amount of energy from $f_K^+$ , a frequency slightly offset from the vortex-shedding frequency $f_K$ . As shown in figure 18(c), the energy transfer terms of $f_{\textit{rc}} = f_K^+$ indicate that the $f_K^+$ -frequency component, in turn, receives most of its energy from the $f_K$ -frequency component. This suggests that a slight shift in the dominant vortex-shedding frequency facilitates energy transfer to the lowest-frequency component of $f_b$ . Similarly, figure 18(d) shows that the $f_K^-$ -frequency component also receives its largest energy input from the $f_K$ -frequency component. However, in this case, there is minimal energy transfer toward the low-frequency range.

We now focus on frequency components that are considered to result from nonlinear interactions involving the lowest-frequency component $f_b$ . As shown in figure 18(e), the energy transfer terms of the $f_{\textit{rc}} = 2f_b$ -frequency component indicate that this component receives energy from the $f_b$ -frequency component indicated by the blue square in the figure, supporting the interpretation that the $2f_b$ -frequency component is a harmonic of the $f_b$ -frequency component. However, the largest amount of energy is transferred from the frequency component of $f_K + 2f_b = f_K^{+} + f_b$ , which is highlighted with a red circle in the figure and denoted as $f_K^{++} = f_K^{+} + f_b$ . The energy transfer terms of $f_{\textit{rc}} = f_K^{++}$ , shown in figure 18(f), reveals that the $f_K^{++}$ -frequency component receives energy from both the $f_K$ - and $f_K^+$ -frequency components. This suggests that the $f_K^{++}$ -frequency component emerges as a result of a slight shift in the vortex-shedding frequency.

Figure 18. Energy transfer terms associated with various $f_{\textit{dn}}$ at $L_z = 3.9D$ : (a) $f_{\textit{rc}} = f_t$ ; (b) $f_{\textit{rc}} = f_b$ ; (c) $f_{\textit{rc}} = f_K^+$ ; (d) $f_{\textit{rc}} = f_K^-$ ; (e) $f_{\textit{rc}} = 2f_b$ ; (f) $f_{\textit{rc}} = f_K^{++}$ .

Finally, we confirm that the interaction relationship between the low-frequency components and the vortex-shedding frequency component is preserved even in a sufficiently large spanwise domain. Figure 19(a) shows the energy transfer terms for various $f_{\textit{dn}}$ when the recipient frequency is set to $f_{\textit{rc}} = 51\Delta f$ . To improve the convergence of the transfer term, the energy transfer was computed after performing a spanwise average every $4D$ . The results indicate that the $f_s$ -frequency component receives energy from the $f_K$ -frequency component and that it also receives energy from the $f_t$ -frequency component through the difference interaction involving $f_t$ and $f_K$ . These interaction patterns are consistent with those observed for the $L_z = 3.7D$ domain.

Figure 19(b) shows the magnitude of the energy transfer for various $f_{\textit{dn}}$ when the recipient frequency is set to $f_{\textit{rc}} = 6\Delta f \approx f_b$ . The results indicate that the $f_b$ -frequency component receives the largest amount of energy from the $f_K^{+}$ -frequency component. Moreover, for frequencies near $f_K$ and $-f_K$ , even slight variations in $f_{\textit{dn}}$ cause frequent sign changes in the transfer direction, implying that $f_b$ and the frequency components around the vortex-shedding frequency exchange energy with one another. This behaviour is consistent with the trend observed for the $L_z = 3.9D$ case, shown in figure 18(b).

Based on this discussion, the low-frequency components can be classified into two groups from the perspective of energy transfer. The first group comprises components such as the $f_t$ - and $f_s$ -frequency components, which receive energy directly from the dominant vortex-shedding mode at frequency $f_K$ . The second group includes components such as the $f_b$ -frequency component, which gain energy from frequency components slightly offset from $f_K$ , such as $f_K^{++}$ and $f_K^+$ . Notably, components in the second group tend to lie at lower frequencies than those in the first group, contributing to the emergence of additional frequency components through harmonic generation and nonlinear interactions with other low-frequency components.

Table 2. Harmonic structure for various $L_z$ values, characterised by the ratios $f_K/g$ and $f_K/s$ , where $s$ denotes the frequency near one-third of $f_K$ , and $g$ is the lowest frequency.

Figure 19. Energy transfer terms associated with various $f_{\textit{dn}}$ at $L_z = 12D$ : (a) $f_{\textit{rc}} = f_s$ ; (b) $f_{\textit{rc}} = f_b$ .

4.1. Harmonic nature between the lowest-frequency mode and Kármán vortex

The BMD results suggest that the $f_t$ -frequency component acts as the third subharmonic of the vortex-shedding frequency for $L_z = 3.7D$ and $3.9D$ . However, further analysis is required to determine whether the $f_b$ -frequency component can be regarded as a subharmonic of $f_K$ , particularly for other spanwise domain sizes. Table 2 presents the dominant Kármán vortex-shedding frequency $f_K$ , the lowest frequency $g$ obtained from the DMD analysis, and the corresponding ratio $f_K/g$ . It also lists the frequency $f'_t(L_z)$ closest to one-third of $f_K$ for each spanwise domain width $L_z$ . As $L_z$ increases, the value of $f_K/f'_t$ gradually deviates from an exact subharmonic relation. Notably, the appearance of a lowest frequency that is not a subharmonic of $f_K$ occurs only when $L_z \geqslant 3.8D$ . For $L_z = 3.8D$ and $3.9D$ , the value of $f_K/f'_t$ remains close to $3$ , with deviations at the decimal level being negligibly small. In contrast, for $L_z = 3.6D$ and $3.7D$ , the decimal part of $f_K/g$ is one order of magnitude smaller than in the cases where the lowest-frequency other than frequency $f'_t$ is present. This indicates that $f_K/g$ cannot be regarded as an integer, especially at larger $L_z$ , suggesting that the $f_K$ -frequency component is not harmonically related to the lowest-frequency component.

Previous studies have shown that fully developed Mode A gives rise to a quasi-periodic flow field. This behaviour is grounded in the bifurcation theory of dynamical systems, where quasi-periodicity arises from the presence of multiple oscillators with incommensurate (i.e. irrationally related) periods. Focusing again on the energy transfer into the low-frequency component $f_b$ , which contributes to the breakdown of periodicity, we observe the appearance of frequencies slightly offset from the dominant vortex-shedding frequency, such as $f_K^+$ . The nonlinear interaction between these components, specifically $f_K^+ - f_K$ , supplies energy to the $f_b$ -frequency component. In this sense, the emergence of the $f_b$ -frequency component suggests the presence of multiple vortex-shedding oscillators.

These oscillators are expected to have incommensurate periods, at least when the spanwise domain is sufficiently wide. Consequently, the $f_b$ -frequency component, maintained through their difference interaction, can be interpreted as an irrational multiple of the vortex-shedding frequency. While any frequency computed numerically must be rational due to finite resolution and truncation, the fact that the values of $f_K/g$ reported in table 2 for $L_z \geqslant 3.8D$ are non-integer supports the plausibility of quasi-periodicity in the system.

4.2. Relationship between the $0$ -frequency component and the $f_b$ -frequency component

Based on the BMD and energy transfer analysis, the $f_b$ -frequency component observed at $L_z = 3.9D$ differs from other components, such as $f_t$ - and $f_s$ -frequency components, in that it does not receive energy directly from the $f_K$ -frequency component. This characteristic highlights $f_b$ as a particularly distinct feature within the low-frequency spectral composition.

Returning to the self-interaction map of the $f_b$ -frequency component shown in figure 17, we observe that the interaction region is localised around the wake recirculation region associated with the $0$ -frequency component. This spatial overlap suggests a potential relationship between the $f_b$ -frequency component and the dynamics of the wake recirculation. Figure 20 presents a cross-section at $z = 0$ for the DMD mode with $f_k^{\textit{DMD}} = 0$ and an isosurface of the DMD mode with $f_k^{\textit{DMD}} = f_b$ , obtained from the $L_z = 3.9D$ case. The symmetric structure of the DMD mode with $f_k^{\textit{DMD}} = f_b$ , concerning $y = 0$ is similar to that of the $0$ -frequency DMD mode. Furthermore, the $f_b$ -frequency DMD mode is distributed along the wake recirculation region of the $0$ -frequency DMD mode. The interaction maps and spatial distribution of the cross-frequency fields from the BMD, as shown in figure 17, further support these spatial correlations. These findings indicate that the $f_b$ -frequency DMD mode ( $f_b$ -frequency component) is closely linked to the distribution of the recirculation region formed behind the cylinder.

Figure 20. Overlay plot at $z = 0$ showing the real part of the streamwise ( $x$ -direction) velocity component of the DMD mode with $f_k^{\textit{DMD}} = 0$ , along with the isosurface of the DMD mode at $f_k^{\textit{DMD}} = f_b$ , obtained for the $L_z = 3.9D$ case.

Numerous studies on the recirculation region in cylinder wakes (Fornberg Reference Fornberg1980; Nakamura et al. Reference Nakamura, Sato and Ohnishi2025) have identified the distance between the cylinder surface and the zero-crossing point of the streamwise velocity component at $y = 0$ as a key parameter characterising the extent of the recirculation region. Figure 21 presents a conceptual diagram of the recirculation region formed behind the cylinder. In this study, the recirculation length, $L_{\textit{recirc}}$ , is defined as the $x$ -coordinate at which the spanwise-averaged streamwise velocity becomes zero, excluding the no-slip region on the cylinder surface.

Figure 21. Definition of $L_{\textit{recirc}}$ .

Figure 22 presents the recirculation length $L_{\textit{recirc}}$ as a function of spanwise length $L_z$ , along with its gradient concerning $L_z$ . The gradient was calculated using central differencing between two reference points, $L_z = L_1$ and $L_2$ , and evaluated at the midpoint $L_z = (L_1 + L_2)/2$ as follows:

(4.1) \begin{equation} \frac {{\rm d}L_{\textit{recirc}}}{{\rm d}L_z}\bigg |_{L_z = \frac {L_1 + L_2}{2}} = \frac {L_{\textit{recirc}}(L_1) - L_{\textit{recirc}}(L_2)}{L_1 - L_2}. \end{equation}

The grey lines in figure 22 highlight that the $f_b$ -frequency component appears for $L_z \geqslant 3.8D$ , consistent with the frequency distribution obtained from the DMD results shown in figure 10.

For $L_z \lt 3.8D$ , the recirculation length $L_{\textit{recirc}}$ increased with spanwise length, indicating that the recirculation region expanded as $L_z$ grew. This trend ceased at $L_z = 3.8D$ , as clearly evidenced by the vanishing gradient of $L_{\textit{recirc}}$ concerning $L_z$ . Beyond this point, the recirculation length remained nearly constant, coinciding with the onset of the $f_b$ -frequency component observed at $L_z = 3.8D$ .

Figure 22. Length of the recirculation region, $L_{\textit{recirc}}$ , and its gradient with respect to $L_z$ . The increase in $L_{\textit{recirc}}$ saturates with the emergence of the $f_b$ -frequency component in the range $3.7D \lt L_z \lt 3.8D$ .

Since the DMD mode at $f_k^{\textit{DMD}} = f_b$ is concentrated around the recirculation region, reconstructing the flow field using only this mode and the $0$ -frequency DMD mode may be more contentious. The reconstructed velocity field using selected DMD modes is given by

(4.2) \begin{align} \begin{split} [b_{k_1}(t_{\!j}), b_{k_2}(t_{\!j}), {\cdots}\, b_{k_s}(t_{\!j})]=\big[\boldsymbol{\varphi }^{\textit{DMD}}_{k_1}(\boldsymbol{x}),\boldsymbol{\varphi }^{\textit{DMD}}_{k_2}(\boldsymbol{x}), {\cdots}\, \boldsymbol{\varphi }^{\textit{DMD}}_{k_s}(\boldsymbol{x})\big]^{\dagger }\boldsymbol{u}(\boldsymbol{x},t_{\!j}),\\ \end{split} \end{align}

where $k_{\!j}$ denotes the indices of the selected DMD modes, and $b_{k_s}(t_{\!j})$ is the coupling coefficient corresponding to mode $\boldsymbol{\varphi }^{\textit{DMD}}_{k_s}$ at time $t = t_{\!j}$ . These coefficients are computed directly from the flow snapshot using

(4.3) \begin{align} \begin{split} [b_{k_1}(t_{\!j}), b_{k_2}(t_{\!j}), {\cdots}\, b_{k_s}(t_{\!j})]=[\boldsymbol{\varphi }^{\textit{DMD}}_{k_1}(\boldsymbol{x}),\boldsymbol{\varphi }^{\textit{DMD}}_{k_2}(\boldsymbol{x}), {\cdots}\, \boldsymbol{\varphi }^{\textit{DMD}}_{k_s}(\boldsymbol{x})]^{\dagger }\boldsymbol{u}(\boldsymbol{x},t_{\!j}).\\ \end{split} \end{align}

The Moore–Penrose pseudoinverse is computed using a QR decomposition for preconditioning. QR refers to the decomposition of a matrix into an orthogonal matrix Q and an upper triangular matrix R.

For the $L_z = 3.9D$ case, the coupling coefficients $b_k(t_{\!j})$ corresponding to the DMD modes with $f_k^{\textit{DMD}} = \pm f_b$ and $0$ were computed. The flow fields were then reconstructed using the DMD modes with $f_k^{\textit{DMD}} = \pm f_b$ and $0$ . The reconstructed flow fields were averaged in the spanwise direction. From these spanwise-averaged fields, the temporal variation of the recirculation length $L_{\textit{recirc}}$ was extracted. Figure 23 shows instantaneous spanwise-averaged fields corresponding to the maximum and minimum values of $L_{\textit{recirc}}$ . The difference in $L_{\textit{recirc}}$ between these two states was approximately $0.1D$ . Since the $0$ -frequency DMD mode contains no temporal fluctuations, the time-varying component of the reconstructed field originates solely from the DMD modes with $f_k^{\textit{DMD}} = \pm f_b$ . This indicates that the $f_b$ -frequency fluctuation is responsible for the expansion and compression of the recirculation region.

Figure 23. Spatial distribution of the spanwise-averaged reconstructed streamwise velocity $u'$ using the DMD modes at $f^{\textit{DMD}}_k = 0$ and $f_b$ , shown at the time of maximum (top) and minimum (bottom) recirculation region length. The bubble pumping phenomenon is characterised by the periodic expansion and contraction of the recirculation region.

The expansion and compression of the recirculation region at frequencies significantly lower than the vortex-shedding frequency have been reported in various body wakes (Berger et al. Reference Berger, Scholz and Schumm1990; Najjar & Balachandar Reference Najjar and Balachandar1998; Ohmichi et al. Reference Ohmichi, Kobayashi and Kanazaki2019; Yokota & Nonomura Reference Yokota and Nonomura2024; Yokota et al. Reference Yokota, Nagata and Nonomura2025), and this phenomenon is referred to as recirculation bubble pumping. Therefore, temporal variations in the recirculation length $L_{\textit{recirc}}$ can be interpreted as manifestations of recirculation bubble pumping.

It is also known that the low-frequency fluctuations originating from this mechanism appear in the spectrum of the drag coefficient (Najjar & Balachandar Reference Najjar and Balachandar1998), and the power spectral density results shown in figure 5 support this interpretation. However, in some complex flows (Ohmichi et al. Reference Ohmichi, Kobayashi and Kanazaki2019; Yokota & Nonomura Reference Yokota and Nonomura2024; Yokota et al. Reference Yokota, Nagata and Nonomura2025), recirculation bubble pumping has been shown to exist even when no distinct and sharp spectral peaks are observed in the low-frequency range. A similar situation arises in the present study at $L_z = 12D$ , where the drag coefficient spectrum exhibits no clear peak in the low-frequency range. Nevertheless, a weak local maximum is observed in the drag coefficient spectrum shown in figure 5(d), and the DMD mode shown in figure 9(d) reveals a structure that is symmetric about $y = 0$ , indicating that it influences the drag force. These observations suggest the presence of recirculation bubble pumping at $L_z=12D$ . Thus, constraining the spanwise domain size has enabled the clear identification of this low-frequency mechanism.

Returning to the relationship between $L_z$ and $L_{\textit{recirc}}$ shown in figure 22, the time-averaged length of the recirculation length $L_{\textit{recirc}}$ tends to saturate beyond $3.8D$ . This trend suggests that recirculation bubble pumping occurs around the $L_{\textit{recirc}}$ value at $L_z=3.9D$ . Since the most stable spanwise wavelength for Mode A is reported to be around $3.8D \lt L_z \lt 3.9D$ (Henderson Reference Henderson1997; Jiang & Cheng Reference Jiang and Cheng2017; Jiang et al. Reference Jiang, Cheng and An2017), the saturation of $L_{\textit{recirc}}$ observed for $L_z \gt 3.8D$ is likely governed by the intrinsic stability of the flow. From this perspective, the stabilisation of the spanwise wavelength leads to a corresponding stabilisation in the size of the recirculation region, around which the wake recirculation bubble pumping develops. Thus, increasing $L_z$ beyond the stability threshold enables the recirculation region to develop fully, providing the necessary conditions for the onset of recirculation bubble pumping.

4.3. Interaction between Kármán vortex-shedding mode and recirculation bubble pumping

We reconstructed the vortex structures at the frequencies $f_K$ , $f_K^+$ and $f_K^-$ using DMD modes corresponding to the $0$ -frequency component and the respective frequencies for the $L_z = 3.9D$ case. The vortices associated with $f^{\textit{DMD}}_k = f_K^+$ and $f_K^-$ are considered to arise from interactions between the DMD modes at frequencies $f_b$ and $f_K$ . In the presence of a fully developed Mode A, a distinct periodic structure emerges in the spanwise direction. Accordingly, the isosurfaces of streamwise vorticity $\omega _x$ characterise the wake vortex pattern, consistent with observations reported by Jiang & Cheng (Reference Jiang and Cheng2017), Jiang et al. (Reference Jiang, Cheng and An2017) and Jiang et al. (Reference Jiang, Cheng, Tong, Draper and An2016b ). The streamwise vorticity $\omega _x$ was computed from the reconstructed flow field as

(4.4) \begin{equation} \omega _x = \frac {D}{U_\infty }\! \left ( \frac {\partial w'}{\partial y} - \frac {\partial v'}{\partial z} \right )\!, \end{equation}

where $w'$ and $v'$ denote the reconstructed spanwise and transverse velocity components, respectively. Figure 24 presents instantaneous isosurfaces of $\omega _x$ for the three frequencies $f^{\textit{DMD}}_k = f_K$ , $f_K^+$ and $f_K^-$ . The vortices at these frequencies are spatially co-located, indicating that they are not independent structures but instead constitute a single vortex street formed through mutual interaction. Since the $f_K$ component corresponds to the most energetic component, the primary vortex structure is captured by the DMD mode at $f^{\textit{DMD}}_k = f_K$ .

Figure 24(d) presents overlay plots of isosurfaces corresponding to the three frequency components: $f_K$ , $f_K^+$ and $f_K^-$ . The vortices associated with $f_K^+$ and $f_K^-$ are distributed around the primary vortex at frequency $f_K$ , suggesting that spatial fluctuations of the dominant Kármán vortex give rise to vortices at slightly different frequencies. According to Jiang et al. (Reference Jiang, Cheng and An2017), increasing $L_z$ permits variations in the spanwise wavelength of Mode A, leading to vortex dislocations. Here, the term vortex dislocation, mainly investigated by Williamson (Reference Williamson1996), refers to the phenomenon in the cylinder wake where the spanwise wavelength undergoes a large variation. These vortex structures with different spanwise wavelengths exhibit slightly different shedding frequencies due to their distinct spatial distributions. Through nonlinear interactions among these frequency components, energy is transferred to the $f_b$ -frequency component, which appears as recirculation bubble pumping. Furthermore, the BMD interaction map shown in figure 17 indicates that these interactions occur near the recirculation region. Specifically, modes with different vortex-shedding frequencies interact in close proximity to the body, serving as an energy source for pumping within the recirculation region. Since the primary frequency components, such as $f_K$ , $f_t$ and $f_s$ , interact with the $f_b$ -frequency component, it is suggested that the recirculation bubble pumping, represented by the $f_b$ frequency, contributes to the maintenance of additional frequency components throughout the flow field. Thus, the emergence of recirculation bubble pumping induces spatial fluctuations in the dominant coherent structure, thereby promoting the appearance of multiple frequency components.

Figure 24. Isosurfaces of streamwise vorticity $\omega _x$ computed from flow fields reconstructed using DMD modes: (a) $f^{\textit{DMD}}_k = 0$ and $f_K$ ; (b) $f^{\textit{DMD}}_k = 0$ and $f_K + f_b$ ; (c) $f^{\textit{DMD}}_k = 0$ and $f_K - f_b$ ; and (d) overlay of all three vorticity isosurfaces. Vortices associated with $f^{\textit{DMD}}_k = f_K \pm f_b$ are distributed around the vortex corresponding to $f^{\textit{DMD}}_k = f_K$ . Due to interaction with the recirculation bubble pumping mechanism, fluctuations appear in the vortex structure at $f^{\textit{DMD}}_k = f_K$ .