Twin vortex computer in fluid flow

Fluids exist universally in nature and technology. Among the many types of fluid flows is the well-known vortex shedding, which takes place when a fluid flows past a bluff body. Diverse types of vortices can be found in this flow as the Reynolds number increases. In this study, we reveal that these vortices can be employed for conducting certain types of computation. The results from computational fluid dynamics simulations showed that optimal computational performance is achieved near the critical Reynolds number, where the flow exhibits a twin vortex before the onset of the Kármán vortex shedding associated with the Hopf bifurcation. It is revealed that as the Reynolds number increases toward the bifurcation point, the input sensitivity of the twin vortex motion also increases, suggesting the modality of information processing within the system. Our finding paves a novel path to understand the relationship between fluid dynamics and its computational capability.


Introduction
Fluids can be found everywhere in nature, exhibiting rich and diverse dynamics, as well as instabilities, which are a source of inspiration for many mathematicians, physicists, engineers, and biologists in the field of nonlinear science. For example, vortices in fluid dynamics play an active role in biological behavior control, such as the swimming of fish [1] or flapping of butterflies [2]. In this study, we focus on the diverse nature of fluid dynamics and make a novel attempt to elucidate its potential as a computational resource.
A vortex shedding is a renowned phenomenon in fluid dynamics that takes place when a flow collides with an obstacle in the direction of fluid travel (figure 1); a typical example is the flow past a circular cylinder. The mechanism generating vortices has attracted the attention of many researchers [3][4][5][6]. The flow in this system is steady and symmetrical at low Reynolds numbers, but at higher Reynolds numbers, a pair of vortices, known as a twin vortex, are generated behind the object. A twin vortex grows in proportion with the Reynolds number; however, when the Reynolds number exceeds a certain value, vibration occurs downstream, resulting in a Kármán vortex, which has two alternating rows of vortices (cf figure 2(A) and supplementary videos (https://stacks.iop.org/NJP/23/063051/mmedia)). In this study, the flow past a circular cylinder was used as a test bed for clarifying the computational performance of fluid dynamics and the types of vortices they form.
Several attempts have been made to implement computations that focus on fluids. The construction of logic circuits using the fluidity of a droplet and the potential applications for soft robots are typical examples [7][8][9]. In this paper, we propose another paradigm of information processing based on fluids, where we managed a machine learning framework called reservoir computing (RC) [10,11]. RC is a framework of recurrent neural network learning and an information processing technique that exploits the input-driven transient behaviors of high-dimensional dynamical systems, called a reservoir. This method can be implemented by adjusting linear and static readout weights from a reservoir. Moreover, an arbitrary Behavior of a typical vortex shedding when the flow is past a circular object, and an illustration of node preparation for the reservoir. The grid behind the cylinder represents the nodes used for the reservoir; let the set of nodes be Ω R . The time series plots are normalized results (u is velocity, and p represent pressure) of the fluid simulation and the input. However, rather than using all such time series in the reservoir, the time data immediately preceding the next input are extracted and used. The total number of output/observation nodes is (10 + 1) × (6 + 1) × 3 = 231. Only several time series from a single spatial location are shown for clarity. All these diagrams are drawn under the conditions of Re = 50 and r k ∈ [0.99, 1.01]. , which is simply denoted by u 2 in this caption. The mean and deviation for max P∈Ω R ,k∈[0,Tr−1] u 2 (P, (k + 1)τ ) and min P∈Ω R ,k∈[0,Tr−1] u 2 (P, (k + 1)τ ) are plotted according to the Reynolds number. (C) The phase portrait (u 2 , ∂u 2 /∂t), which was constructed at (x 1 , x 2 ) = (5.5, 0.0) and Re = 10, 40, 50, and 100. time series can be accurately generated when the reservoir exhibits sufficient nonlinearity and memory [12]. A recent development in RC is the use of the dynamics of a physical system as a reservoir, which can be referred to as physical reservoir computing (PRC) [13]. Examples of this implementation can be found in a range of physics fields, from mechanical [14][15][16] to photonics [17][18][19], spins [20][21][22][23][24][25], and quantum [26][27][28], where each platform exhibits a particular computational property intrinsic to its respective spatio-temporal scales. It is important to note that the first PRC was implemented using the dynamics of a water surface [29], where the mechanical inputs are injected to the surface as perturbations and the resulting shock waves propagate and interact each other to process input information. We build on this seminal approach and investigate the computational account of vortex types in fluid flow, which has not been clarified to date. Similar studies have recently been proposed to implement a recurrent neural network in physical wave systems [30] and to exploit nonlinear waves, such as rogue waves, dispersive shocks, and solitons, as a reservoir for neuromorphic devices [31,32]. Because vortices can be found universally in nature, if we could exploit them for our computational purposes, then we can outsource the computational loads to physical dynamics, which results in energy-efficient information processing.
In this study, we exploited the dynamics of flow past a circular cylinder as a reservoir to solve temporal machine learning tasks, and thereby systematically analyzed the fluid computational capabilities. The reservoir dynamics was simulated through Navier-Stokes equations expressed by a set of partial differential equations (PDE). First, we investigated the bifurcation structure of the reservoir driven by the input stream to understand its dynamic property according to the Reynolds number and show that it exhibits Hopf bifurcation. Second, we analyzed the echo state property (ESP) [33], which provides a prerequisite for dynamics to function as a successful reservoir. We will demonstrate the relationship between the ESP and the stability of the fluid dynamics associated with the bifurcation structure. It was found that the ESP breaks down when the Reynolds number reaches the critical bifurcation point at the onset of the Kármán vortex shedding associated with symmetrical breaking of the twin vortex. Third, we focus on the morphological behavior of the twin vortex before the critical point and reveal that the long-diameter oscillation of the twin vortex is sensitive to input stream, suggesting a mechanism for the processing of input information in the system. Next, we implemented two typical benchmark tasks in the machine learning context, using the reservoir to demonstrate the range of computational capabilities and the dependence of the computational performance on the Reynolds number. The optimal computational performance was achieved near the critical Reynolds number, where the flow exhibits a twin vortex before the onset of Kármán vortex shedding. Finally, we demonstrate a prediction task of missing fluid variables, using PRC as an illustration of the potential application scenario of our scheme. To the best of our knowledge, this study is the first to discuss the computational capability of vortices generated using fluid dynamics.

Results
We start by introducing our reservoir settings in the next subsection. Subsequently, we move to analyze the reservoir property from dynamical systems point of view, including bifurcation structures and input-driven phenomena. Then, we investigate the information processing capability of our reservoir through a number of simulation experiments and show how the dynamical property of the system affects the information processing capability.

Flow past a circular cylinder as a reservoir
As for our reservoir system, we considered flows past a circular cylinder governed by the two-dimensional Navier-Stokes equations with inputs and solved it numerically using the stabilized Lagrange-Galerkin (LG) method [34,35]. The details of this problem, this method, and the input settings are provided in appendix and in the supplementary materials (SM). Here, we provide the definition of our reservoir using the numerical solution (u h , p h ) = ((u h1 , u h2 ) T , p h ), which approximates the pair of velocity u = (u 1 , u 2 ) T and pressure p of fluid flow, where h is a mesh size in the computation and the superscript T means the transposition. We used two parameters, Δt and τ , for time in the flow simulation. Δt is a time-step size in the simulation and τ is a transient time for the reservoir that determines its timescale. Although the sizes of h and Δt are related to the 'error' between the numerical solution (u h , p h ) and the exact solution (u, p), we describe the setup of our reservoir here without Δt to avoid confusion and with subscript h to indicate only that (u h , p h ) is numerical solution, and the set of values of (u h , p h ) at point P = (P 1 , P 2 ) T ∈ R 2 and time kτ (k ∈ N) is denoted by (u h , p h )(P, kτ ) = (u h (P, kτ ), p h (P, kτ )). The flow passes a unit circular cylinder located at the origin. Let Ω be a bounded computational domain in R 2 considered for the reservoir node preparations and be a set of points that, for the sake of simplicity, were rewritten as figure 1). Let τ = 5 (∈ R) be a transient time in the flow simulation, T f = 10500 (∈ R) be the final time of the flow simulation, T r = T f /τ = 2100 (∈ N) be the total time in our reservoir, and N r = 3S = 231 (∈ N) be the total number of computational nodes in our reservoir. The values of numerical solution at time (k + 1)τ , correspond to the kth input of reservoir, r k ∈ R (k ∈ N ∪ {0}), used most of the time interval [kτ , (k + 1)τ ]. Throughout this study, a one-dimensional random number sequence, {r k } T r k=0 ⊂ R, was used for the input stream z(t) and imposed as the (inflow) Dirichlet boundary condition, u(x, t) = z(t) = (z(t), 0) T on the inflow boundary, in the system, cf the SM for the construction of z(t) from {r k } T r k=0 . The range of r k was fixed at [0.99, 1.01] in the studies discussed in the main text, and multiple different conditions in the input range are investigated in the SM. Using the numerical solution, (u h , p h ), we defined our reservoir {x i k ∈ R; i = 1, . . . , N r , k = 0, . . . , T r − 1} (N r = 3S = 231, T r = 2100) by for = 1, . . . , S(= 77). In short, we employ the state of fluid flow at time (k + 1)τ as our reservoir corresponding to the kth input at time kτ . Then, the outputŷ k was computed withŷ where x 0 k = 1 is a bias and w i out is the readout weight of the ith computational node x i k . As usual, in the framework of RC, the target function, y k = f(r k , r k−1 , . . .), was learned by adjusting the linear readout weights. The ridge regression, known as L 2 regularization [36], was employed to obtain the optimal weights. Details of the training procedures for the weights are given in appendix. The performance of the system output is evaluated by comparingŷ k with the target output y k in the evaluation phase.

Dynamical property of the flow reservoir
In this subsection, we first analyze the bifurcation structure of our reservoir according to the change of Reynolds number and then investigate how the response property to external inputs is regulated according to this bifurcation. Furthermore, the vortex structures, emerging before the onset of the bifurcation point, and its behavior to the external input are investigated in detail.

Hopf bifurcation of the cylinder wake with inputs
The bifurcation structure of the cylinder wake without inputs was first studied theoretically and numerically by Dušek and Le Gal [37], who revealed it to be a Hopf bifurcation. To date, however, the bifurcation structure of a system with inputs has not been clarified. show the behavior of circular cylinder flow using a max/min-velocity diagram and a phase portrait, respectively. We then find that velocity is the periodic type solution at Re = 50 and 100, which means our system has non-chaotic behavior, and the Hopf bifurcation occurs near the critical Reynolds number (around Re = 45). This result is similar to the system without inputs.

Echo state property and common-signal-induced synchronization
When a signal-driven dynamical system is used as a successful reservoir, the dynamical system must accomplish common-signal-induced synchronization [38], which is an outcome of ESP [33]. We analyzed two measures, E m and λ, in each Reynolds number, where E m represents an indicator of synchronization and asymptotically converges to 1 if synchronization occurs, and λ is the amplification factor of perturbation and represents the stability/instability of the flow field (see appendix for the definitions of E m and λ). Figure 3 shows the outcome of the averaged |1 − E m | and λ with each Reynolds number. Both |1 − E m | and λ tend to increase depending on the Reynolds number, suggesting a critical point (∈ [40,50]) where synchronization breaks down during the transition from the twin vortex to the Kármán vortex. The synchronization phenomenon is unlikely to occur under the condition of Re 50 (i.e. the Kármán vortex is generated). However, in the region where synchronization is collapsed (Re ∈ (45, 100]), periodic behavior occurs (figures 2(B) and (C)). Furthermore, the onset of flow instability near the critical Reynolds number has been known to be a manifestation of Hopf bifurcation. Note that the critical point differs from that of chaos in the flow field (chaotic behavior appears at a larger Re (cf, e.g. [37,39])). In our setup, z uses an input with a small swing, such as [0.99, 1.01]; thus, λ is almost the same for systems without input (figure 3). Interestingly, the behaviors of |1 − E m | and λ correspond to the task results presented in figure 5, which will be explained in detail in a later section.

Twin vortex behavior
The twin vortices are known to grow large according to the Reynolds number. We calculated the long diameter of the twin vortices using the stream function ( figure 4(A)). The behavior of this long diameter becomes a brief indicator of the input effects because the long diameter is observed to be oscillating in response to the input stream (see the supplementary videos). To clarify the vortex internal structures of our  system, we focus on the twin vortices (i.e. Re ∈ [10,45]) and, here, analyze the input-sensitivity of the twin vortices. Primarily, the twin vortices are described, but as a comparison, the Kármán vortex cases are also described. Let us introduce the long diameter of the twin vortices: where ψ is the stream function satisfying u 1 = ∂ψ/∂x 2 and u 2 = −∂ψ/∂x 1 . In figure 4(B), the averaged (γ +γ)/2 and max/min γ orγ for time are plotted according to the Reynolds number. We observed the oscillation of the long diameter by the numerical simulation, and then the long diameter was calculated. Subsequently, we determine the expanded long-diameter range with the input effects ( figure 4(B)). Next, we calculated the cross-correlation between inputs r n and the stream function ψ n These coordinates represent the , the plots and error bar show the mean and standard deviation, respectively. We clarified the involvement of the vortex in computational capability by comparing LR, no-object, and ESN with those of the system. We determined that the peak of the computational performance was approximately Re = 40, immediately before the bifurcation point. The background color shading indicates whether the flow was stable (yellow) or unstable (purple). coordinate near the cylinder where the vortices are generated and the approximate coordinates along the long diameter at Re = 45, respectively. In figure 4(C), we used the indicator Cor i (d), i = 1, 2 as follows: where cov(ψ k+d i , r k ) and s(ψ k+d i ) represent the covariance between ψ k+d i and r k and the standard deviation of ψ k+d i , respectively, and the delay d = 0, . . . , 100. In the beginning, there was no significant difference between the input-sensitivity at x min according to the Reynolds number Re ∈ [10,45], but higher Reynolds numbers are associated with slower Cor i (d) convergence. Moreover, the input sensitivity at x max was more significantly different depending on the Reynolds number Re ∈ [10,45]. This indicates that the system had wider effective vortices with respect to the input sensitivity, in proportion to the Reynolds number Re ∈ [10,45], which is an essential property of the vortices reservoir. In the case of the Kármán vortex Re ∈ [50, 100], the input sensitivity was exceedingly low at both points, suggesting that the ESP had collapsed.

Analysis of computational performances
We systematically investigated the computational performance of our reservoir through a number of experiments. As for comparisons, we discuss the results obtained by the flow simulated without the cylinder (called no-object), which can clarify the contribution of vortices to the reservoir dynamics in some cases. Furthermore, we compare the performance of our system with a simple linear regression (LR) model, y k = w 1 LR r k + w 0 LR , which bypassed the reservoir part and directly connected the input to the readout, clarifying the contribution of the reservoir, and with an echo-state network (ESN) [10], which is the standard recurrent neural network (details for the setting of ESN are provided in the SM). For these models, training data, conditions, and procedures are kept the same with our system for fair comparisons.

Evaluation of the memory capacity
Here, we evaluated the memory capacity of the system by investigating whether the system can reproduce previous inputs and nonlinearly process them using its current states. In this study, we applied Legendre polynomials for each time step expressed as y k = P n (r k−d ), and d represents delay given a value, d = 0, 1, . . . , 50: In each plot, the performance of the LR system, system, and no-object are compared with the target. We confirmed that the system can trace accurately through the involvement of a vortex.
The finite products of Legendre polynomials for each time step were used in [12]; however, for simplicity, the nonlinear degree of Legendre polynomials for each delay was changed and treated as a target in the present study, for example, if n = 3, P 3 = 1 2 (5r 3 k−d − 3r k−d ). In this task, we used the memory function and memory capacity as follows: where σ 2 (y k ) represents the variance of y k . The objective of this task was to quantify whether the system could regenerate previous inputs and emulate the nonlinear functions of previous inputs. The quantified values are referred to as linear and nonlinear memory capacity, respectively. The memory capacity showed a high value if the system outputs successfully emulated the target outputs for each target delay. Figure 5(A) shows the results of MC in the evaluation phase. The most striking feature is that MC increased as the system approached the critical Reynolds number and suddenly decreased if the Reynolds number exceeded the critical value. In other words, MC increased as the state of the dynamics neared the move from the twin vortex to the Kármán vortex and significantly decreased after the transition. At Re = 40, immediately before the bifurcation point, the system had the most memory, with MC Re=40 ≈ 9.579 as the actual value of linear memory. Moreover, the extensive twin vortices had more input sensitivity (cf figure 4(C)), and we infer that the implementation of greater computational capability occurred via greater input sensitivity. We confirmed that the highest MC in the system showed a similar value for linear memory as the results for ESN10 (which has ten nodes; this expression applies hereafter), while system performance in the case of nonlinear memory was comparable to ESN100. Furthermore, our system can represent both odd and even functions, whereas some reservoirs, such as ESN, can only demonstrate one of the functions because of the constraint of the setting of activation functions [12].

Emulating nonlinear dynamical systems
We here evaluate whether our system is capable of emulating nonlinear dynamical systems using nonlinear autoregressive moving average (NARMA) models (see appendix for details). The objective of this emulation task was to test the system capability to process input history with nonlinearity. The essential point of the NARMA model was to include the long-term dependencies of the system with nth time lag, where the corresponding NARMA model is called NARMA n model in this study. Figure 5(B) shows the task performance of our system with respect to error in the evaluation phase. The averaged normalized mean squared error (NMSE) was calculated in each NARMA task separately based on each Reynolds number as NMSE = k (y k −ŷ k ) 2 / k y 2 k . We compared the system with no-object and observed that vortices work effectively, depending on the task difficulty and the Reynolds number. We also conducted comparisons of the system based on ESN or LR and observed that the fluid dynamics worked effectively. The results showed that the Reynolds number at optimal performance coincided with the critical Reynolds number of the system ( figure 5(B)). In NARMA2, NARMA3, and NARMA4, system performance at Re = 40 was much better than the result of ESN with 600 nodes. For example, system performance at Re = 10 and 100 was inferior to the performance level at Re = 40 ( figure 5 and the supplementary videos). The plots in figure 6 are representative examples of the actual performance of NARMA2, NARMA3, NARMA4, and NARMA5 at Re = 40 and demonstrate that the system can clearly trace the target model when compared with the LR system (see the supplementary videos). We confirmed that tracking accuracy was improved by the involvement of a vortex, in comparison with that observed in the no-object case. We can infer that nonlinear processing capability was low at Re = 10 because the input sensitivity of vortices was low even though the synchronization occurred and that the computational capability was low at Re = 100 because synchronization was broken.

Time series predictions of missing fluid variables
RC enables highly accurate predictions of time series data and can estimate the Lyapunov exponents of high dimensional spatio-temporal chaotic systems [40,41]. The behaviors of fluid variables [42] or the two-dimensional turbulent Rayleigh-Bénard convection flow [43] has also been predicted/modeled using RC. This phase aimed to implement the prediction of fluid variables (velocity and pressure) using the own dynamics of the system based on the framework of PRC. Dynamics other than the variable to be predicted are used as the reservoir to predict the unknown variable. Thus, if the target output is u 1 , such as, y k = u 1 (P , kτ ), the reservoir {x i k ∈ R; i = 1, . . . , N rr , k = 1, . . . , T r } is defined by for = 1, . . . , S and N rr = 2S = 154. Figure 7 depicts a typical example of an error, with the target output and actual performance for each target. The color map in figure 7 shows that errors increased as the Reynolds number increased; however, the actual performance successfully overlapped the target (figures 7(A)-(C)). As shown in figure 7(D) for each target, although u 2 was moderate, errors of u 1 and p increased remarkably according to the Reynolds number, especially when the Reynolds number exceeded a critical number. This observation indicates that making predictions under conditions where a Kármán vortex is generated is relatively more difficult than when a twin vortex is generated. Figure 7(D) shows that time series prediction was easy with twin vortices in all cases. Time series prediction using the Kármán vortex was difficult with LR and ESN but comparatively easier using the PRC system. The ESP of the system was collapsed in the region Re 50. If the reservoir does not satisfy ESP, its dynamics cannot be uniquely determined by the input signal. Therefore, even if we used optimal weights trained in the training phase, errors increased in the evaluation phase. We can infer that prediction performance was degraded in this Reynolds number region because of flow instability, as the performance of the externally connected ESN model was significantly degraded. Moreover, in the aforementioned memory capacity and NARMA tasks, computational performance may decrease in this Reynolds number region, and we can consider that prediction performance decreases when using PRC. In the region Re 50, the ESN did not work, and the prediction capability decreased when using other variables; however, the prediction capability of the system was better than using the ESN that constructs an input-only function because the target variable was both a function of inputs and a function of other variables that contain some of the target variable components, including its history.

Discussion
In this study, we first investigated the bifurcation of our system and demonstrated that Hopf bifurcation occurs, followed by the synchronization phenomenon to the input stream. We investigated ESP in terms of direct calculation using reservoir dynamics and stability analysis and determined that the synchronization phenomenon gradually collapsed as the phase shifted from a twin vortex to a Kármán vortex. However, we note that this correspondence between the stability analysis and the Reynolds number held true only when input intensity was low and did not necessarily hold true when we used a high input intensity (cf SM). The input sensitivity of vortices was investigated, and the details of its internal structure satisfying ESP were discussed, suggesting that the input sensitivity of dynamics is useful in estimating computational capability. Moreover, we demonstrated that flow dynamics including vortices can be used as a resource for temporal machine learning tasks and that computational performance is maximized in the phase just before synchronization failure. We also showed that in the Reynolds number band where ESP is collapsed, using the remainder of the variables in the same system is more effective for making predictions about the behavior of a missing variable than constructing external prediction models from scratch. This result suggests, in a sense, fluid problems are better solved by fluid computers.
There are still many free parameters that should be investigated in this system. One of the crucial free parameters that affect the information processing capability of the system would be the number of sampling points or the spatial configuration of the sampling points. By altering the sampling locations, according to the vortex structure of the dynamics, the information processing capability would change even if the number of computational nodes is fixed. The choice of the reservoir states also would be important. In our current study, we used (u h , p h ) in each spatial point, but we have found that u h is more important than p h , which is shown in SM. If the number of sampling points (that is related to the number of sensors to be utilized in a physical experiment) is limited, then we can use u h as a much more effective parameter. Several parameters would alter the property of the dynamics drastically. For example, the scaling parameter of the input intensity drastically changes the bifurcation structures of the entire system, as discussed earlier (cf SM). Furthermore, the duration of time τ assigned to each input, which is now fixed to τ = 5, would also affect the dynamics and thus alters the information processing capability of the system when varied. Another important parameter in this context is the shape of the object. This parameter will alter the spatial structure of the vortex and would be crucial for the information processing capability of the system. Although we have focused on the vortex shedding that takes place when a fluid flows past a bluff body, we expect that our scheme presented in this study, investigating the ESP and the common-signal-induced synchronization, can be easily applied to various types of fluid systems and is useful.
Another important next step would be to prove the universality of the proposed system. The universality of the RC system has been studied since the original model of the liquid state machine was proposed [11]. Because RC is targeting temporal tasks rather than static tasks in general, universality in this context means to uniformly approximate any filter with a fading memory property. It was shown that if the system consists of a filter bank, containing a number of filters with fading memory, with a point-wise separation property and of a readout with a universal approximation property, then this system can approximate any filter with a fading memory property in any precision. This result was induced based on the work by Boyd and Chua [44], and recently, there have been interesting theoretical results suggesting that the ESN with a linear and static readout can also give rise to a universal approximator of filters having fading memory in some conditions [45][46][47]. Although it is out of scope of this paper, it should be a fruitful direction to explore in the future.
Considering that the vortices in fluid flow can be found universally in nature and technology, the potential of our findings cannot be overestimated. We expect our findings to be verified in real physical experiments as a future step. Experiments using particle image velocimetry (PIV) would be one of the realistic options for this purpose. Toward the real implementation of the PRC with the fluid flow, it is important to check the effect of noise. We have investigated it in the NARMAn task with n = 2, 3, 4, 6 for a fixed Reynolds number, Re = 40, where the input variance σ 2 takes five values up to 20% (cf SM). The results are natural, that is, the effect of noise becomes large when the noise increases. In addition to the numerical experiments for the noise effect, we also investigated the computational performance of the twin vortex system by using just velocity (or just pressure) in the NARMAn (n = 2, 3, 4, 6) and the memory capacity tasks (cf SM), as the PIV can measure velocity only and the measurement of pressure might not be easy. It is observed that the computational performance of the system with just velocity in the NARMA task is similar to the one with both velocity and pressure as presented and that the pressure contribution is not significant. These additional numerical experiments provide the possibility of physical implementation of the velocity-based PRC via the PIV.
where u = (u 1 , u 2 ) T is the velocity, p is the pressure, σ(u, is the strain-rate tensor, I ∈ R 2×2 is the identity tensor, n : ∂Ω → R 2 is the outward unit normal vector, n ⊥ : ∂Ω → R 2 is the unit tangential vector, z : Γ 1 × (0, T f ) → R 2 is a given boundary velocity, and u 0 : Ω → R 2 is a given initial velocity.
In our system, we set for L = 7.5, Thus, Ω is the computational domain, Γ 1 is an inflow boundary on the left side, Γ 2 is a slip boundary on the top and bottom sides, Γ 3 is a stress-free boundary on the right, and Γ 4 is the no-slip boundary on the circle. We also set z = (z, 0) T for an input function z = z(t) ∈ R, where the input z takes a random value at each fixed time interval (cf SM). The range of z is set as [0.99, 1.01], and experiments with the other ranges are reviewed in the SM. We solved problem (A1) numerically using the stabilized LG method [34,35], a stabilized finite element method combined with the method of characteristics. Moreover, the material derivative was discretized along the trajectory of a fluid particle as where Δt is a time increment for the flow computation, and N T = T f /Δt , t n = nΔt, and u n = u(·, t n ), p n = p(·, t n ) for n ∈ {0, . . . , N T }. Let h > 0 be a representative mesh size. The numerical solution to be obtained by the LG method is a series of piece-wise linear continuous functions defined in Ω, {(u n h , p n h ) : Ω → R 2 × R; n = 1, . . . , N T }. Note that (u n h , p n h ) approximates (u n , p n ). For the fully discretized scheme, please refer to the scheme in the LG method section of the SM.

Appendix B. Training the readout weights using ridge regression
We will briefly review the framework of RC. Assume a reservoir mapF : ∈ N), and a readout maph : R M → R, an infinite discrete inputz = (. . . ,z −1 ,z 0 ,z 1 , . . . ) ∈ (R m ) Z , and an output signalȳ ∈ R Z , then a reservoir state vectorx k ∈ R M and an output signal are determined bȳ x k =F(x k−1 ,z k ) andȳ k =h(x k ), respectively. To exploit a reservoir map as a filter, it is preferable to have its state be a function of the previous input sequence,x k−1 = E(z k−1 ,z k−2 , . . . ), which is related to the concept of ESP, explained in the main text. These expressions are a generalization of the definitions of a reservoir computer, and the readout maph is a linear map. This map, which is expressed as linear and static readout weights, is adjusted to emulate the target outputs.
We used ridge regression [36] following L 2 regularization in the process of learning. The detailed steps are as follows. Ridge regression was used to minimize the residual sum of squares between system and target outputs by penalizing the size of the weights in the training phase. In our setup, the training phase is defined as time domain 200 k 1399 and the optimal weights are w i out , the problem to be solved is where λ R ∈ R 0 is the ridge parameter. By collecting a training data set of 1200 timesteps, we obtained a 1200 × N r matrix X, where N r is the number of nodes in the reservoir. Furthermore, the target output corresponding to 1200 timesteps was determined by y = [y 200 , . . . , y 1399 ] T . Then we obtained the optimal weights, W out = [w 0 out , . . . , w N r −1 out ] T , by We must determine the appropriate λ R , for which the degrees of freedom (df) is crucial. The value of df is determined by where {λ i } are the eigenvalues of X T X. Because df is a one-to-one correspondence with λ R and its maximum value is N r , we can obtain a unique λ R by increasing the df from 1 to N r by one, and the corresponding W out can also be obtained. The residual sum of squares RSS = 1399 k=200 (y k −ŷ k ) 2 can be calculated using W out , and we can obtain Akaike's information criterion (AIC) as follows: where M is the training size of 1200 used in our research. By calculating AIC according to each df, the optimal λ R and W out that minimize AIC can be obtained.

Appendix C. Descriptions of an indicator of synchronization and the amplification factor
In common-signal-induced synchronization, the states of two different initial condition reservoirs driven by the same input sequence approach the same value over time. Intuitively, it indicates that the reservoir asymptotically washes away information related to the initial conditions. This condition indicates that if a certain input sequence is injected at any time, it will return the same response, which expresses the minimum characteristics necessary for the reservoir to implement a reproducible input-output relationship and calculation.
A concept similar to common-signal-induced synchronization is ESP, in which the current network statē r(t) is expressed as a function of only the previous input seriesz(t), independent of the initial valuer(0). Thus, there exists E such thatr(t) = E(. . . ,z(t − 1),z(t)) [33]. In this study, we directly investigated the degree of synchronization using two reservoir states corresponding to two identical input sequences with different initial values and using the amplification factor of perturbation [5], where, to get different initial values, different initial input sequences are used. Let x k andx k (∈ R N r ) be two reservoir states corresponding to two identical input sequences with different initial values; furthermore, we calculated E m as follows: for the norm x k 2 (Ω R ) = { N r i=1 (x i k ) 2 } 1/2 . E m asymptotically converges to 1 if x k andx k systems are synchronized with each other.
Then, the stability analysis of the fluid solution using perturbation was used to investigate whether the perturbation given as an initial value increased or decreased with time. The amplification factor of the perturbation was obtained by numerically solving the perturbation equation as described in the SM, and we obtained the following λ (letũ h = {ũ h } N T n=0 be the perturbation for the velocity obtained by solving the perturbation equation in the SM): , which is the amplification factor of the perturbationũ h . If λ > 0, the perturbation increases exponentially, and the solution is unstable in Ω in this sense. Because this stability analysis is usually performed for steady solutions, we analyzed our setup and the steady solution for comparison, conducting the analysis at z = 1 (called no-input). Notably, the amplification factor of the perturbation, obtained by stability analysis of the fluid solution and the Lyapunov exponent (an indicator used to quantitatively evaluate initial value sensitivity containing a feature of chaos), differ in derivation.

Appendix D. Descriptions of NARMA models
NARMA models used in this study were proposed by Atiya and Parlos [48]. The essential point of the NARMA model was to include the long-term dependencies of the system with nth time lag, which requires learning systems to posses certain amount of nonlinearity and memory to emulate it. We initially introduce the NARMA model of a second-order nonlinear dynamical system as follows: y k+1 = 0.4y k + 0.4y k y k−1 + 0.6r 3 k + 0.1, where y k and r k are the state and input of the model at timestep k, respectively. For the input r k , the same uniform random real value used to the reservoir system was linearly scaled to the range [0, 0.2] and injected.
In this study, this system was called the NARMA2, and the NARMA system such that the nth order nonlinear dynamical system can be written as follows: y k+1 = 0.3y k + 0.05y k ⎛ ⎝ n−1 j=0 y k−j ⎞ ⎠ + 1.5r k−n+1 r k + 0.1.
We have varied n from 3 to 6 and tested the performance. Accordingly, these tasks were called NARMA3, NARMA4, NARMA5, and NARMA6 in the main text.