Modeling synchronization in forced turbulent oscillator flows

Periodically forced, oscillatory fluid flows have been the focus of intense research for decades due to their richness as a nonlinear dynamical system and their relevance to applications in transportation, aeronautics, and energy conversion. Recently, it has been observed that turbulent bluff-body wakes exhibit a subharmonic resonant response when excited with specific spatial symmetries at twice the natural vortex shedding frequency, which is hypothesized to be caused by triadic interactions. The focus of this paper is to provide new physical insight into the dynamics of turbulent oscillator flows, based on improved mechanistic models informed by a comprehensive experimental study of the turbulent wake behind a D-shaped body under periodic forcing. We confirm for the first time the role of resonant triadic interactions in the forced flow by studying the dominant components in the power spectra across multiple excitation frequencies and amplitudes. We then develop an extended Stuart-Landau model for the forced global wake mode, incorporating parametric and non-harmonic forcing. This model captures the system dynamics and reveals the boundaries of multiple synchronization regions. Further, it is possible to identify model coefficients from sparse measurement data, making it applicable to a wide range of turbulent oscillator flows. We believe these generalized synchronization models will be valuable for prediction, control, and understanding of the underlying physics in this ubiquitous class of flows.


I. INTRODUCTION
Fluid flows that display unsteadiness characterized by a well-defined frequency and that are insensitive to lowlevel external noise are known as oscillator flows [1,2]. These flows have been the focus of research efforts for over 75 years [3,4], in part because of their rich physics, and also because of their relevance to numerous applications where aerodynamic forces and mixing play a significant role, such as transportation, aeronautics, and energy conversion [5]. Models that capture the evolution of dominant fluid coherent structures are of the utmost importance for prediction, control, and understanding of the underlying physical processes that drive these flows [6]. The wake past a bluff body is one example of an oscillator flow, where self-sustained periodic vortex shedding arises after an increase in the Reynolds number renders the flow incapable of maintaining a steady state. This scenario unfolds when a supercritical Hopf bifurcation takes place, where disturbances associated with a spatial structure, known as the global mode of the flow, become linearly unstable, leading to exponential growth of the mode amplitude A, followed by nonlinear saturation onto a stable limit cycle [7,8]. The Stuart-Landau model [3,4] has been widely used to explain this nonlinear oscillator behavior of the wake past bluff bodies [7][8][9][10][11][12]. Sipp and Lebedev [7] formally derived this model from the Navier-Stokes equations by means of a rigorous asymptotic expansion close to the Hopf bifurcation. When periodically forced with certain frequencies, the wake past a bluff body has been observed to adjust its * benherrm@uw.edu natural vortex shedding rate to some rational multiple of the forcing frequency, as first reported by Provansal et al. [13] for the cylinder flow. This is synchronization -the spontaneous emergence of rhythmic oscillatory dynamics -an inherently nonlinear phenomenon that is abundant in natural and engineering systems such as chemical reactions, electric circuits, structural vibrations, cardiac cells, spiking neurons, and the locomotion of animals and robots [14][15][16][17][18][19][20]. In the context of fluid dynamics, the recent work of Taira and Nakao [21] was the first to study the synchronization properties of the cylinder flow using phase-reduction analysis, a technique commonly used for biological and chemical systems [22]. Periodic forcing continues to present an appealing flow control strategy for a wide range of applications, as it has been shown to effectively reduce bluff body drag [23], increase lift of airfoils [24], and enhance mixing in heat exchangers [25]. Understanding synchronization in the context of unsteady aerodynamics is key to leverage periodic flow control and explain the mechanisms that lead to the performance improvements observed in these success stories.
Recent experimental work by Barros et al. [26] and Rigas et al. [27] reported synchronization in the harmonically forced turbulent wake past an Ahmed body and an axisymmetric blunt body, respectively. Both studies linked the spatial symmetry properties of the forcing mode to the type of response observed. They found the presence of a 1 : 2 subharmonic resonance for symmetric disturbances, where the vortex shedding frequency synchronizes to half of the forcing frequency. This lockon phenomenon was attributed to resonant wave-triads, where the spatial structure of the interacting forcing and response modes are constrained by the triadic consistency condition [28][29][30]. This hypothesis suggests that studying the synchronization properties of forced flows past bluff bodies is a promising way to learn about the dominant nonlinear mode interactions present in the wake.
In previous work, the response of oscillator flows to periodic disturbances has been modeled using the Stuart-Landau equation with the addition of a forcing term [13,27,[31][32][33]. In the work of Sipp [32], this model was derived analytically by extending the weakly nonlinear analysis of laminar globally unstable flows to include the effects of external forcing. Rigas et al. [27] used an eddy viscosity closure and the phase-averaged Navier-Stokes equations to extend the analysis to the turbulent regime. They validated the resulting model in the proximity of a subharmonic resonance by comparing against experiments of the turbulent wake past an axisymmetric body. Hence, weakly nonlinear analysis [27,32] provides a theoretically-based structure for a model of the form where g(A, F ) captures the excitation induced by the coupling between the global mode, with amplitude A, and the periodic forcing mode, with amplitude F . However, in practice, the computation of g(A, F ) requires high-fidelity numerical simulations. Alternatively, recent data-driven techniques are enabling the identification of models directly from data [34][35][36][37][38][39][40][41][42][43]. Our goal in this work is to obtain new physical insight into the dynamics of forced turbulent oscillator flows by deriving improved mechanistic models from more comprehensive experimental data than has been previously reported. As an example of this class of systems, we study the turbulent wake behind a D-shaped bluff body, subject to periodic Coanda blowing, using wind tunnel experimental data, as shown in Fig. 1 and described in §II. In §III we present our experimental dataset characterizing the wake response to periodic forcing for a wide range of excitation parameters. We elucidate the role of resonant triadic interactions for the first time by studying the finely resolved variations of the power spectral density of the global mode response with the excitation frequency. In §IV we derive an extended Stuart-Landau model for the evolution of the forced global mode that is amenable to analysis and explain how its coefficients can be identified directly from data. In §V, we show that our model accurately captures multiple resonances and frequency lock-on regions observed experimentally. Furthermore, model predictions reveal the boundaries of the synchronization regions for the forced D-shaped bluff body wake. Concluding remarks are offered in §VI.

II. EXPERIMENTAL SETUP
A major contribution of this work is the exhaustive experimental investigation of the turbulent wake response to periodic blowing for a broad range of excitation fre-quencies, excitation amplitudes, and forcing configurations. The experiments are conducted in the "Leiser Niedriggeschwindigkeitswindkanal Braunschweig" (LNB) wind tunnel at the institute of fluid mechanics of the Technische Universität Braunschweig. The LNB is a continuous atmospheric Eiffel type open return wind tunnel with room recirculation and a closed test section. It has a Burger-type nozzle with a contraction ratio of 16 : 1. To reduce turbulence, incoming air is guided through a 30mm thick fleece mat, a 133mm thick honeycomb, and a fine woven screen. The resulting turbulence level is below 0.1% at 10m/s. The test section has a width of 400mm, a height of 600mm, and a length of 1500mm. To compensate for boundary layer growth on the walls, the test section has a horizontal opening angle of 1 degree. The flow is driven by a 9-blade fan at the end of the diffuser.
The experimental model is a D-shaped bluff body with a blunt trailing edge; it has a height of 53.4mm, a length of 190.6mm, and a width of 390mm. The model is horizontally mounted in the wind tunnel and held by one steel tube on each side. The model nearly spans the entire width of the test section. Zigzag tape is applied to the upper and lower side at about 9% body length to trip the boundary layer and to prevent the formation of a laminar separation bubble. A sketch of the model is presented in Fig. 1(a). The model is equipped with two Coanda actuators at the trailing corners, each fed by four plenum chambers. The Coanda surfaces have a 9.4mm radius, which was determined by numerical optimization [44]. The jet slit height is set to 0.2mm. The uniformity of each jet is verified with a fish mouth probe to be within 10% of the mean exit pressure.
Unsteady actuation is enabled through eight Festo MHJ9-QS4-MF monostable 2/2-way valves with an operating pressure range of 0.5 to 6bar. The valves can be operated at maximum frequency of 1kHz. Time-resolved pressure signals are acquired by five Honeywell SLP pressure sensors distributed along the mid-span of the rear face of the model. The sensors have a measurement range of ±1000Pa differential pressure, a repeatability of 0.5% of the full scale, and a response time of 100µs. Plenum pressure is monitored by two Kulite pressure sensors with a range of ±3.5 · 10 4 Pa differential pressure and an accuracy of ±0.1% of the full scale. All differential pressures are measured relative to the static pressure in the free stream. The instantaneous jet velocity is estimated from the pressure measurements in the plenum chambers.

III. EXPERIMENTAL DATASET
The dynamics of turbulent bluff body wakes exhibit self-sustained oscillations. The dominant coherent structure characterizing the oscillatory motion is known as the global vortex shedding mode. We use the antisymmetric average of the five base pressure sensors x(t) to capture the amplitude of the global vortex shedding mode, as shown in Fig. 1(a). For the unforced flow, x(t) has a power spectral density with one distinct peak located at the fundamental shedding frequency ω 0 , corresponding to a Strouhal number based on the free-stream velocity and body height of St = ω 0 H/2πU ∞ = 0.23, as shown in Fig. 1(b). We study the long-term response of the global vortex shedding mode to periodic forcing, for a broad range of parameters, including the forcing configuration, and the excitation amplitude and frequency. Two forcing configurations are considered: in-phase blowing through the top and bottom slits, leading to the excitation of a spatially symmetric flow structure, and 180 • out-of-phase blowing, exciting a spatially antisymmetric structure. The effect of forcing amplitude is investigated for three blowing intensities quantified non-dimensionally by the momentum coefficient as where U jet is the root mean square value of the jet velocity calculated from the plenum pressure measurements, h = 0.2mm is the slot height, H = 53.4mm is the body height, and the factor 2 accounts for the number of actuators. We define the excitation amplitude ε as the average of c µ over all excitation frequencies at a given tank pressure.
For each forcing configuration and amplitude, the base pressure is recorded for 8s at a sampling rate of 5kHz for blowing frequencies between 1Hz and 210Hz with 1Hz intervals. For all cases, the conditions are maintained for 5s to ensure steady state behavior before the measurements are taken. The complete dataset is conformed by a total of 1260 time-series of x(t), each consisting of 40000 samples. These time-series are used to characterize the mean frequency and amplitude of the wake response. Each time-series x(t) is low-pass filtered using a 5 th -order Butterworth filter with a cutoff frequency of 1.3ω 0 . Its Hilbert transform x H (t) is then computed to build an analytic signal for the complex global mode amplitude A(t) = (x + ix H )(t), which is a common practice when studying oscillatory dynamics from data [20]. Once we have the complex time-series A(t), we compute its mean amplitude r, and its mean frequency ω from the mean of the time derivative of its instantaneous phase.
The mean wake response provides information about resonances and lock-on regions that has already been discussed is in previous studies [26,27]. A deeper insight into the underlying nonlinear interactions is obtained by examining the power spectral density (PSD) of the global mode response at different forcing frequencies. The normalized PSD of x(t) is computed using Welch's method [45], splitting the time-series into 10 segments with 50% overlap and tapered by a Hanning window. The resulting PSD obtained for each excitation frequency are stacked as columns in a heat map, shown in Fig. 1(c). The large number of test cases yields a fine resolution of the effect of ω f on the PSD of the global vortex shedding mode, resulting in a highly interpretable visualization that we refer to as the power spectrum response. An enlarged version of the power spectrum response for the case of symmetric forcing is shown in Fig. 2(a). Identification of the dominant frequency components of the response, Ω, provides insight into the nature of the interactions between forcing and response modes. The highest energy component coincides with the observed mean vortex shedding rate, Ω = ω. In addition, Fig. 2(a) displays components following sharp straight lines, each with a physical interpretation, as shown in the panels of 2(b). As expected, the response exhibits a strong component at the excitation frequency Ω = ω f . Moreover, triadic interactions are readily seen in this visualization as side-bands to the forcing frequency component, i.e., Ω = ω f ± ω. All other relevant components of the power spectrum response for this case are explained as higher harmonics of the forcing and the triadic interactions with these higher harmonics. Fig. 2(c) shows the result of our low order model, which we develop in the next section.
For antisymmetric forcing, the power spectrum response is quite different, as shown in Fig. 3(a). In this case, the dominant frequency components in the power spectrum response correspond to the vortex shedding rate Ω = ω and to harmonics of the forcing Ω = mω f , with m being an integer, as shown in Fig. 3(b). The absence of components mω f ± ω induced by nonlinear interactions is expected because the forcing mode and the pair of conjugate response modes have spatial symmetries that are incompatible with the triadic consistency condition, i.e. their wavenumbers do not sum to zero. The other present frequency components, correspond to interactions between the forcing and higher harmonics of the response, which are orders of magnitude weaker in power and do not play a significant role in the wake dynamics.

IV. MODIFIED STUART-LANDAU MODEL
We now summarize the derivation of our mechanistic model, its use to analyze the dynamics of the system, and how to obtain its coefficients directly from data.

A. Derivation of the Model
We derive a parsimonious model for the evolution of the amplitude of the global vortex shedding mode in the turbulent wake behind a bluff body. We begin with the forced Stuart-Landau model in Eq. (2) obtained from weakly nonlinear theory [7]. Our goal is to obtain an interpretable model for the excitation function g(A, F ) that explains the interactions with the periodic forcing F (t) = F (t + 2π/ω f ). We posit a simple structure for g and proceed by systematically increasing its complexity until the resulting power spectrum response exhibits the same features observed in Fig. 2(b). For easier comparison we evaluate each model in this hierarchy at the same actuation frequencies as our experimental data using the same number of sample records with the same length and sampling rate. Measurement noise is added to the numerical solutions before computing the power spectrum response to represent background turbulence [24]. Zero- mean white Gaussian noise is used with a standard deviation proportional to the corresponding response amplitude for each ω f . This enables unbiased comparison between the model response and experimental data, as in Fig. 2(c) and Fig. 2(a), respectively.
We first show the power spectrum response in the absence of forcing (g = 0). As expected, the spectral power is localized at the natural shedding frequency, as shown in the first panel of Fig. 2(c). Next, we consider the case of harmonic forcing at the known actuation frequency, F = εF e iω f t , where ε is the forcing amplitude andF is a constant that depends on the spatial structure being excited in the flow. This harmonic forcing acts linearly on the dynamics, i.e. g = F , as shown in the second panel in Fig. 2(c). As a third model in this hierarchy, we add terms involving quadratic interactions of the forcing mode with the global mode and with its complex conjugate, g = AF o +F 1 +A * F 2 , as shown in the third panel in Fig. 2(c), where we allow for the possibility that different flow structuresF i are excited by each term. Finally, we allow for non-harmonic external forcing, which can be expanded as a Fourier series whereF im depends on the spatial structure being excited in the flow by the m th harmonic of the forcing. The resulting model for the forced complex amplitude of the global mode is a modified Stuart-Landau equation: The coefficientsF 0m ,F 1m andF 2m , for integer m, parametrize the direct and parametric excitation terms induced by the periodic forcing. These differ depending on the forcing configuration, symmetric or antisymmetric, and are identified from data based on the analysis presented in the following section. The model power spectral response remarkably resembles the results obtained from the turbulent flow experiments, as shown in the last panel in Fig. 2(c). When carrying out the analogous procedure for antisymmetric forcing, we find that the same model generalizes to both forcing modes. Nevertheless, in the latter case the termsF 0m andF 2m are negligible due to the absence of triadic interactions, as shown in Fig. 3(c). In the next section, we will show that under nearly resonant forcing, our extended Stuart-Landau model simplifies to the model of Rigas et al. [27], generalized to non-harmonic excitation. In addition, our model also explains the response observed away from resonances in our extensive experimental data set.

B. Model Based Analysis
We now demonstrate the value of our modified Stuart-Landau model. We exploit the fact that the dynamical system given by Eq. (5) is amenable to classical analysis of nonlinear oscillators. When there is no actuation, the solution to Eq. (5) for A = re iθ exhibits an unstable fixed point at r = 0 and a stable limit cycle of radius r 0 = σ r /l r with frequency ω 0 = σ i − l i r 2 0 ; the subscripts r and i denote the real and imaginary parts of the complex constants. In the presence of periodic excitation, we are interested in finding frequency locked solutions, where the phase of the complex amplitude rotates with a mean angular frequency θ = ω, which is a rational multiple of the forcing frequency ω f . For this purpose, we introduce the change of variables A = re iθ = re i(φ+ωt) =Ãe iωt , whereÃ = re iφ has a slowly varying amplitude r and phase φ. Substituting into Eq. (5) and rearranging yields which describes the evolution of the slow complex am-plitudeÃ. Up to this point, no approximations have been made. We have managed to turn the search for frequency locked solutions of A into the search for fixed points ofÃ. Nevertheless, Eq. (6) is non-autonomous, meaning that the dynamics depend explicitly on t. By applying the method of averaging [46], also known as the Krylov-Bogoliubov method, we can further simplify this equation to an autonomous dynamical system. In practice this is achieved by integrating over one period of the complex amplitude T = 2π/ω and neglecting the changes of the slow variables over that time horizon. The terms on the right hand side of Eq. (6) that are kept after the averaging procedure are those that cancel out the fast oscillations ∼ e iωt on the left. Therefore, the resulting expression for the slow dynamics depends on the actuation frequency. The three possible cases are analyzed below: non resonant forcing, harmonic resonance nω f ≈ ω, and subharmonic resonance nω f ≈ 2ω, for integer n. Away from any resonance of the system, nω f = ω and nω f = 2ω, the slow dynamics are governed by the balance between the left hand side of Eq. (6) and the zeroth harmonic of F 0 , as follows Notice thatF 00 changes the eigenvalue of the linear dynamics of the global mode at the origin to σ = σ + εF 00 . As a consequence, the radius and frequency of the stable limit cycle are modified according to r 2 0 = r 2 0 + εF 00 /l r and ω 0 = ω 0 − εF 00 l i /l r .
After identifying σ and l from the unforced dynamics, we can solve forF 00 from Eq. (8) using measurements of the long-term unforced response, r 0 and ω 0 , along with measurements of the modified response away from resonances, r 0 and ω 0 . The values identified forF 00 in the cases of symmetric and antisymmetric forcing of the turbulent D-shaped bluff body wake are shown in Tab. I. In the proximity of an order n harmonic resonance, i.e., nω f = ω, the slow dynamics are governed by the balance between the left hand side of Eq. (6) and the terms on the right hand side that includeF 00 ,F 1n , and F 22n . Furthermore, if the linear interaction of the global mode with the n th harmonic of the forcing dominates over the nonlinear interaction with its 2n th harmonic, thenF 1n F 22nÃ * , resulting in where σ includes the contribution of the interaction with F 0 . Recasting the system in polar form we obtain dr dt = −l r r r 2 − r 2 0 + εF 1n cos(φ), (10a) This system of equations has no explicit time dependence, hence we can search for fixed points by setting dr dt = dφ dt = 0, which represent frequency locked solutions of the form A = re inω f t . Furthermore, equating the r 2 − r 2 0 terms and rearranging yields Thus, a frequency locked solution exists only if the frequency detuning is in the range where r 0 approximates r evaluated at the synchronization boundary. Therefore, Eq. (12) represents the bounds for the n : 1 synchronization region. Furthermore, having already identified l from the unforced dynamics, this expression can be used to identify the coefficientsF 1n from data. This is achieved using experimental measurements of the mean global mode frequency ω, and finding the values of ω f that delimit the corresponding n : 1 frequency lock-on regions, i.e., where |ω − nω f | < δ is satisfied, with δ being a small threshold. Subsequently, we subtract the upper and lower bounds to compute the width of the frequency lock-on regions and substitute them into Eq. (12), allowing the computation of one of these model parameters for every harmonic resonance considered. In the present experiment we observe three harmonic resonances for both symmetric and antisymmetric forcing. The identified coefficients are presented in Tab. I. We now investigate the cases where a subharmonic resonance of the type nω f = 2ω takes place, excluding even values of n as those were accounted for in the previous scenario. In the proximity of a subharmonic resonance, the slow dynamics are governed by the balance between the left hand side of Eq. (6) and the terms on the right hand side that includeF 00 , andF 2n , as follows where σ includes the contribution of the interaction with F 0 . Recasting Eq. (13) in polar form we obtain dr dt = −l r r r 2 − r 2 0 + εF 2n r cos(2φ), As in the previous case, these equations have no explicit time dependence, hence we search for fixed points by setting dr dt = dφ dt = 0, which represent frequency locked solutions that are now of the form A = re inω f /2t . Again, equating the r 2 − r 2 0 terms and rearranging yields for which a solution exists only for the range of the frequency detuning given by This inequality determines the bounds for the n : 2 synchronization region. Moreover, having already identified l from the unforced dynamics, this expression can be used to identify the coefficientsF 2n from data. Similar to the case of harmonic resonance, this is achieved using experimental measurements of the mean global mode frequency ω, and finding the values of ω f that delimit the corresponding n : 2 frequency lock-on regions, i.e., where |ω−nω f /2| < δ is satisfied, with δ being a small threshold. Subsequently, we calculate the width of the frequency lock-on regions and substitute them into Eq. (15), allowing the computation of one of these model parameters for every subharmonic resonance observed. For the present flow, we observe four subharmonic resonances for symmetric forcing, and one for antisymmetric. The identified coefficients are shown in Tab. I.

C. Identifying the Stuart-Landau Coefficients
In the absence of forcing, the proposed model, given by Eq. (5), reduces to the classic Stuart-Landau equation that is parametrized by the complex coefficients σ, and l. Therefore, these unknown parameters can be identified using measurements of transients of the system when the forcing is not active. To characterize the unforced dynamics, we drive the system away from its long-term behavior using steady blowing and using resonant periodic forcing. Starting from each of these conditions, the forcing is turned off and we record the evolution of the global mode. These transients are low-pass filtered using a 5 th -order Butterworth filter with a cutoff frequency of 1.3ω 0 and then phase-averaged over an ensemble of 59 realizations. Figures 4(a) and 4(b) show the evolution of the antisymmetric pressure average starting from larger and smaller amplitude oscillations, respectively. In both cases it takes on the order of 10 shedding cycles for the wake to relax onto the limit cycle describing its long-term behavior.
Based on the data for the unforced transients, we identify σ and l through a constrained least squares regression. The phase-averaged and filtered time-series x(t) and its Hilbert transform x H (t) are used to build an analytic signal for the complex global mode amplitude A(t) = (x + ix H )(t). Once we have the instantaneous complex amplitude, we can compute its instantaneous amplitude r(t) and phase θ(t), as well as the respective time derivatives using 2 nd -order central finite differences. In addition to the transients, we measure the mean amplitude and frequency, r 0 ω 0 , from the magnitude and position of the peak of the PSD of the long-term unforced time-series. In accordance with the Stuart-Landau equation, we know that these measurements must satisfy the solutions for the radius r 0 = σ r /l r and frequency ω 0 = σ i − l i r 2 0 of the stable limit cycle, where the subscripts r and i denote the real and imaginary parts of the complex constants. Therefore, the regression problem is formulated by recasting the Stuart-Landau equation in polar form and deriving linear equality constraints from the known limit cycle solution, as follows Here, the vector b and the matrix Θ are constructed by vertically stacking the time-series data for r,ṙ, andθ and using null vectors of the same length. The matrix C and the vector d are also computed from experimental data, since we have measured r 0 and ω 0 . Therefore, the coefficients ξ can be identified from which is a convex optimization problem that can be solved using readily available software; we use CVXPY developed by Diamond and Boyd [47]. The identified model is simulated and compared against the experimental data in Fig. 4. The identified values for σ and l are shown in Tab. I.

V. RESULTS AND DISCUSSION
In this section we present a comparison between predictions from our proposed model and experimental data, and discuss the synchronization properties of the forced turbulent wake past a D-shaped body. Using our model with the coefficients shown in Tab. I, we compute the long-term amplitude and frequency response of the global vortex shedding mode under periodic forcing as a function of the excitation frequency. A comparison against experimental data for symmetric and antisymmetric forcing at three excitation amplitudes is shown in Figs. 5(a), (b), (d ), and (e). The modified Stuart-Landau model captures the experimental behavior very well, displaying multiple resonances and frequency lock-on regions. Moreover, the same model structure allows for the description of the long-term response with various forcing amplitudes and different forcing configuration, which translates into different dominating resonances; harmonic or subharmonic.
In accordance with the work of Barros et al. [26], when the forcing is symmetric, the largest amplification of the response occurs when the system is excited near twice the fundamental frequency ω f /ω 0 ≈ 2, as shown in Fig. 5(a). As pointed out by Rigas et al. [27], this subharmonic resonance arises due to triadic interactions between the global mode, its complex conjugate and the forcing mode. This can only occur when the spatial characteristics of these modes satisfy the triadic consistency conditions, meaning that, if expressed using a travelling wave ansatz, their wavenumbers and frequencies must add up to zero [ 28,29]. This requirement is clearly not satisfied for the case of antisymmetric forcing, where the largest amplification of the response occurs at the harmonic resonance ω f /ω 0 ≈ 1, as shown in Fig. 5(d ).
We now investigate which combinations of the actuation frequency and amplitude cause the turbulent wake to synchronize with the external excitation. These regions in ω f -ε space are known as Arnold tongues [16]. From experimental data, the mean frequency of the response ω is used to find the excitation frequencies that delimit the n : 1 harmonic and n : 2 subharmonic lock-on regions. These are marked with blue crosses in Figs. 5(c) and (f ) for both symmetric and antisymmetric forcing, respectively. The model based Arnold tongues which are shaded in Figs. 5(c) and (f ) correspond to the regions bounded by Eq. (12) and Eq. (15). As the figure shows, our modified Stuart-Landau model captures the observed tongues and predicts the position of the transition boundaries for the range of excitation amplitudes studied. In this particular example, significant drag reduction is observed outside the Arnold tongues, which highlights the relevance of modeling synchronization as a key enabler of effective periodic flow control.

VI. CONCLUSIONS
In this work, we have leveraged a uniquely comprehensive experimental dataset of a periodically forced turbulent wake to develop a modified Stuart-Landau model for the response of the global vortex shedding mode. The breadth and quality of our dataset reveals previously unobserved resonances and frequency lock-on regions. Moreover, it enables the construction of heat maps showing the power spectral density of the wake response as a function of the excitation frequency. This novel visualization exposes triadic interactions as dominant frequency components, providing the first confirmation of their previously conjectured role in the forced wake response [26,27].
We derive a parsimonious model for the evolution of the forced global mode that generalizes previous models to non-harmonic excitation and extends their applicability to a broader range of parameters. From this model we gain new insights into the nonlinear mechanisms underlying synchronization of the flow response to external excitation. We find excellent agreement between the model and experiments, building a foundation for future investigations into effective periodic flow control. Low-order models that capture synchronization mechanisms are essential for the design of feedback controllers to manipulate periodic coherent structures that govern lift, drag and mixing in oscillator flows. Another key advantage of the model is that its unknown coefficients can be identified directly from data. Our approach to expose the synchronization properties of the system in an extensive parameter space using only a few measurements is promising, as it can be generalized to a large class of forced oscillator flows.