Proper orthogonal decomposition and recurrence map for the identification of spatial–temporal patterns in a low-Re wake downstream of two cylinders

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. © 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/1.5144978., s


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) = Σiai(t)Mi(x), where ω 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 x. 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. 3,[5][6][7][8][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, Mi(x), are stationary distributions containing no temporal information. To investigate spatial-temporal variations, Schmid proposes a dynamic-mode-decomposition (DMD) algorithm. 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), . . ., an(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][23][24] In the early development of recurrence network methods, the state X(t) is usually constructed using retarded single-point measurements. [25][26][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.

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 bodyconforming 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 ∂u ∂t 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. The center of the upstream cylinder is located 8D from the inlet flow boundary to reduce the inlet effect and 30.4D away from the outlet boundary to ensure unstrained vortex wake development. The center-tocenter 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 Δx = Δ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 11D 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 ith POD mode is given by λi/ ∑ N e k=1 λ k , which is also referred to as the energy of each mode, where Ne is the total number of eigenvalues. The cumulative energy up to the ith POD modes can be represented by 40 After the simulation reaches 100 oscillating cycles, the vorticity results of 70 oscillation cycles (from 101T to 170T, 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.

C. Recurrence network construction
As discussed before, the phase space is constructed using the POD coefficients X(t) = [a 1 (t), a 2 (t), . . ., an(t)] T , where a 1 -an 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 ti as the ith state Xi. 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 Xi and Xj at the time ti and tj, 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 Rij = H(ε − ∥Xi − Xj∥), 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, Rij = Rji. This matrix can be conveniently plotted using binary images, i.e., regarding Rij 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, Rii = 1.

FIG. 3.
Recurrence plot with ε = 60. The black and white pixel is inverted, i.e., 1 -R is shown. A black stroke is a connected region. Its length is l.

ARTICLE scitation.org/journal/adv
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 ∥Xi − Xj∥ 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) | Rij = 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 Eε = −Σp(l)log 2 (p(l)).
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

ARTICLE
scitation.org/journal/adv studied in statistical physics. [41][42][43][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.3T. 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 Xj. 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 = 13T = 10T ′ ), 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 dasheddotted 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 η-ξ frame
Our analysis begins with the simplest pattern--the solid lines, which are located at ξ = {0, 245, 835, 1080, 1675, 1918, 2506, 2747, 3344, 3587 Note that 4T ≈ 3T ′ , 13T = 10T ′ , 17T ≈ 13T ′ . . .. 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 4T, 13T and 17T, 26T and 30T, . . .). The interval between each pair is always 4T, i.e., another solid line is located after a 4T period following the first recurrence. A solid line means that the system always revisits a previous state after time ξ, i.e., ∥Xi − Xj∥ ≤ ε 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 kth mode is S jk = a k (ti) − a k (tj) (ξ 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 ∥Xi − Xj∥ 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 nonrecurrent ξ 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  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 ∥Xi − Xj∥ 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 kth mode is measured using an average percentage P k = Σ N j=1 (S 2 jk /Σ k S 2 jk )/N, where N is the total number of index j and Σ k S 2 jk is ∥Xj+ ξ − Xj∥ 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, ∥Xj+ ξ − Xj∥ − ⟨∥Xj+ ξ −Xj∥⟩, where j = 1, . . ., N and the correlation is normalized. A large positive W k means that the kth mode significantly contributes to ||Xj+ ξ − Xj|| 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 ARTICLE scitation.org/journal/adv 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 a5 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 Ws 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-ofevolution of the system. Taking the ADA partition from ξ = 835 -1080 as an example, we can select any initial state Xj. After γ = 835 steps, a definite repetition occurs, ||Xj+γ − Xj|| < ε; after δ = 1080 steps, another definite repetition occurs ||Xj+ δ − Xj|| < ε. In the phase space, the loop from Xj+γ to Xj+ δ forms a quasi-closed loop with a gap ||X j+δ − Xj+γ|| ≤ ||X j+δ − Xj|| + ||Xj+γ − Xj|| ≤ 2ε. For some initial state Xj, this loop has a twist as the system revisits the initial state between Xi+γ and Xi+ δ , 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.