Flow decomposition methods provide systematic ways to extract the flow modes, which can be regarded as the spatial distribution of a coherent structure. They have been successfully used in the study of wake, boundary layer, and mixing. However, real flow structures also possess complex temporal patterns that can hardly be captured using the spatial modes obtained in the decomposition. In order to analyze the temporal variation of coherent structures in a complex flow field, this paper studies the recurrence in phase space to identify the pattern and classify the evolution of the flow modes. The recurrence pattern depends on the time delay and initial condition. In some cases, the flow system will revisit a previous state regardless of the initial state, and in other cases, the system’s recurrence will depend on the initial state. These patterns are determined by the arrangement and interactions of coherent structures in the flow. The temporal order of the repetition pattern reflects the possible ways of flow evolution.

## I. INTRODUCTION

Flow fields of jets and wakes usually contain coherent structures that possess various spatial and temporal scales. The co-existence and interactions of these structures lead to patterns, such as the well-known Karman vortex street. In the past few decades, many techniques have been developed to identify coherent structures. For example, vortices can be captured using the Q-criterion, swirling strength, or λ_{2} method.^{1,2} These methods are widely used to study turbulence and engineering flows. However, there are other types of coherent structures than just vortices. To extract general flow structures from a complex flow without relying on any conditional criterion, the proper orthogonal decomposition (POD) is used in fluid dynamics.^{3,4} Briefly, the POD extracts flow modes based on spatial correlation. These modes are the orthogonal basis of a vector space. They can be regarded as the spatial distribution of coherent structures. The flow field can be decomposed using the modes, *ω*(** x**,

*t*) = Σ

_{i}

*a*

_{i}(

*t*)

*M*

_{i}(

**), where**

*x**ω*can be any distribution in the field

**,**

*x**a*is the decomposition coefficient that varies with time, and

*M*is the mode that depends only on

**. The field can be reconstructed using low order modes that contain a larger portion of the field energy; hence, the POD provides a framework for reduced order modeling.**

*x*^{3,5–9}The high order modes are usually considered to be incoherent. The POD has been used to study the structures in pipe flow, shear layer, jets, and boundary layers, among others. Nevertheless, the POD modes,

*M*

_{i}(

**), are stationary distributions containing no temporal information. To investigate spatial–temporal variations, Schmid proposes a dynamic-mode-decomposition (DMD) algorithm.**

*x*^{10}The dynamic modes approximate the eigenvectors of the Koopman operator.

^{11}Each mode is associated with a complex eigenvalue,

*ρe*

^{iθ}. This eigenvalue attributes a time-evolution factor to each mode. If

*ρ*> 1, the related mode becomes infinite as the time increases; if ρ < 1, the corresponding mode reduces to 0; the modes with ρ = 1 form a set of oscillators. The DMD is also widely used to extract flow structures,

^{11,12}and the usefulness of POD and DMD is sometimes compared in analyzing the same flow.

^{13}

Because the POD modes possess no direct temporal information and the DMD essentially regards a flow field as a linear system, neither of them is ideal to represent the temporal evolution of coherent structures in a complex flow, which is highly nonlinear. Therefore, researchers have developed other decomposition and analysis technologies, some of which include the temporal aspect of the flow modes. Bagheri^{14} analyzed the wake dynamics of a flow over a stationary cylinder and established a relationship between Koopman modes with the DMD algorithm. It is shown that the reduced model obtained by Koopman modes reflects the nonlinear global mode in Hopf bifurcation under a certain condition. Glaz *et al.*^{15} found a new normal form model based on Koopman mode decomposition to analyze the quasiperiodic intermittency phenomenon in the nonlinear field. Noack *et al.*^{16} combined the features of POD and DMD to analyze flow structures that evolved coherently in both space and time, which lead to modes that contain a single frequency and are spatially orthogonal to all other modes at all frequencies. Clainche *et al.*^{17} proposed the higher order dynamic mode decomposition and applied sequentially in time and space to obtain a spatiotemporal expansion of the flow field. This method can analyze complex signals and noisy data. Towne *et al.*^{18} developed the spectral POD, in which POD was done in the frequency domain of the flow. This method has been used to analyze the turbulent jet flow and shown better results than POD.

In fact, the temporal variation of modal coefficients can be more important than the mode, as the coefficients reflect the evolution of a system in the phase space. The purpose of this paper is to show that the recurrence of modal coefficients in the phase space can be used to identify and classify the temporal evolution of coherent flow structures. A wake flow downstream of two cylinders will be used as an example. The phase space is constructed using POD results. Due to the development of large data and computational power, the recurrence network method has recently been widely used to study the nonlinearity with high dimensionality,^{19,20} and the trajectories of the system in the phase space reveal global features of the flow. The concept of recurrence originates from the pioneering work of Poincaré.^{21} The state of a system is *X*(*t*) = [*a*_{1}(*t*), *a*_{2}(*t*), …, *a*_{n}(*t*)]^{T}, where the entries are measures of the system. In a finite state space, the system usually revisits the neighborhood of a previous state. In other words, there exist *ε* and *τ*, which satisfy ∥*X*(*t* + *τ*) − *X*(*t*)∥ < *ε*, where ||⋯|| is the norm of a vector. Recurrence is a fundamental property of a deterministic system. Intuitively, the re-visitation implies that the system finishes a loop of evolution and returns to the original state, within the tolerance of ε. In the following time, the system will start a new evolution. Therefore, by studying the way a system repeats itself, its dynamics can be analyzed and predicted to a certain extent. The recurrence network has been applied to study the phase space geometry, characterize the temporal patterns, and identify sudden changes in a dynamical system.^{20,22–24} In the early development of recurrence network methods, the state *X*(*t*) is usually constructed using retarded single-point measurements.^{25–27} Recently, multiplex recurrence networks are developed to analyze multi-variant time series,^{28} but these variables are not necessarily related to any coherent structure. In this study, we use the POD modes to construct the state space, i.e., the entries of the state vector *X*(*t*) are composed of POD coefficients. The state vector reflects the temporal evolution, and it is also related to specific flow modes. Hence, the so-constructed network carries both spatial and temporal information of a complex flow field.

In the following, the numerical simulation of the wake flow downstream of two cylinders will be presented, construction of the recurrence map will be described, and the repetition patterns and their relationship with each POD mode will be discussed.

## II. NUMERICAL SIMULATION AND NETWORK CONSTRUCTION

### A. Computation setup

The flow studied in this paper is the wake downstream of two cylinders with the same diameter (Fig. 1). The upstream cylinder is fixed, and the downstream cylinder oscillates transversely in the flow. The immersed-boundary method is employed to simulate the flow, as it is effective in dealing with moving structures and no body-conforming mesh is required around objects in motion. Specifically, the improved direct-forcing immersed-boundary method^{29} is used, in which the solid cylinders are considered as a porous medium^{30} with a large resistivity. The model equations are the Navier–Stokes equations for incompressible flow, with the modified Zwikker–Kosten (Z–K) equation^{31} for flow inside the cylinder body. All the equations and variables are nondimensionalized with the free stream velocity and the diameter of the cylinder. The modified governing equations for incompressible, unsteady, viscous fluid flow are written as

where ** f** is the body-force term representing the virtual boundary force, and

*Re*is the Reynolds number defined as

*ρU*

_{∞}

*D*/

*μ*, where

*U*

_{∞}is the free stream velocity and

*D*the diameter of the cylinder. In order to restrict the flow to two dimensional, the simulation is performed at a Reynolds number of 100, which is the same as that in Zhang and Zheng.

^{32}The flow, both outside and inside the cylinders, can be simulated by the same format of the above governing equations, with the definition of the forcing term as

where *σ* is the dimensionless flow resistivity of the cylinders in the flow and *V* is the moving velocity of the objects. For the upstream cylinder, *V* is zero; for the downstream cylinder, *V* is the oscillating velocity of the cylinder. Thus, the original Navier–Stokes equations are used for fluid flow region and the modified equations, with a source term in Eq. (3) applied, are used inside the cylinders. The equations are solved on a nonmoving staggered Cartesian grid.

The simulation domain and the configuration of the problem are shown in Fig. 1, where a pair of cylinders in tandem arrangement is placed in a uniform, incompressible, viscous flow. The computational domain size for the two-dimensional study is [0, 38.4] × [0, 25.6], where the length and width of the domain are nondimensionalized by the cylinder diameter, *D*. The center of the upstream cylinder is located 8*D* from the inlet flow boundary to reduce the inlet effect and 30.4*D* away from the outlet boundary to ensure unstrained vortex wake development. The center-to-center distance between the cylinders is 6, by which the flow is in the vortex formation regime.^{33} The size of the uniform, nonmoving grid used in this study is $\Delta x=\Delta y=0.025$ for the computational case. The check for grid-independent solution has been carried out in our previous work.^{29} The Dirichlet-type boundary condition is employed at the inlet for velocity; the symmetry boundary condition is used for both upper and lower boundaries; and the outlet is specified with the Neumann-type boundary condition. In addition, the computational scheme uses a low-storage third-order Runge-Kutta scheme for time,^{34} the fifth order WENO-Z method for convection,^{35} and the second-order central difference for viscous terms. The incompressibility condition is satisfied by solving a Poisson equation for pressure correction using MUDPACK.^{36} For the transversely oscillating downstream cylinder, the sinusoidal motion is specified as

where *h* = 0.15 is the heaving magnitude, *f*_{c} = 0.211 is the oscillating frequency, and *t* is time. The computational scheme has been verified and validated with numerous computational and experimental data.^{37,38}

### B. POD for the vorticity field

A 3.75 × 12.5 domain is selected for POD and recurrence network analysis (see Fig. 1). This location is chosen because the vortex shedding pattern changes around 11*D* away from the downstream cylinder center.^{33} In this study, vorticity is used to analyze the flow dynamics. Using POD, the vorticity field can be decomposed into orthogonal modes with the corresponding time coefficients.^{39} The energy contained in the modes is measured using the sum of the eigenvalues, and the energy percentage captured by the ** i**th POD mode is given by $\lambda i/\u2211k=1Ne\lambda k$, which is also referred to as the energy of each mode, where

*N*_{e}is the total number of eigenvalues. The cumulative energy up to the

**th POD modes can be represented by $\u2211k=1i\lambda k/\u2211k=1Ne\lambda k$.**

*i*^{40}

After the simulation reaches 100 oscillating cycles, the vorticity results of 70 oscillation cycles (from 101*T* to 170*T*, where *T* = 1/*f*_{c}) are chosen for POD and recurrence map analysis. There are 2048 time steps in each oscillating cycle. To ensure the convergence of the POD modes, sampling at two different time intervals (*T*/64 and *T*/128 between two successive snapshots) is tested. The *L*^{2} norm difference for the first twelve modes is calculated and the detail of the calculation is explained in Ref. 40. The results show that the maximum difference is below 0.5% for time interval equal to *T*/64. Thus, the flow information in the selected domain with time interval *T*/64 is used. In this way, we totally obtained 4480 time series data.

The first twelve eigenvalues, corresponding to the first twelve energetic modes, are listed in Table I with their percentages of energy contribution and cumulative energy. Clearly, the first two modes dominate and comprise more than 55% of the total energy of the motion. The first twelve modes collectively contain more than 95.5% of the energy. Thus, within these modes we have almost captured the entire spatial structure of the flow field. The vorticity contours of modes 1 through 6 and the time histories of the first two modal coefficients are shown in Fig. 2. The first four modes show symmetric structures with respect to the wake centerline, while the fifth and sixth modes are anti-symmetric. In addition, the coefficients of the first two modes illustrate no exact periodicity.

Modes . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . | 9 . | 10 . | 11 . | 12 . |
---|---|---|---|---|---|---|---|---|---|---|---|---|

Energy of each mode (%) | 32.95 | 23.93 | 7.64 | 7.50 | 6.90 | 5.90 | 4.41 | 1.86 | 1.79 | 1.24 | 0.88 | 0.85 |

Cumulative energy (%) | 32.95 | 56.88 | 64.52 | 72.02 | 78.92 | 84.82 | 89.23 | 91.09 | 92.88 | 94.12 | 95 | 95.85 |

Modes . | 1 . | 2 . | 3 . | 4 . | 5 . | 6 . | 7 . | 8 . | 9 . | 10 . | 11 . | 12 . |
---|---|---|---|---|---|---|---|---|---|---|---|---|

Energy of each mode (%) | 32.95 | 23.93 | 7.64 | 7.50 | 6.90 | 5.90 | 4.41 | 1.86 | 1.79 | 1.24 | 0.88 | 0.85 |

Cumulative energy (%) | 32.95 | 56.88 | 64.52 | 72.02 | 78.92 | 84.82 | 89.23 | 91.09 | 92.88 | 94.12 | 95 | 95.85 |

### C. Recurrence network construction

As discussed before, the phase space is constructed using the POD coefficients *X*(*t*) = [*a*_{1}(*t*), *a*_{2}(*t*), …, *a*_{n}(*t*)]^{T}, where *a*_{1}–*a*_{n} are the modal coefficients. For this study, n = 12, i.e., the first twelve modes that contain more than 95% energy. For convenience, we refer to the state at time *t*_{i} as the *i*th state *X*_{i}. The recurrence in the phase space can be represented using the recurrence network, which is obtained in the following way. Given the states of a system *X*_{i} and *X*_{j} at the time *t*_{i} and *t*_{j}, two nodes are used to represent these states. If the distance between two states is smaller than a threshold *ε*, recurrence happens and the corresponding nodes in the network can be connected. In other words, the network’s adjacency matrix is *R*_{ij} = *H*(*ε* − ∥*X*_{i} − *X*_{j}∥), where *H* is the Heaviside function and ||⋯|| is the norm of a vector. The standard *L*^{2} norm is used for this research though other norms yield similar results. The *R* matrix is symmetric, *R*_{ij} = *R*_{ji}. This matrix can be conveniently plotted using binary images, i.e., regarding *R*_{ij} as a pixel (*i*, *j*) on an image. Figure 3 shows the recurrence plot using ε = 60. The determination of the value of ε will be discussed in the following paragraphs. The black and white pixels are inverted for clarity, i.e., 1 – *R* is shown in the figure. Hence, a black pixel (*i*, *j*) indicates that the nodes *i* and *j* are connected. Note that in the recurrence plot, the vertical axis *j* is pointing upward, which is different from the normal arrangement of pixels on an image. Since a state is always similar to itself, *R*_{ii} = 1. Consequently, the secondary diagonal line must be black. In the recurrence plot, both indices *i* and *j* indicate a time direction. The lines parallel to the secondary diagonal represent the recurrence after a certain time interval *i* − *j*. When studying the temporal variation, we can fix one index and vary the other. These line patterns will be discussed in detail in Sec. III.

The primary issue in the network construction is to determine the parameter ε. Figures 4(a) and 4(b) show the recurrence matrices obtained using ε = 10 and 100. The ε threshold defines the similarity threshold between two states. Since a chaotic system never exactly repeats itself, a tiny ε excludes the recurrences that actually characterize the dynamic system. As a result, fewer recurrences can be recognized [Fig. 4(a)]. On the other hand, a large ε allows too many states to be regarded as recurrent states. Figure 4(e) shows the magnitude of ∥*X*_{i} − *X*_{j}∥ using a contour line. The strokes become fatter as ε increases. At ε = 100, the line patterns are so swollen that they cease to exist. The black regions [Fig. 4(b)] extend so much that they get connected to each other, i.e., they percolate through the entire recurrence plot.

To select a proper ε, we need to employ the recurrence quantitative analysis. As illustrated in the inset of Fig. 3, we first define a connected line structure or a stroke as the set *C* = {(*i*, *j*) | *R*_{ij} = 1 and all path-connected elements are included}. In other words, the stroke is a path-connected black region on the recurrence plot, and there are many such strokes on the plot. The stroke length is defined as *l* = 1 2[max (*i* + *j*) − min (*i* + *j*)], for any (*i*, *j*) ∈ *C*. The collection of all structures allows us to calculate the mean structure length ⟨*l*⟩_{ε} and the probability *p*(*l*). Consequently, the Shannon entropy of this distribution can be computed as

The Shannon entropy is a proper way to examine the selection of ε because it considers both the total number of structures on the recurrence plot and their length distribution. If ε is too small and there are few structures, the Shannon entropy is small; on the other hand, if ε is too big, the probability concentrates at strokes of a large *l*, and the entropy is also small. The property of the Shannon entropy and the related probability distribution are extensively studied in statistical physics.^{41–44} Figure 4(c) shows the variation of ⟨*l*⟩_{ε} vs ε. The ⟨*l*⟩_{ε} remains at a low level and suddenly increases at ε > 70. The *E*_{ε} − ε curve also has a significant change at ε > 70 [Fig. 4(d)]. The transition at ε ∼ 70–90 reflects the percolation of swollen line strokes on the recurrence plot. In addition, the value of entropy is also small at ε < 40, as there are fewer strokes on the recurrence map and the probability *p*(*l*) concentrates at a few *l*’s. As ε increases, more structures can be identified and the entropy increases. Based on the above observation, a proper ε should be chosen so that it is slightly smaller than the percolation transition. In this paper, ε = 60 is used.

## III. RESULTS AND DISCUSSION

For the convenience of analysis, we use the following transform to rotate and stretch the recurrence plot: η = (*i* + *j*) and ξ = (*j* − *i*). Clearly, ξ is the time delay between *i* and *j* states, and ξ = 0 corresponds to the secondary diagonal of the recurrence plot (Fig. 3). The resulting image is in Fig. 5. Since the adjacency matrix is symmetric, only ξ ≥ 0 part is plotted. The plot is largely white with many horizontal lines. These recurrence structures form a layered pattern that looks like an Egyptian pyramid. The thickness of each layer is nearly a constant. Note that there are two major time scales in the flow: the oscillation period *T* and the vortex shedding period *T*′. In our flow, *T*′ = 1.3*T*. For a flow system to repeat a state, the boundary condition should be the same. After a period of *mT*, where *m* is an integer, the oscillating cylinder returns to the same position with the same velocity. Therefore, the boundary condition restores and this facilitates the repetition of vorticity distribution. However, this time may not be an integer time of *T*′. As a result, the inception of recurrence depends on the initial flow state *X*_{j}. This is manifested by the discontinuous horizontal lines in the recurrence plot. If the time delay is an integer time of both *T* and *T*′ (e.g., *t* = 13*T* = 10*T*′), the system has a much higher chance to repeat.

Figure 5 reveals four types of horizontal line patterns: solid line (pattern A), dashed line with longer stroke (pattern B), dashed line with shorter stroke and dots (pattern C), and staggered dashed-dotted line (pattern D). Note that the dashes and dots in pattern C do not have the same ξ value, and pattern D can be considered as a composition of two dished-line structures. In the ξ direction, these patterns seem to appear sequentially. In the following, we will analyze these patterns and their appearance order.

### A. Horizontal line patterns in the $\eta $–$\xi $ frame

Our analysis begins with the simplest pattern—-the solid lines, which are located at

Note that 4*T* ≈ 3*T*′, 13*T* = 10*T*′, 17*T* ≈ 13*T*′…. All these numbers are also close to integer times of *T*′, where the vortex shedding pattern should repeat itself. In addition, the repetitions appear in pairs (0 and 4*T*, 13*T* and 17*T*, 26*T* and 30*T*, …). The interval between each pair is always 4*T*, i.e., another solid line is located after a 4*T* period following the first recurrence. A solid line means that the system always revisits a previous state after time ξ, i.e., ∥*X*_{i} − *X*_{j}∥ ≤ *ε* for a fixed ξ = *I* − *j*, where *j* is arbitrary. To show how the recurrence happens, a component-wise analysis is carried out to study each mode. We choose the ξ = 835 line as an example and examine the difference between each individual mode. The coefficient difference of the *k*th mode is *S*_{jk} = *a*_{k}(*t*_{i}) − *a*_{k}(*t*_{j}) (ξ is fixed at 835). The coefficient differences of all modes are illustrated in Fig. 6(a), where the modes are listed along the vertical axis and the time index *j* is arrayed along the horizontal direction. The color represents the value of *S*_{jk}. The extreme values of *S*_{jk} of the first 9 modes are on the order of ±10, and their magnitudes are much larger than those of the other 3 modes. The temporal variation of the norm ∥*X*_{i} − *X*_{j}∥ is also shown in Fig. 6(a). It fluctuates around 10 and stays far below the ε = 60 threshold. Along the *j* axis, the extreme values of *S*_{jk} appear at a regular pace. Note that the first two modes contain more than 55% cumulative energy (Table I), and the extreme values of the coefficients *a*_{1} and *a*_{2} are on the order of ±60. Therefore, in terms of the extremes, the coefficient difference *S*_{j1} and *S*_{j2} is an order of magnitude smaller than *a*_{1} and *a*_{2}, respectively. In other words, the energetic modes repeat themselves very well.

As a comparison, the modal coefficients of non-recurrence cases are studied, a horizontal line at ξ = 900 is randomly chosen. The coefficient differences are shown in Fig. 6(b). The magnitudes of *S*_{j1} and *S*_{j2} are on the order of 100, which are nearly twice of *a*_{1} (or *a*_{2}) and are much larger than the other *S*_{jk} (for *k* > 2). This is in contrast to the *S*_{jk} behavior in the above pattern A. Examination of other non-recurrent ξ values shows a similar result. Therefore, recurrence only happens at these ξ values where the low order energetic modes repeat themselves. As discussed before, the oscillating cylinder moves back to a previous position only at certain ξ values, creating a repeated boundary condition. At other ξ values, the boundary conditions are different, the large-scale flow profile is not the same, and no repetition occurs. Hence, recurrence only occurs at a set of discrete ξ values.

Pattern B (dashed line with long strokes) indicates that, at a certain ξ, the recurrence is switched on and off, depending on the starting time. We choose ξ = 1320 as an example to study the modal behaviors along pattern B. Figure 7(a) shows the magnitudes of *S*_{jk}. Generally speaking, the magnitude of *S*_{jk} with a small *k* is larger than that with a big *k*. The extreme values of *S*_{jk} for the first 6 modes are on the order of ±50, which is close to the extreme values of coefficients *a*_{1} and *a*_{2}. The norm ∥*X*_{i} − *X*_{j}∥ fluctuates around the ε threshold [Fig. 7(a)], so it is not surprising to see the on-and-off pattern. In order to find out which mode causes the on-and-off of recurrence, we define the significance of each mode as a product of weight and correlation, *W*_{k} = *P*_{k} ⋅ *Q*_{k}. The weight for the *k*th mode is measured using an average percentage $Pk=\Sigma j=1N(Sjk2/\Sigma kSjk2)/N$, where *N* is the total number of index *j* and $\Sigma kSjk2$ is ∥*X*_{j+ξ} − *X*_{j}∥^{2}. *P*_{k} measures the contribution of a mode to the norm. The correlation *Q*_{k} of mode *k* characterizes the synchronization between the mode’s variation, *S*_{jk} − ⟨*S*_{jk}⟩, and that of the overall norm, ∥*X*_{j+ξ} − *X*_{j}∥ − ⟨∥*X*_{j+ξ}−*X*_{j}∥⟩,

where *j* = 1, …, *N* and the correlation is normalized. A large positive *W*_{k} means that the *k*th mode significantly contributes to ||*X*_{j+ξ} − *X*_{j}|| and is positively correlated with *S*_{jk}. Figure 7(b) shows the significance factor *W*_{k} of each mode. Mode 1, containing the largest energy, is the most significant. It is interesting that the second important mode is actually mode 5, not mode 2 or 3 though their modal energy is larger (Table I).

Pattern B shows that a higher order mode can be more important than a lower order mode in terms of recurrence dynamics. In fact, modes 5 and 6 are crucial because they are the first pair of anti-symmetric modes, as shown in Fig. 2(a). Note that the vorticity field of a cylinder wake (e.g., the Karman vortex street) is nearly anti-symmetric after shifting along the wake axis. Hence, the coefficients of modes 5 and 6 are closely related to the stream-wise position of vortices. For example, if the modal coefficient *a*_{5} is large, a strong vortex is located at the center of the selected window; if *a*_{6} is large, the vortices are near the edge of the window. The repetition of the field requires the vortex to be at the same location; hence, the coefficients of anti-symmetric modes should be the same. If the vortex shedding at time *i* and *j* are not in-phase, the recurrence is broken. Therefore, pattern B becomes a partial recurrence pattern. These observations show that anti-symmetric modes are crucial to the description of wake flow structures; even these modes are of high-order and possess less energy than the low-order symmetric modes.

The above modal analysis method is used to analyze pattern C (dashed-dotted line). Here, the pattern near ξ = 1176 is selected as an example. For the dashes, the temporal variation of *S*_{jk} looks similar to that of pattern B, but the modal significance has a different order: the 5th mode is significant but less important than mode 2 [Fig. 8(a)]. Regarding the small dots between the dashed lines, they have a ξ = 1222 that slightly differs from the dashed lines. The mode significance analysis along the dots indicates that the first two modes are crucial to the appearance of the dots in pattern C, but mode 5 has no significance [Fig. 8(b)]. Therefore, the comparison between patterns B and C (including both dashes and dots) shows that the first mode is always the most important mode for partial recurrence, and the importance of higher-order modes is related to the length of the stroke in the recurrence plot: a higher importance is related to a longer stroke. For the dots, the *W*s of high-order modes (especially for *k* ≥ 5) are tiny.

Pattern D is composed of two rows of strokes that are shorter than dashes but longer than dots in pattern C. According to the previous observation about the modal significance, we anticipate that the first two modes are the dominant modes and modes 3 and 4 have small *W*_{k}. In addition, the significance of modes 5 and 6 should be between those of dashes and dots in pattern C. The results in Fig. 9 confirm the expectations. Based on the above results, we here, discuss the relationship between the mode energy and its significance factor in recurrence. In general, low-order modes are more important since they possess higher energy. However, their significance factor is largely independent of the energy in the sense that, if the energetic mode repeats itself very well, it contributes nothing to recurrence breakdown. In the partial recurrence patterns, the low-order modes become less significant because they are in phase, hence the high-order modes determine the inception of recurrence, i.e., they are significant.

### B. Flow evolution and patterns in the vertical direction

The above horizontal patterns reflect the behavior of all modes as a collection. In this section, we will examine the appearing order of the patterns in the vertical ξ direction, which reflects the way all modes evolve in time. Pattern A indicates definite recurrence, regardless of the initial time. Hence, we partition the vertical order of all patterns according to the appearance of pattern A. In our case of low-Re wake behind two tandem cylinders, the first partition (from ξ = 0 to 245) is ADA, followed by ABBBCA, ADA, ACBBCA, ADA, …, as illustrated in Fig. 10, which is the center portion of Fig. 5. The fact that some letter sequences (e.g., ADA) keep reappearing suggests that these patterns have temporal coherence.

These repeated pattern sequences in fact reflect the way-of-evolution of the system. Taking the ADA partition from ξ = 835–1080 as an example, we can select any initial state *X*_{j}. After *γ* = 835 steps, a definite repetition occurs, ||*X*_{j+γ} − *X*_{j}|| < *ε*; after *δ* = 1080 steps, another definite repetition occurs ||*X*_{j+δ} − *X*_{j}|| < *ε*. In the phase space, the loop from *X*_{j+γ} to *X*_{j+δ} forms a quasi-closed loop with a gap $||Xj+\delta \u2212Xj+\gamma ||\u2264||Xj+\delta \u2212Xj||+||Xj+\gamma \u2212Xj||\u22642\epsilon $. For some initial state *X*_{j}, this loop has a twist as the system revisits the initial state between *X*_{i+γ} and *X*_{i+δ}, which is dubbed partial repetition in the previous discussion. The combination of partial and definite repetition can be used to classify the system’s temporal evolution. In our simulation, the most commonly appeared evolution is ADA and ACBBCA.

## IV. SUMMARY

In this paper, we examine the POD modes and the recurrence of flow dynamics in the phase space to analyze the temporal features of coherent structures. The POD modes are spatial distributions of vorticity, which do not have time variation. Although the wake flow in our study is chaotic, as discussed in Ref. 33, the recurrence results show that the system keeps revisiting a previous state after a longer period of time. Here, the recurrence is not exact, and the tolerance is selected to be slightly smaller than the percolation threshold. For the studied low-Re wake flow, its recurrence is determined by two factors. The first is the oscillation period of the downstream cylinder, and the second is the vortex-shedding period. The system has a high chance to revisit a previous state once the flow boundary condition is similar and the vortex shedding is in accordance. The first factor indicates that the recurrence plot should occur at nearly integer times of the oscillation period (and perhaps the vortex shedding period), which is shown as the patterns parallel to the secondary diagonal in the recurrence plot. The second factor results in partial recurrence, which is reflected by discontinuous lines in the plot. If the flow boundary condition is distinctive, the large-scale vorticity profiles are considerably different. As a result, the first few POD modes, which are the most energetic and describe the large-scale field profile, cannot repeat themselves and recurrence is generally impossible. However, the modal energy is not the only critical parameter that defines the role of a mode in the recurrence. If the energetic modes repeat themselves very well after a certain period of time, the recurrence may still break down due to high-order modes. As in the partial recurrence cases, higher-order modes (especially the first anti-symmetric mode 5) play a key role. The organization of full and partial recurrence contains a few sequences of patterns, which have a fixed order such as ADA. They reveal the coherence in the time direction. The analysis methodology presented in this study can be applied to investigate the quasi-periodic flow structures that appear intermittently in a turbulent flow.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

The authors would like to thank the support from the University of Kansas GO program and General Research Funding. H.W. is also grateful to Justin Phommachanh for his help with the manuscript preparation.