1. Introduction
Dynamical systems with rhythmic unsteadiness are omnipresent in engineering problems such as structural vibrations, electric circuits, chemical reactions and robotic locomotion (Kuramoto Reference Kuramoto1984; Pikovsky, Rosenblum & Kurths Reference Pikovsky, Rosenblum and Kurths2001). The nonlinear dynamics of such systems is characterized by a self-sustained oscillation, when the state of the system resides sufficiently close to its limit-cycle attractor. The rich physics of limit-cycle oscillators can abundantly be found in periodic fluid flows as well (Huerre & Monkewitz Reference Huerre and Monkewitz1990). These unsteady flows have long been the subject of extensive research for the purpose of prediction, control and extraction of coherent fluid structures (Holmes, Lumley & Berkooz Reference Holmes, Lumley and Berkooz1996). In particular, the study of phase-locking between the flow oscillation and harmonic actuations is essential in revealing the underlying synchronization mechanisms for modifying the dominant oscillation frequency of the flow in a mathematically rigorous and computationally inexpensive fashion (Kawamura & Nakao Reference Kawamura and Nakao2013, Reference Kawamura and Nakao2015; Taira & Nakao Reference Taira and Nakao2018; Iima Reference Iima2019). Within the scope of the present work, synchronization (phase-locking) is defined as the adjustment of the natural frequency of the flow to the frequency of the actuation. In the context of aerodynamic performance of flyers, altering the flow frequency can in turn lead to the enhancement of lift-to-drag ratio, e.g. by advancing or delaying the vortex shedding from the body (Nair et al. Reference Nair, Taira, Brunton and Brunton2020). These objectives require the development of low-dimensional models describing the nonlinear interaction between all oscillators in the fluid system in a simplified fashion (Nair, Brunton & Taira Reference Nair, Brunton and Taira2018).
Most control or reduction techniques conventionally employed for fluid flows rely on the linearization of the flow around a steady or quasisteady state, which makes them underperform when the deviation from the bifurcation point is relatively large (Bewley Reference Bewley2001). Approaches such as Floquet theory (Herbert, Bertolotti & Santos Reference Herbert, Bertolotti, Santos, Dwoyer and Hussaini1987) or phase-reduction theory (Kuramoto & Nakao Reference Kuramoto and Nakao2019) enable the extension of analysis to periodic flows, as they pave a path towards active flow control with respect to time-varying base flows.
The present investigation provides a linear framework for the phase-based description of temporal evolution of time-periodic flows. In particular, we analyse the response of a two-dimensional (2-D) circular cylinder wake to harmonic actuations. The Reynolds number of interest is adequately large so that the von Kármán vortex street forms, but it is also low enough to ensure that the wake remains laminar. This analysis proposes a novel and computationally tractable way to uncover the regions in the flow best suited for the placement of external oscillatory forcing inputs, to achieve synchronization between the vortex shedding and these actuations, and, consequently, to modify the shedding frequency. These regions are obtained by surveying the phase properties of the global near-wake region, significantly expanding the local study performed by Taira & Nakao (Reference Taira and Nakao2018). Furthermore, this methodology is capable of identifying the most sensitive regions of the flow to control inputs at any given time, enabling an energy-optimal control strategy, without any explicit knowledge of the governing dynamics. Additional insights are gained by comparing the findings of the present phase-based analysis against those of other modal- and adjoint-based methods, especially by inspecting the spatial maps of synchronization at different harmonics. The paper is organized as follows. The phase-based, low-dimensional description of the flow is derived in § 2. Section 3 describes the numerical solver used to conduct the direct numerical simulations (DNS) of the flow. Section 4 presents the results of the phase-based model, specifically the regions optimal for synchronization and energy-efficient control, and validates the accuracy of the model predictions for synchronization against DNS results. This section also investigates the potential commonalities and the characteristic differences between the findings of the present model and the widely-used dynamic mode decomposition (DMD), while highlighting the desirable capabilities of the phase-based analysis. Section 5 summarizes the main results and discusses the outlook for subsequent studies.
2. Phase-based analysis of periodic flows
Consider a periodic flow whose dynamics is governed by
\begin{equation} \dot{\boldsymbol{q}} = \boldsymbol{F}(\boldsymbol{q}), \end{equation}
where 
$\boldsymbol {q}$ is the state vector of the flow with a stable limit cycle 
$\boldsymbol {q_0}$ such that 
$\boldsymbol {q_0}(\boldsymbol {x}, t+T) = \boldsymbol {q_0}(\boldsymbol {x}, t)$. Here, 
$T$ denotes the periodicity, through which the natural frequency 
$\omega _n$ can be calculated as 
$\omega _n = 2{\rm \pi} /T$. The periodic nature of the flow allows for phase 
$\theta$ to satisfy
\begin{equation} \dot{\theta} = \omega_n,\quad \theta \in [-{\rm \pi}, {\rm \pi}]. \end{equation}
Note that 
$\boldsymbol {q_0}(\theta (t))$ describes the full flow field at time 
$t$ corresponding to phase 
$\theta (t)$. In the vicinity of the limit cycle, and in the basin of the attractor, the phase dynamics can be captured via the phase function 
${\varTheta }(\boldsymbol {q})$, so that 
$\theta = {\varTheta }(\boldsymbol {q})$. Consequently, we find that
\begin{equation} \dot{\theta} = {\dot{\varTheta}}(\boldsymbol{q}) = \boldsymbol{\nabla}_{\boldsymbol{q}} {\varTheta}(\boldsymbol{q}) \boldsymbol{\cdot} \dot{\boldsymbol{q}} = \boldsymbol{\nabla}_{\boldsymbol{q}} {\varTheta}(\boldsymbol{q}) \boldsymbol{\cdot} \boldsymbol{F}(\boldsymbol{q}) = \omega_n, \end{equation}
which provides a one-dimensional, linear framework for describing the originally nonlinear and high-dimensional dynamics of (2.1) with respect to the limit cycle 
$\boldsymbol {q_0}$.
The flow dynamics of (2.1) can then be weakly perturbed in the following fashion:
\begin{equation} \dot{\boldsymbol{q}} = \boldsymbol{F}(\boldsymbol{q}) + \varepsilon \boldsymbol{f}(t), \end{equation}
where 
$\varepsilon \ll 1$ and 
${\left \lVert \,\boldsymbol {f} \right \rVert }_2 = 1$. This equation can be expressed in terms of the phase variable and function to arrive at
\begin{equation} \dot{\theta} = {\dot{\varTheta}}(\boldsymbol{q}) = \boldsymbol{\nabla}_{\boldsymbol{q}} {\varTheta}(\boldsymbol{q}) \boldsymbol{\cdot} \dot{\boldsymbol{q}} = \boldsymbol{\nabla}_{\boldsymbol{q}} {\varTheta}(\boldsymbol{q}) \boldsymbol{\cdot} [\boldsymbol{F}(\boldsymbol{q}) + \varepsilon \boldsymbol{f}(t)]. \end{equation}Equations (2.5) and (2.3) can then be combined to obtain
\begin{equation} \dot{\theta} = \omega_n + \varepsilon \boldsymbol{\nabla}_{\boldsymbol{q}} {\varTheta}(\boldsymbol{q})\big|_{\boldsymbol{q} = \boldsymbol{q_0}}\boldsymbol{\cdot}\boldsymbol{f}(t). \end{equation}
Here, the higher-order terms have been neglected, assuming that the perturbation is sufficiently small. Hereafter, we denote 
$\boldsymbol {\nabla }_{\boldsymbol {q}} {\varTheta }(\boldsymbol {q}) |_{\boldsymbol {q} = \boldsymbol {q_0}}$ by 
$\boldsymbol {Z}(\theta )$, and refer to it as the phase-sensitivity function.
We find the phase-sensitivity function via a direct method by applying weak impulse perturbations in the form 
$I \delta (t - t_0) \hat {\boldsymbol {e}}_{\boldsymbol {j}}$ at different phase values over a period by changing 
$t_0$, where 
$I$, 
$\delta (t-t_0)$ and 
$\hat {\boldsymbol {e}}_{\boldsymbol {j}}$ represent the perturbation amplitude, the Dirac delta function centred at 
$t_0$ and the unit vector specifying the impulse direction, respectively. In this research, 
$\delta (t-t_0)$ is modelled as a narrow Gaussian function given by
\begin{equation} \delta(t-t_0) = \frac{1}{\sqrt{2{\rm \pi}} \sigma} \exp\left[ - \frac{1}{2}\left(\frac{t - t_0}{\sigma}\right)^2 \right], \end{equation}
where 
$\sigma$ is taken as 
$10 {\rm \Delta} t$, with 
${\rm \Delta} t$ indicating the DNS time step (see § 3).
The introduction of impulse perturbation results in an asymptotic phase advancement or delay. This asymptotic phase shift, known as the phase-response function 
$g(\theta; I\hat {\boldsymbol {e}}_{\boldsymbol {j}})$, can be measured experimentally or numerically for a certain phase value and location at which the perturbation is added. Given that 
$I \ll 1$, the phase-sensitivity function can be evaluated as
\begin{equation} Z_j(\theta) = \lim_{I \rightarrow 0} \frac{g(\theta; I\hat{\boldsymbol{e}}_{\boldsymbol{j}})}{I} \approx \frac{g(\theta; I\hat{\boldsymbol{e}}_{\boldsymbol{j}})}{I}. \end{equation}
The phase-sensitivity functions can be computed for various phase values and locations, to obtain the 
$2{\rm \pi}$-periodic phase-sensitivity functions of the entire flow field. One can alternatively integrate the adjoint equations of the dynamical system backward in time to obtain 
$\boldsymbol {Z}$ (Ermentrout & Terman Reference Ermentrout and Terman2010; Nakao Reference Nakao2016). This adjoint-based approach relies on the calculation of Jacobi matrix of 
$\boldsymbol {F}(\boldsymbol {q})$ at each time, which requires explicit knowledge of the governing equations and, as a result, cannot be implemented experimentally.
Once 
$\boldsymbol {Z}$ is available, the linear model of (2.6) can be used for designing optimal open- or closed-loop controls (Nair et al. Reference Nair, Taira, Brunton and Brunton2020). It can also be utilized to study the synchronization of the periodic flow to an external actuation such as 
$\varepsilon \boldsymbol {f}(t)$ with the frequency 
$\varOmega$, assuming that 
${\left \lvert \varOmega - \omega _n \right \rvert } \ll 1$. We define the relative phase between the actuation and the flow oscillation as 
$\phi (t) = \omega _n - \varOmega t$. This definition can be combined with (2.6) to calculate the temporal rate of change of 
$\phi$,
\begin{equation} \dot{\phi} = \varepsilon[\varDelta + \boldsymbol{Z}(\phi + \varOmega t) \boldsymbol{\cdot} \boldsymbol{f}(t)], \end{equation}
where 
$\varDelta = (\omega _n - \varOmega )/\varepsilon = {O}(1)$.
The synchronization or phase-locking between the flow oscillation and harmonic forcing occurs when 
$\dot {\phi } \rightarrow 0$. As a consequence, and after many cycles, the flow adjusts its frequency, so that it oscillates with the new frequency 
$\varOmega$. Since 
$\varOmega t$ changes much more rapidly than the relative phase (
${\left \lvert \varOmega - \omega _n \right \rvert } \ll 1$), the right-hand side of (2.9) can be estimated by taking its average over one period (Kuramoto Reference Kuramoto1984; Ermentrout & Terman Reference Ermentrout and Terman2010). Hence, (2.9) can be simplified to
\begin{equation} \dot{\phi} = \varepsilon[\varDelta + {\varGamma}(\phi)], \end{equation}
where the 
$2{\rm \pi}$-periodic phase-coupling function 
${\varGamma }(\phi )$ is formulated as
\begin{equation} {\varGamma}_m(\phi) = \frac{\varOmega}{2{\rm \pi}} \int_{0}^{2{\rm \pi}/\varOmega} \boldsymbol{Z} (\phi + \varOmega t'/m ) \boldsymbol{\cdot} \boldsymbol{f}(t')\,\mathrm{d}t' = \frac{1}{2{\rm \pi}} \int_{0}^{2{\rm \pi}} \boldsymbol{Z}(\phi + \varphi/m) \boldsymbol{\cdot} \boldsymbol{f}(\varphi) \,\mathrm{d}\varphi, \end{equation}
while 
$\varOmega = m\omega _n$ with 
$m \in \mathbb {N}$. Stable phase-locking between the external forcing and the periodic flow can be achieved if and only if
\begin{equation} \varepsilon {\varGamma}_{{min}} < \varOmega - \omega_n < \varepsilon {\varGamma}_{{max}}, \end{equation}
where 
${\varGamma }_{{max}}$ and 
${\varGamma }_{{min}}$ are respectively the maximum and minimum values of the phase-coupling function at a certain point. The region of synchronization in the frequency–amplitude plane, also called the Arnold tongue (Arnold Reference Arnold1997), is thus bounded by two lines with the slopes 
$-1/\varGamma _{{min}}$ and 
$1/\varGamma _{{max}}$, both passing through the point 
$(\varOmega /\omega _n, \varepsilon ) = (1, 0)$ (e.g. see figure 5). Evidently, as the difference between 
${\varGamma }_{{max}}$ and 
${\varGamma }_{{min}}$ becomes larger, the area of the synchronization region grows as well. We therefore define the parameter synchronizability 
$S \equiv {\varGamma }_{{max}} - {\varGamma }_{{min}}$ as a measure for the convenience of phase-locking. While the above derivations assume 
$\varepsilon \ll 1$, the synchronization criterion (2.12) often holds for moderate values of 
$\varepsilon$.
3. Numerical simulations
Let us consider the application of the phase-reduction analysis to characterize the synchronization properties of the circular cylinder wake. The governing equations for this 2-D incompressible flow are
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \cdot \boldsymbol{\nabla}{\boldsymbol{u}} = \frac{1}{Re} \nabla^2 \boldsymbol{u} - \boldsymbol{\nabla} p + \varepsilon \boldsymbol{f} \quad \textrm{and} \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{equation}
where 
$\boldsymbol {u} = u \hat {\boldsymbol {e}}_{\boldsymbol {x}} + v \hat {\boldsymbol {e}}_{\boldsymbol {y}}$, 
$p$ and 
$\boldsymbol {f}$ represent dimensionless velocity, pressure and external forcing, respectively. Equations (3.1a,b) are non-dimensionalized by taking the cylinder diameter 
$d$ and the freestream velocity 
$U$ as the characteristic length and velocity (figure 1a). Here, we study the flow at a Reynolds number of 
$Re \equiv Ud/\nu = 100$, where 
$\nu$ denotes the kinematic viscosity. For the cylinder flow under consideration, the state of the flow is determined by measuring the lift coefficient 
$C_L$ and its temporal derivative 
$\dot {C}_L$ (figure 1b), whose asymptotic changes in response to the impulse perturbations give the phase response and, subsequently, the phase-sensitivity function of the wake flow.
(a) Set-up of impulse perturbations and harmonic actuations used for calculating the phase properties of the fluid flow. The quantities 
$\alpha$ and 
$r$ denote the angular position of the perturbation and its distance from the centre of the cylinder, respectively. Horizontal and vertical black arrows indicate the position and orientation of impulse perturbations used for developing the linear, phase-based model, while the red arrows show example arrangements of top and bottom actuators used to study the synchronization properties of the wake flow. (b) Limit-cycle attractor of the flow depicted in the 
$\dot {C}_L$–
$C_L$ plane, along with the definitions of the phase variable 
$\theta$ and the phase function 
${\varTheta }$.
A filter in the form of the three-point discrete delta function of Roma, Peskin & Berger (Reference Roma, Peskin and Berger1999) is multiplied by the forcing 
$\boldsymbol {f}$ of (3.1a,b) to ensure its spatial compactness. This is enforced in the DNS solver by the following integral over the entire domain:
\begin{equation} \boldsymbol{f}(\boldsymbol{x}, t) = \sum_{j = 1}^{n_f} \int_{\boldsymbol{x^{\prime}}} \delta_x(\boldsymbol{x^{\prime}} - \boldsymbol{x}_j)\hat{\boldsymbol{e}}_{\boldsymbol{j}} T(t) \,\mathrm{d}\kern0.6pt \boldsymbol{x^{\prime}}, \end{equation}
where 
$n_f$, 
$\delta _x$ and 
$\boldsymbol {x}_j$ respectively show the number of forcings, the discrete delta function and the location at which the forcings are applied. The temporal part of the forcing 
$T(t)$ represents the impulse perturbations 
$I \delta (t - t_0)$ of § 2 or the periodic forcings of the next section. Such impulsive, point forces are placed at different angular positions 
$\alpha$ and distances 
$r$ from the centre of the cylinder, to map out the phase-sensitivity function for the entire flow field (figure 1a). The perturbation amplitudes are selected so that 
$I \in [0.01, 0.02]$, making them adequately small to satisfy the linearity assumption.
The DNS of governing equations (3.1a,b) is carried out using the immersed boundary projection method of Taira & Colonius (Reference Taira and Colonius2007). The far-field boundary conditions are handled via the multi-domain technique of Colonius & Taira (Reference Colonius and Taira2008), while each domain is discretized uniformly. The computational domain is chosen as 
$x/d \in [-39.5, 40.5]$ and 
$y/d \in [-40, 40]$ with the cylinder located at the origin. The first domain, which is the nearest to the cylinder, lies within the area 
$[-2, 3]\times [-2.5, 2.5]$ and has the grid spacing 
${\rm \Delta} x/d = {\rm \Delta} y/d = 0.025$. The time step 
${\rm \Delta} t$ is set to satisfy 
${\rm \Delta} t < 0.5{\rm \Delta} x/U$. When the impulse perturbations are located too close to the boundaries of the first domain or fall outside of it, the domain is sufficiently enlarged. The natural frequency, the time series of lift and drag coefficients, and the velocity and vorticity fields rendered by the present computational approach are in agreement with those of earlier studies such as Munday & Taira (Reference Munday and Taira2013) and Taira & Nakao (Reference Taira and Nakao2018), validating the overall set-up.
4. Results
In this section, we discuss the results regarding the maps of phase-sensitivity functions and synchronizability for the entire flow field. Such findings reveal valuable information about designing optimal control strategies for the frequency or phase modification of the wake dynamics. The accuracy of the model predictions for the synchronizability are validated against DNS results. We furthermore study the connections and differences between the present method and the DMD-based analysis of the flow, in order to illustrate the distinctive features of the phase-reduction analysis.
The spatial maps of phase-sensitivity functions demonstrated in figure 2 exhibit a symmetry condition between the upper and lower halves of the flow, so that for horizontal impulses 
$Z_x(\theta; x/d, \alpha ) = Z_x(\theta + {\rm \pi}; x/d, -\alpha )$, and for the vertical ones 
$Z_y(\theta; x/d, \alpha ) = -Z_y(\theta + {\rm \pi}; x/d, -\alpha )$. In accordance with the findings of Luchini, Giannetti & Pralits (Reference Luchini, Giannetti, Pralits, Braza and Hourigan2009) and Iima (Reference Iima2019), we observe that the sensitivity to the forcing inputs is typically the largest downstream of the cylinder and in its vicinity. In addition, a horizontal perturbation in the positive direction and placed upstream of the cylinder either does not change the phase 
$\theta$ (
$Z_x = 0$) or causes phase advancement (
$Z_x < 0$), as also reported by Iima (Reference Iima2019). The maps of figure 2 indicate that an energy-efficient control strategy for shifting the phase of the flow requires placing the control inputs around the separation point (
${\rm \pi} /3 \lesssim {\left \lvert \alpha \right \rvert } \lesssim {\rm \pi}/2$ and 
$0.6 \lesssim x/d \lesssim 0.8$) when they are oriented horizontally, and much closer to the wake (
${\rm \pi} /12 \lesssim {\left \lvert \alpha \right \rvert } \lesssim {\rm \pi}/4$ and 
$0.6 \lesssim x/d \lesssim 0.8$) when they are vertical. The stated regions correspond to the peaks of 
$Z_x$ and 
$Z_y$.
Phase-sensitivity function at different phases 
$\theta$, given by impulse perturbations in (a) 
$\hat {\boldsymbol {e}}_{\boldsymbol {x}}$ and (b) 
$\hat {\boldsymbol {e}}_{\boldsymbol {y}}$ directions.
In this study, we consider the synchronization of the wake to two sinusoidal forcing inputs in the form 
$f = 0.5\varepsilon [1 + \sin (\varOmega t)]$, placed at the top and bottom of the cylinder, as shown in figure 1. Note also that there is a phase difference 
${\rm \pi}$ between the two actuators. The symmetric distributions of synchronizability 
$S$ with respect to the wake centreline are depicted at different harmonics in figure 3, when 
$\boldsymbol {f}$ is directed horizontally and vertically. As can be seen in this figure, for the odd-harmonic forcings, i.e. when 
$\varOmega = (2k-1)\omega _n$ (
$k \in \mathbb {N}$), horizontal actuations give much larger values for 
$S$ in comparison with vertical actuations, suggesting that for these forcings, phase-locking to horizontal actuations can take place with much lower forcing amplitudes. For the even-harmonic forcings, however, the opposite holds, meaning that synchronization to vertical actuations at these harmonics requires lower values for 
$\varepsilon$. Furthermore, the regions with large synchronizability in the first harmonic mainly overlap with those of the associated phase-sensitivity map. In fact, the location of maximum synchronizability for horizontal (vertical) actuations at 
$\varOmega = \omega _n$ is found to be at 
$x/d = 0.6$ and 
$\alpha = 5{\rm \pi} /12$ (
$x/d = 0.6$ and 
$\alpha = {\rm \pi}/4$), marked by the green point of the corresponding map. For vertical periodic forcings, the general shape of synchronizability maps does not change substantially for different harmonics, and the maximum value of 
$S$ appears at the same location for all. In contrast, when the forcings are horizontal, the pattern of 
$S$ can significantly differ among different harmonics. Specifically, when 
$\varOmega = 2k\omega _n$, 
$S$ in the centreline of the wake increases from approximately zero at the first harmonic to some fairly large values in 
${O} (1)$. This can be attributed to the 
${\rm \pi}$-periodicity of 
$Z_x$ along the centreline of the wake, leading to the orthogonality of the phase-sensitivity function to 
$\boldsymbol {f}$, and nearly zero values for both 
${\varGamma }_{{max}}$ and 
${\varGamma }_{{min}}$ when 
$\varOmega = \omega _n$. At 
$\varOmega = 2\omega _n$ or other higher harmonics, this orthogonality no longer holds, as can be understood from the formulation for the phase-coupling function (2.11), resulting in relatively large values for 
$S$ at 
$\alpha = 0$. In fact, the maximum value of 
$S$ occurs at 
$\alpha = 0$ and 
$x/d = 1.2$ (
$x/d = 1.4$) for 
$\varOmega = 2\omega _n$ (
$\varOmega = 4\omega _n$). The synchronization of actuations placed at the coloured points shown in the synchronizability maps of the first harmonic are compared against DNS results in figure 5.
(a) Synchronizability fields at different harmonics for horizontal periodic actuators. The green point in each panel indicates the position at which 
$S$ for a certain value of 
$\varOmega /\omega _n$ is maximized. (b) Same as (a), except harmonic actuators are directed vertically.
The comparison between the findings of the present phase-based approach, particularly for synchronizability maps, and the DMD analysis of the wake flow underlines the characteristic differences between the two methods. Most distinctively, the phase-based model finds the optimal location for the synchronization fairly close to the cylinder, while the peaks of DMD modes fall at least two diameters away from the cylinder, as can be seen for the velocity modes of figure 4. We should emphasize that the phase-reduction and DMD analyses are employed for different objectives. The former identifies the optimal location for an oscillatory control input capable of modifying the phase or oscillation frequency of the flow, whereas the latter aims to extract the most dynamically relevant structures of the flow field. The comparison also exhibits some commonalities between the two. Most noticeably, it shows that the patterns identified by DMD modes are typically strong along the same directions in which 
$S$ is large. For instance, the structures in DMD modes of the horizontal velocity field are vanishing at 
$\alpha = 0$ for the first and third modes, but are the strongest for the second and fourth modes. In addition, for the odd modes of vertical velocity, it is observed that the dominant structures emerge closer to the wake centre in comparison with those of 
$u$-modes. These findings are in line with the results of the present analysis for synchronizability fields.
Magnitude of DMD modes for (a) horizontal and (b) vertical velocity fields, with the frequency of mode 
$m$ being 
$m\omega _n$.
The present phase-based model reveals the necessary conditions for the synchronization of the vortex shedding to external actuations. The theoretically calculated Arnold tongues corresponding to the points A to D of figure 3 (solid lines in figure 5) are compared against the numerically obtained boundaries of synchronization in the frequency–amplitude plane for the same cases (blue crosses in the same figure). The numerical boundaries for a prescribed value of 
$\varepsilon$ are found by locating the largest value of 
${\left \lvert \varOmega - \omega _n \right \rvert }$ for which the difference between the oscillation frequency of the flow and 
$\varOmega$ is below 
$1\,\%$. Except for the case with horizontal actuation at 
$x/d = 0.6$ and 
$\alpha = {\rm \pi}/2$ (point C in figure 3a), in which model predictions slightly deviate from DNS results, all other cases show excellent agreement between the two. The linearity of the phase-based model requires a careful parametric study for the perturbation amplitude 
$I$, to ensure that 
$I$ is adequately small so that its further decrease does not change the phase-sensitivity function noticeably. This step has significant impact on the overall accuracy of the model predictions. Furthermore, as the amplitude of the excitation grows, the state of the flow moves farther from its stable limit cycle, and as a consequence, nonlinear effects become gradually more pronounced; however, the results of figure 5 evidence that the model remains accurate for a broad range of forcing amplitudes. Consistent with the observations of Rigas, Morgans & Morrison (Reference Rigas, Morgans and Morrison2017) and Herrmann et al. (Reference Herrmann, Oswald, Semaan and Brunton2020), we also find that subharmonic synchronization (i.e. when the vortex shedding oscillates with the frequency 
$\varOmega /m$ instead of 
$\varOmega$) can commonly arise at higher harmonics. However, since the smaller values of 
$S$ at these harmonics (especially at the even harmonics of the horizontal synchronizability fields) make them difficult to achieve, phase-locking at higher harmonics is not shown here.
Synchronization of vortex shedding to harmonic actuations for (a) horizontal and (b) vertical forcing inputs, located at points A to D of figure 3. The phase-based analysis suggests that only cases falling inside the areas bounded by the solid black lines meet the theoretical criteria for phase-locking. The numerically found boundaries of synchronization are also marked by the blue crosses.
Once the maps of 
$Z_x$ and 
$Z_y$ are available, the phase-sensitivity function for any arbitrarily-oriented impulse perturbation at any point can be calculated as 
$Z(\beta ) = Z_x \cos (\beta ) + Z_y \sin (\beta )$, where 
$\beta$ is the angle between the forcing direction and 
$\hat {\boldsymbol {e}}_{\boldsymbol {x}}$. Subsequently, other phase properties such as 
$\varGamma$ and 
$S$ can also be computed for all possible locations and directions of forcing inputs, following the formulations of § 2. The optimal orientation of actuations leading to the maximum synchronizability for different points in the flow field are shown in figure 6. The length of each vector is proportional to the magnitude of 
$S$, when a harmonic forcing in the same direction and at the same location is introduced. The optimal orientation for the periodic actuations in the near wake (
$x/d \lesssim 1.5$) is found to be nearly tangential to the cylinder. This is of practical significance, since tangential forcings can result in a more streamlined and stabilized wake, causing substantial drag reduction (Munday & Taira Reference Munday and Taira2013). It is also seen that actuations tangential to the cylinder and placed at 
$x/d = 0.6$ and 
$\alpha = \pm 5{\rm \pi} /12$ (red arrows of figure 6a) yield the largest value possible for 
$S$. Therefore, a harmonic actuation at the same location as the red vectors and parallel to them results in phase-locking with the smallest amplitude 
$\varepsilon$, in comparison with all other actuations at any other locations and with any other orientations. The synchronizability field corresponding to periodic forcings parallel to the top red arrow of figure 6(a), and with 
$\varOmega = \omega _n$, is remarkably similar to the synchronizability map of horizontal forcings at the same harmonic (compare figures 6b and 3a), except the values of 
$S$ in figure 6(b) are around 
$5\,\%$ larger.
(a) The vector field represents the optimal orientation of forcing inputs at different locations in the flow field, in order to acquire the maximum synchronizability. The length of each vector is proportional to the synchronizability achieved by applying the perturbation in the same direction and at the same location as the vector. The red vectors correspond to the actuations leading to the maximum value of 
$S$. (b) Synchronizability field at 
$\varOmega = \omega _n$, when the periodic actuations are aligned with the top red vector of panel (a).
The capabilities of the present phase-based analysis in modifying the oscillation frequency of the flow can give insights into the physics of the vortex shedding process, as well as the interaction mechanisms of vortical structures in the wake flow, such as vortex splitting and merging. Furthermore, this method proposes a novel approach for controlling the unsteady dynamics of periodic flows through changing their oscillation frequency or phase by finding the optimal strategy for the placement of forcing inputs. Such a strategy can in turn enable substantial drag reduction and/or lift enhancement. The current framework of the phase-reduction theory holds promise for its application to a variety of engineering problems, such as biological swimmers and flyers, vortex-induced vibration and fluid-structure interactions.
5. Concluding remarks
We performed a phase-based analysis on the 2-D periodic flow past a cylinder at 
$Re = 100$. Weak impulse perturbations were introduced at different locations in the flow field and at various times over a period, in order to uncover the important phase properties of the periodic wake. The present direct method does not need a priori knowledge of the underlying dynamics, is computationally tractable, and can be straightforwardly implemented in experimental set-ups. The maps of the phase-sensitivity function and synchronizability identified by the phase-based model of (2.6) reveal optimal strategies for altering the wake dynamics via modification of the phase or frequency of the flow. The present method accurately predicts the necessary conditions for the synchronization of the wake flow to external harmonic actuations, as confirmed by comparing the model predictions against DNS results. Furthermore, the comparison between the phase- and DMD-based analyses of the flow shows that, despite some interesting similarities, the optimal positioning of forcing inputs required for changing the frequency or the phase of the flow is better determined using the present phase-based approach. This should not come across as surprising, as the two approaches have been developed for vastly different objectives. In subsequent studies, the authors aim to investigate the commonalities and discrepancies between the two methods in detail, consider more complex physics such as the interactions between the periodic wakes of multiple bodies, and extend the current approach to quasiperiodic flows, e.g. via cluster-based techniques.
Acknowledgements
We gratefully acknowledge the support of the US Air Force Office of Scientific Research (grant no. FA9550-16-1-0650; programme managers: Drs G. Abate and D. Smith) and the Army Research Office (grant no. W911NF-19-1-0032; programme manager: Dr M. J. Munson).
Declaration of interests
The authors report no conflict of interest.
