Steady-state statistics, emergent patterns and intermittent energy transfer in a ring of oscillators.

Networks of coupled nonlinear oscillators model a broad class of physical, chemical and biological systems. Understanding emergent patterns in such networks is an ongoing effort with profound implications for different fields. In this work, we analytically and numerically study a symmetric ring of N coupled self-oscillators of van der Pol type under external stochastic forcing. The system is proposed as a model of the thermo- and aeroacoustic interactions of sound fields in rigid enclosures with compact source regions in a can-annular combustor. The oscillators are connected via linear resistive coupling with nonlinear saturation. After transforming the system to amplitude-phase coordinates, deterministic and stochastic averaging is performed to eliminate the fast oscillating terms. By projecting the potential of the slow-flow dynamics onto the phase-locked quasi-limit cycle solutions, we obtain a compact, low-order description of the (de-)synchronization transition for an arbitrary number of oscillators. The stationary probability density function of the state variables is derived from the Fokker-Planck equation, studied for varying parameter values and compared to time series simulations. We leverage our analysis to offer explanations for the intermittent energy transfer between Bloch waves observed in acoustic pressure spectrograms observed of real-world gas turbines.

Introduction

Thermoacoustic instabilities in can-annular combustors

In a rigid enclosed volume, the heat release rate response of an compact unsteady flame to acoustic perturbations forms a feedback loop with the acoustics of the enclosure. If, at a given condition, this interaction between the sound field and the flame exceeds the dissipation due to radiation losses and viscous effects in the fluid, occurring mainly at the boundary of the enclosure, it can give rise to so-called thermoacoustic instabilities. The study of this physical phenomenon goes back to the work of Rayleigh [1]. When an instability is insufficiently damped, it can lead to high-amplitude acoustic pressure oscillations in the enclosure. In industrial machines, such as stationary gas turbines, the resulting pulsations induce high-cycle fatigue cracks in the metal parts surrounding the enclosure, which cause down-time, incurring fees for the manufacturer, and can sometimes lead to catastrophic system failures if a broken part flies through the combustor and collides at high speed with the vanes and blades of the turbine. A review of high-cycle fatigue in gas turbines is given in Ref. [2]. Applied studies on Helmholtz dampers to suppress thermoacoustic instabilities in gas turbine combustors are found in Refs. [3, 4]. Suppression of a thermoacoustic instability in two coupled Rijke tubes in the presence of noise is discussed in Ref. [5]. Sensitivity and nonlinear dynamics of thermoacoustic oscillations are reviewed in [6]. A complex systems approach to modeling and mitigating thermoacoustic instability in turbulent combustors is presented in [7]. Synchronization phenomena and nonlinear energy pumping in two coupled thermoacoustic oscillators are studied experimentally in [8].

Since the turn of the millennium, the modeling, prediction and suppression of thermoacoustic instabilities have drawn renewed interest due to internationally imposed emission limits on power generation systems and the resulting increased demand for lean-premixed combustion. In a lean-premixed combustor, the entire air flow passes through the burner to ensure a sufficiently lean mixture and to avoid rich reactant pockets, which produce significantly more pollutants. There are two consequences to this concept: First, lean premixed flames are much more sensitive to flow perturbations than non-premixed flames, which increases critically their response to acoustic forcing. Secondly, in a premixed setting, one does not need to dilute the combustion products to temperatures that are acceptable for the available cooling of the turbine parts, and hence there are no dilution holes, which provide strong acoustic damping in non-premixed combustors. These two factors lead to combustion instability problems in lean-premixed gas turbines similar to those encountered in rocket engines. In both applications, these instabilities remain a key problem in the design process [9, 10].

Traditionally, most literature on modeling thermoacoustic instabilities focused on silo-type (one single can) or annular combustors (see, for instance, [11-14] and [15-17], respectively). However, as energy demand continues to rise and energy networks are changing, larger gas turbines are required. For these types of gas turbines, which produce more than 750 MW in combined cycles at more than efficiency, and which can be supplied with sustainably produced (hydrogen) that is blended with natural gas, the silo-type and annular combustor designs become uneconomical due to several engineering trade-offs. Therefore, modern high-efficiency H-class gas turbines exclusively feature can-annular combustor architectures. In this type of system, combustion takes place in a number (typically 12 or 16) of cans evenly distributed along the turbine annulus. The cans are thermodynamically decoupled from each other, but the annular turbine inlet is shared by all can outlets and allows for acoustic cross-talk between neighboring cans. Theoretical, numerical and experimental studies performed at Siemens [18-21], General Electric [22-24] and Ansaldo Energia Switzerland [25, 26] are testament to the increased interest of industry in the thermoacoustics of can-annular combustors. In Ref. [25], so-called Bloch modes, or rotating waves of the acoustic pressure field along the turbine annulus, are discussed. In Fig. 8 therein, the authors provide indirect experimental evidence of Bloch modes in a real-world gas turbine by decomposing measured acoustic pressure signals from different cans into rotating wave components of different azimuthal order using the discrete Fourier transform. Direct evidence of Bloch modes occurring in a four-can system, showing the synchronized state and a wave-like phase pattern along the annulus, is shown in Figs. 5 and 6 of Ref. [23], respectively. We note that the oscillators in the latter reference are coupled in a different way (by long tubes) than in the former reference (small gap near the turbine inlet), and that, as discussed in Sect. 1.4, our study focuses on a scenario closer to the setting of the former reference, which is typical for industrial can-annular gas turbines. The theoretical studies in Refs. [27, 28] use periodic (Bloch-type) boundary conditions to simplify the analysis of reduced-order models of can-annular combustors. The effect of resistive (damper-like) and reactive (mass- or spring-like) coupling on the linear stability of a ring of thermoacoustic oscillators is investigated in Ref. [29]. We also mention that nonlinear phenomena in can-annular configurations, such as amplitude death and quenching, have drawn increased interest in recent years [30-33].

Fig. 8

Analysis of a time trace of the slow-flow system given by Eqs. (5) and (6) with () over 5000 cycles. a Time traces of the amplitudes . The oscillators are indicated by blue, red and cyan color, respectively. The limit cycle amplitude , given by Eq. (27), is shown in green. b Time traces of the phase difference . Time traces of other oscillators appear similar. c Phase distribution along the ring at the end time. d Spectrogram of the reconstructed fast oscillating signal , showing the estimated signal power over time t as a function of the normalized frequency . The average of the signal power over the shown frequency and time domains, denoted by , is also indicated in the inset. e Comparison of the analytical PDF of the phase differences (red) with a histogram of (blue) from the time trace shown in (b). f Histograms of and , represented using the geometric analogy described in Fig. 7. Spectro-and histograms of other oscillators appear similar

Aeroacoustic coupling between the cans

The cans are acoustically coupled through the apertures at the turbine inlet, leading to acoustic-hydrodynamic effects which are strongly influenced by the mean axial bulk flow speed of the combustion products : When a sound wave with frequency f from inside the can is reflected at the coupling aperture of streamwise width W, the acoustic (compressible) velocity field inside the can interacts with the hydrodynamic (incompressible) fluid motion in the aperture. As a consequence of this interaction, depending on the value of the nondimensional Strouhal number , either destructive or constructive interference can occur.

On the hydrodynamic side, the acoustic pressure fluctuations can lead to shear layer flapping in the turbulent wake between the cans (see Fig. 1c), which causes vorticity fluctuations that generate sound according to Howe’s acoustic energy corollary [34]. For a study on sound production by shear layer flapping in a laboratory-scale aeroacoustic system (a whistling beer bottle under grazing flow over the bottle top), the reader is referred to Ref. [35]. In extreme cases, at high pressure amplitudes, the shear layer can roll up and form discrete vortices which are shed periodically from the upstream edge of the aperture and produce sound when they meet the downstream edge. Numerical studies of sound generation by vortex shedding in side-branch apertures can be found in Refs. [36, 37].

Fig. 1

a Sketch of the modeled can-annular combustor, from a radial view. The dimensions are not true to scale. Shown are the control volume , the acoustic pressure at the can outlet , the aperture width W, cross-section area and depth d, respectively, the can length L and cross-section area A, respectively, as well as the coherent heat release fluctuations of the flame , which characterize the flame response to acoustic perturbations. b Abstraction of the system illustrated in (a), viewed from an axial perspective, for cans. Shown are the amplitude and phase of the acoustic pressure at the can outlet , respectively, the linear resistive coupling and the reactive coupling . For simplicity, reactive coupling is neglected in this study. c Sketch of the turbulent wake in the aperture between neighboring cans. Shown is a typical profile of the mean axial velocity, the mean axial bulk flow speed of the combustion products and the turbulent wake, which is bounded by two shear layers. The acoustic-hydrodynamic interaction in the aperture leads to the resistive and reactive coupling between the oscillators shown in (b)

The main consequence of the acoustic-hydrodynamic interaction described above is the following. For fixed aperture width W and frequency f, depending on the mean flow speed , the sound field may lose or gain acoustic energy from the interaction with the grazing mean flow over the aperture. Experimental and numerical studies highlight also the role of the nonlinear saturation of the aeroacoustic response of a shear layer [38, 39]: When the acoustic forcing amplitude is increased, the response of the shear layer, in terms of its kinetic energy or its acoustic energy production at the fundamental frequency of excitation, decreases. A theoretical model of this nonlinear shear layer response, based on Howe’s formulation, was derived and validated against experiments in [40]. The same saturation mechanism leads to self-excited aeroacoustic instabilities, which occur, for example, when we whistle. The first modern study of this type of instability was conducted by Sondhauss in the nineteenth century [41]. Readers interested in the role of aeroacoustic instabilities in whistling and musical instruments are referred to the classic experiments of Wilson [42] as well as the more recent review of Fabre et al. [43]. Aeroacoustic instabilities can also occur in industrial machines, where they cause noise pollution and fatigue damage of components [44]. A discussion of the aeroacoustic instability mechanism from an industrial perspective, with a focus on mitigation measures and design rules, is found in [45].

Literature review

The study of synchronization and rotating waves in rings of limit cycle oscillators is not novel. We give below a review of related analyses and discuss how our work differs from these previous efforts. A pair of coupled oscillators is a limiting case of a ring (see, for instance, [46-50]). For compactness, we do not include these works in the review below.

We find the first modern study of rings of coupled van der Pol oscillators in the work of Endo and Mori from the 1970s [51], wherein the authors classify the limit cycle solutions in a symmetric ring of a large number of oscillators with resistive coupling. Deterministic averaging is performed to identify limit cycle amplitudes. All limit cycle solutions are found to be linearly stable, and this finding is confirmed by experiments on an electrical circuit with four elements. Their study is continued in [52], where the effect of delay in the coupling term is studied. It is found that the addition of such a delay can lead to instabilities of certain modes, which is confirmed experimentally. A group-theoretical approach to the classification of phase-locked solutions in rings of coupled oscillators is taken in [53].

In Ref. [54], averaging, linear stability analysis, simulations and experiments are used to study the effect of differing parameters in members of rings of three and four oscillators. In [55], theory and experiments are combined to study the effect of external forcing and non-uniform parameter distribution on the synchronization of a four-member ring of van der Pol oscillators. Following up on their earlier work in Ref. [56], in [57], the authors study the effect of a non-uniform distribution of the reactive coupling parameter along the ring on the system stability by transforming the system into a Hill equation. They show that a harmonic variation of the reactive coupling parameter along the ring can affect the linear stability of limit cycles. We mention also Ref. [58], where it is shown that particular asymmetric parameter distributions can favor the stability of the perfectly symmetric, synchronized solution in a ring of nonidentical limit cycle oscillators.

An approach based on Floquet analysis of Hill’s equation is taken by [59] to study the stability of synchronized states in a large ring of van der Pol oscillators. Therein, the authors put emphasis on linking their results to intercellular dynamics in biology. Concepts from neuroinformatics are invoked in Ref. [60], focusing on synchronization phenomena in a ring of van der Pol oscillators coupled by time-varying resistors.

In Ref. [61], the authors begin by characterizing limit cycle solutions in a ring of (linearly) resistively coupled limit cycle oscillators. They perform a linear stability analysis and derive a quasi-potential (a candidate function for Lyapunov’s second method of stability [62]) to determine the nonlinear stability of these solutions. After a change to amplitude-phase coordinates, deterministic averaging is applied to eliminate the fast oscillating terms. Subsequently, external noise forcing is added to the phase equation to numerically investigate stochastic effects on the synchronization between the oscillators. Using a small-noise approximation, the authors obtain analytical formulas for the mean time of a state to leave an attractor basin. Because the phase noise is added phenomenologically after deterministic averaging is carried out on the original, fast system, it is unclear how Eqs. (11) and (12) relate to Eq. (1) in Ref. [61].

We note that instabilities of azimuthal waves in a discrete, ring-shaped fluid-dynamical system are also encountered, for instance, in Refs. [63-66], where theoretical and experimental methods are combined to study of the dynamics in a ring of bouncing droplets. Amplitude death and phase-flip bifurcations between in-phase and anti-phase synchronization in a system of coupled chemical oscillators are investigated in [67].

Aim of this work

In this work, we combine analytical and numerical tools to study a low-order model of a can-annular combustor. Our goal is to gain a qualitative understanding of some of the system’s dynamic features that are independent of the number of cans N. In particular, we are interested in how, in the steady state, the amplitude and phase statistics of the acoustic pressure in the cans are related to the intermittent energy transfer between different Bloch modes observed in real-world gas turbines. To this end, we propose a symmetric model where all reactive effects of the flame response and the aeroacoustic coupling are neglected, and only resistive effects are taken into account. We show that this simple model is able to describe a wide variety of possible emergent patterns, including synchronization, rotating waves as well as quasi-steady superpositions of clockwise (CW) and counterclockwise (CCW) waves, and that it reproduces the intermittent energy transfer between different Bloch mode components of the acoustic pressure observed in real-world gas turbines. These patterns are characterized and differ from each other by the change of the phase of the acoustic pressure along the ring. For example, a wave-like pattern is described by a sinusoidal variation of the phase and the in-phase synchronized state by a uniform phase distribution. A more detailed discussion is given in the next section after the governing equations of the system have been derived.

A sketch of the physical system we consider is shown in Fig. 1a from a radial view. Shown are the control volume , the acoustic pressure at the can outlet , the aperture width W, cross-section area and depth d, respectively, the can length L and cross section area A, respectively, as well as the coherent (at the acoustic frequency) heat release fluctuations of the flame , which characterize the flame response to acoustic perturbations. Figure 1b shows an abstraction of the same system for cans from an axial perspective, the amplitude and phase of the acoustic pressure at the can outlet , the linear resistive coupling and the reactive coupling . For simplicity, reactive coupling is neglected in this study. We use the convention that a positive increment in j implies a CW rotation by an angle of along the ring, and that corresponds to the oscillator at the 12/N o’clock position. The turbulent wake in the aperture between neighboring cans is sketched in Fig. 1c. Shown is a typical profile of the mean axial velocity in the aperture, the mean axial bulk flow speed of the combustion products and the turbulent wake, which is bounded by two shear layers.

From an acoustic perspective, we restrict our analysis to the case of nearly closed, weakly coupled cavities whose area ratio is smaller than the ratio of effective aperture depth to can length , respectively: . In this regime, the effects of the coupling can be studied as small perturbations of the dynamics of N decoupled cavities. Elementary first-principles calculations quantifying this limit case are provided in [29]. As discussed in this reference, anti-phase synchronization can occur also in a different setting when the cavities are strongly coupled. Therefore, when a phase pattern is referred to, for example, as the “anti-phase synchronized state,” this should not be understood in an exclusive sense, but in the context of the parameter range considered in this work. In other words, the presently considered coupling scenario is not the only situation in which periodic phase patterns along the ring can arise.

In our model, the flame drives a single natural (longitudinal) eigenmode of the can with corresponding eigenfrequency . In a first approximation, we assume that the mode shape of is unperturbed by the thermo- and aeroacoustic interactions and that the acoustic pressure signal is close to harmonic. These are often reasonable assumptions in practice [68, 69]. Following the low-order modeling approach presented in the latter two references as well as Ref. [70], the thermoacoustic instability in each can is described by a van der Pol oscillator with linear growth rate and nonlinear saturation . For simplicity, all oscillators are assumed to be identical. For , where and f is the frequency of oscillation, the acoustic pressure field in the jth can is approximated by the unimodal projectionwhere is the dominant modal amplitude. It is the dynamic variable of the van der Pol oscillator. Following [71], we assume a pressure antinode at the turbine inlet (see Fig. 1a). In this case, is equivalent to the acoustic pressure at the can outlet.

The oscillators are driven by white noise, which represents broadband, quasi-random acoustic perturbations induced by turbulent fluctuations of the mean axial flow in the cans.

The thermoacoustic growth rate can be further decomposed into a sum of the acoustic gains from the flame response and dissipative losses at the boundary of the enclosure, each of which may be measured separately [72]. The parameter describes the saturation of the coherent response of the flame to acoustic perturbations: for , a small acoustic perturbation is amplified by the flame and grows, until the response saturates at a point where the acoustic gain from the flame equals the dissipative losses, and a limit cycle is reached.

The acoustic coupling is incorporated in the oscillator model by linear and van der Pol-type nonlinear resistive terms. The linear terms are analogous to damper and resistor elements in mechanical and electrical oscillators, respectively.

This linear resistive coupling describes the small-amplitude aeroacoustic response of the turbulent wake in the aperture to pressure differences between neighboring cans. Consistent with the low-order model of an aeroacoustic instability validated against experiments in Ref. [38], we include the nonlinear saturation of the aeroacoustic response with increasing pressure amplitude by adding a quadratic term in the pressure difference to the linear coupling . When the nonlinear coupling term is neglected, the governing equations of our model describe a network of nonlinear oscillators with linear coupling. Such systems are ubiquitous in classical mechanics [73-76], chemistry [77-79] and biology [80-82].

We go further than previous studies by accounting for nonlinearity in the (purely resistive) coupling, in combination with additive stochastic forcing of the fast variables, which is consistently included in the averaging procedure. To the authors’ knowledge, this case, which leads to novel and unexpected results, has not been investigated so far.

Due to the symmetry of the system considered in this work, we find that different emergent patterns (synchronized state, rotating waves, etc.) occur at similar amplitudes and are mainly distinguished by different phase dynamics. In regard to this, we mention Adler’s equation for phase injection locking [83-85] and Kuromoto’s classic model [86], which has become a paradigmatic example of phase synchronization in coupled oscillators [87-89].

Overview

In Sect. 2, we present the model and derive the potential governing the amplitude and phase dynamics. In Sect. 3, numerical experiments on the deterministic, noise-free system are performed and analyzed. In Sect. 4, the steady-state statistics of the noise-driven system are derived and studied. Selected results and topics for future research are discussed in Sect. 5. In Sect. 6, we summarize our conclusions.

Theoretical model

In an isolated can, the interaction of a single, dominant acoustic mode with the flame can be represented by a self-excited (van der Pol) oscillator [68-70]. Such low-order oscillator models of thermoacoustic instabilities in silo-type combustors are well-understood in terms of their accuracy compared to more detailed modeling approaches and experimental results (see, for instance, Fig. 17 in Ref. [90]). As shown in the latter reference, the linear growth rate and the nonlinear saturation can be adjusted to obtain a qualitative approximation of the experimentally measured phase space spanned by the acoustic pressure p and its time derivative . Our model combines N of such oscillators to approximate the acoustic pressure dynamics in a can-annular combustor.

Following Ref. [38], we account for the aeroacoustic interaction between neighboring cans with the linear resistive coupling and its nonlinear saturation .

Stochastic forcing of the acoustic field by broadband noise production in the turbulent flow in the cans is modeled by additive, zero-mean white noise signals of equal intensity with the autocorrelation functions , where denotes the time integral over the real line and is the Dirac delta function.

Based on the above discussion on low-order models of thermo- and aeroacoustic instabilities, we propose the following model of N coupled van der Pol oscillators for the acoustic pressure dynamics in the can-annular combustor:for and , . The sign before is positive because the coupling is amplifying if and we assume that as in Ref. [40], the constructive aeroacoustic coupling saturates at high amplitudes (monotonic decay of constructive feedback strength). For simplicity, this work focuses on situations where the nonlinear aeroacoustic coupling can be considered as a small perturbation of the thermoacoustic oscillations in the cans, so that the saturation associated with the thermoacoustic instability is significantly larger than the saturation of the aeroacoustic coupling: . We exclude the case of vanishing resistive coupling in our study, because a vanishing linear response with nonzero saturation is not physical. The parameter values used in the numerical examples throughout this work are listed in Table 1. If multiple values are used, the default value is shown in bold. It is indicated below if values different from the default ones are used.

Table 1

Parameter values used in the numerical examples throughout this work. If multiple values are used, the default value is shown in bold

Parameter	Meaning	Value
ω_0/2π	Frequency of dominant acoustic mode	255 Hz
ν	Thermoacoustic growth rate	{25,50} s^-1
κ	Thermoacoustic saturation	0.48 s^-1 Pa^-2
G=κΓ/8|ν|^2ω02	Normalized white noise intensity	2.88×{10^-1,10^1,10^3}
Λ=2λ/ν	Normalized linear resistive coupling	∈[-1,2]
K=16ϑ/κ	Normalized coupling nonlinearity	{0,1/6,1/3}
N	Number of oscillators	{3,12}
b	Bloch wavenumber	{0,±1,…,±floor(N/2)}
β	Continuous Bloch wavenumber	∈[-N/2,N/2]
B=πβ/N	Normalized Bloch wavenumber	∈[-π/2,π/2]

We define the coordinate change to amplitude-phase variables , , as follows:where , and is the two-argument arctangent function. Assuming that (A) the variation of and over an acoustic period is negligible and (B) the oscillators perform quasi-sinusoidal motion, deterministic and stochastic averaging (see Refs. [91, 92] and [93, 94], respectively) can be performed on Eq. (2) to eliminate the fast-oscillating terms. This averaging is done in Appendix A and yields the slow-flow dynamics in amplitude-phase form:andfor , where is the phase difference and and are 2N uncorrelated white noise sources with the same noise intensity , respectively. The above definitions of and were used in the derivation of the slow-flow system, given by Eqs. (5) and (6). For the presentation of our results in Sects. 3 and 4, we set and . These operations, which leave Eqs. (5) and (6) unchanged, are equivalent to adding integer multiples of to and so that the sums lie in the domains and , respectively.

The slow flow system (5) and (6) is the main focus of this work. It has the form of a 2N-dimensional Langevin equation [95]:where and the potential is defined as follows:Equations (7) and (8) imply that trajectories , , are stochastically perturbed and attracted toward lower values of the potential defined in Eq. (9).

In the deterministic case, trajectories of the fast and averaged systems are pointwise close to each other over time, and the degree of closeness is determined by how strongly the averaging assumptions (A) and (B) are satisfied [91, 92]. For a linear damped harmonic oscillator, the averaged system is equivalent to the fast system. In contrast to averaging in deterministic systems, due to the lack of smoothness, we do not expect trajectories of the stochastically averaged system to be pointwise close to those of the fast system, but for the joint probability density function (PDF) , , computed from the fast and the averaged system, given by Eq. (2) and Eqs. (5) and (6), respectively, to converge in the steady state, i.e., for .

Several theorems on stochastic averaging exist [93, 94]. It is beyond the scope of this work to elaborate on these mathematically intricate results, which would require introducing additional notation and theory to understand under which conditions they apply to the fast system (2). Instead, it is accepted that stochastic averaging is justified by a combination of physical and mathematical arguments, and careful numerical validation is performed to ensure that, for the parameter range considered, the joint PDFs of the fast system (2) and the slow amplitude-phase system (5) and (6) indeed coincide over long time spans , , where is the relaxation time ( time is takes to reach the steady state) and is the correlation time of the signal “a”. The correlation time is defined as , where denotes the one-sided limit from above (see Ref. [93], p. 65, Eq. (3.61) with ). This validation ensures that any qualitative statistical observations from long time series of the averaged system carry over to the fast system.

We note that in real-world experiments and in numerical simulations, the signals , and will generally have nonzero correlation times. However, because we model , and as white noise signals with (theoretically) vanishing correlation times, we assume throughout this work that , i.e., that the time to reach the steady state significantly exceeds the correlation time of the external stochastic forcing.

Potential landscape of the deterministic system

In the absence of noise, . Under the averaging assumptions (A) and (B), the limit cycles of the deterministic part of the fast system (Eq. 2 with ) are given by fixed points of the slow-flow system (5) and (6) with . We find that these fixed points describe solutions with uniform amplitude along the ring and equal phase difference between neighboring oscillators:The periodicity condition imposes the following additional restriction on the phase difference:We call the solutions given by Eqs. (10)–(12) Bloch modes and b the Bloch wavenumber because of the present result’s similarity to Felix Bloch’s classic theory of electron dynamics in perfectly periodic crystal lattices [96]. For compactness, we also introduce the continuous and normalized Bloch wavenumbers and , respectively.

Solutions with appear as rotating waves along the ring [61]. We show below that under the current assumptions, these solutions are pairwise degenerate, i.e., that the linearization around the solutions with positive and negative b have the same exact eigenvalues. At the limit cycle, and , so that, for , under a negative increment in j (CCW shift along the ring) and a suitably chosen positive increment in t, the acoustic pressure remains constant. This implies that positive values of b () correspond to CCW rotating waves and negative values to CW rotating waves. In the absence of explicit asymmetries, both CW and CCW rotating waves have the same amplitude and frequency. The case is the in-phase synchronized solution, for which all oscillators reach maximum and minimum pressure simultaneously. For even N, we call the solution with the anti-phase synchronized solution (also referred to as the “push-pull mode” in thermoacoustics literature), which features a phase difference of between neighboring oscillators, so that neighboring oscillators reach maximum and minimum alternatingly.

We note that whereas the spinning direction the nodal lines of an apparent wave is determined by the sign of b, due to the finite number of members in the oscillator ring N, if , there appear to be fewer nodal lines than expected (a Bloch mode with wavenumber b has exactly |b| nodal lines) and they appear to spin in the opposite direction than expected from the sign of b. This visual phenomenon is a consequence of a spatial analogue of the wagon-wheel effect [97] and is discussed more in Sect. 5.

For vanishing coupling , the system (2) corresponds to N decoupled van der Pol oscillators, whose weak-nonlinearity limit cycle solution , , is asymptotically stable [70]. By the hyperbolicity of this limit cycle and the implicit function theorem [98], linear stability persists for small enough coupling parameters , which act as small perturbations of the decoupled system. Indeed, for small enough , all fixed points are linearly stable. This can be seen by linearizing Eqs. (5) and (6) around the fixed points defined in Eqs. (10)–(12) by writing and , where and are small perturbations of the limit cycle amplitude (10) and the phase difference (12). Applying periodic boundary conditions to the linearization (see Refs. [25, 27, 28]) yieldsThis implies that the eigenvalues of the linearization of the system (5) and (6) around the degenerate rotating waves given by (10) and (12) are . For given b and , the linear stability of the corresponding limit cycles, demonstrated by Eq. (13), shows that all Bloch modes are local minima of the potential given by Eq. (9). Therefore, for small enough linear resistive coupling , trajectories may converge to any of these solutions, dependent on their initial conditions.

For oscillators, the limit cycle solutions corresponding to and are shown in Fig. 2. Shown are the values of the phase difference and the phase distribution along the ring.

Fig. 2

Illustration of the limit cycle solutions corresponding to and for oscillators. Shown are the values of the phase difference and the phase distribution along the ring. A positive value of b corresponds to a CCW rotating wave and a negative value to a CW rotating wave. The case is the synchronized solution, for which all oscillators are in phase

Substituting Eqs. (10)–(12) into the slow-flow system (5) and (6) yields the limit cycle amplitudeswhere we have defined is the normalized linear resistive coupling and the normalized coupling nonlinearity.

Limit cycles for a given value of the Bloch wavenumber b can only exist ifotherwise the amplitude defined in Eq. (15) does not exist. Indeed, as evidenced by Eq. (13), when exceeds this value for a given b, the corresponding rotating waves are linearly unstable. This occurs when the acoustic energy losses from the response to small amplitude perturbations of the turbulent wake in the apertures between the cans exceed the gain of acoustic energy from the coherent response of the flame.

In the following, we frequently make use of the assumption that all oscillators converge to quasi-limit cycles so that the dynamics take place in the vicinity of a two-dimensional submanifold of the 2N-dimensional phase space , wheredefined by , where . These quasi-limit cycle solutions include the case of true limit cycles, which are fixed points of the slow-flow system (5) and (6), for . However, we keep k general because, as we show below, quasi-steady solutions with uniform amplitudes and uniform absolute phase differences play a crucial role in the nonlinear dynamics of the oscillator ring. The projection of the potential defined by Eq. (9) onto the submanifold readsNormalizing this result, and using the fact that for a thermoacoustic instability, , yieldswhere is the normalized projection of the potential onto the phase-locked quasi-limit cycle solutions. Note is independent of the signs in k and is symmetric with respect to positive and negative B. We show below that provides a compact description of the (de-)synchronization transition in the parameter space which is independent of the number of oscillators.

The Taylor expansion of the projected potential , given by Eq. (21), around (synchronized state) isand the expansion around isFor (dissipative coupling) and , we consider Eq. (22), which shows that the global minimum of the potential is at , i.e., the synchronized state. Conversely, for (strong amplifying coupling) and , we consider Eq. (23) which shows the global minima of are at (or for odd N), corresponding to anti-phase synchronization.

For , which we call the “transition region” in the following, the minima of the projected potential , given by Eq. (21), lie between and . This parameter range is further investigated below. By computing the zero level set of the slope of with respect to B, we obtain the curve parameterizing the minima of in the transition region as a function of :Numerical examples were performed and analyzed to explore the projected potential , given by Eq. (21), by launching a number of realizations of the slow-flow system (5) and (6) from random initial conditions. A part of the results is reported in Fig. 3, where, for oscillators and , realizations were computed for , respectively, and realizations were computed for oscillators, and . Trajectories of the oscillators and 3 are shown in blue, red and cyan, respectively. For , only the first three oscillators are shown, but trajectories of other oscillators appear similar. In Fig. 3a–c, the distribution of the projected potential as a function of and is shown. These insets illustrate how a nonzero coupling nonlinearity leads to a transition region for , bounded by the dotted lines, which separates the domains where the in-phase (In-phase) and the anti-phase synchronized state (Anti-phase) are the global minimum of the potential , given by Eq. (9), respectively. The solid black curve indicates the location of the minimum of with respect to 2B in the transition region. The blue, black and red crosses in 3c mark the location of the minimum of the projected potential , given by Eq. (21), for , and , respectively. The dashed black curve demarcates the domain within which Bloch modes are linearly unstable (LUW: linearly unstable waves). Figure 3d, g and j shows the variation of the projected potential as a function of 2|B|. In Fig. 3e, f, (), h, i (), it is shown how, in a ring of oscillators with dissipative coupling, all trajectories converge to the synchronized state and for strong amplifying coupling with , all trajectories converge to rotating wave solutions, with a symmetric distribution between CW () and CCW () rotating waves. The number of realizations that converge to either solution class is indicated above the -curves. As shown in Fig. 3k, l, for and , out of all trajectories that converge to solutions with positive b (97 out of 200 total realizations), the majority (52 out of 97) converge to the limit cycle with , but a significant amount of them do not, and instead converge to the nearby solutions with (43 out of 97) and (2 out of 97). The small inset in Fig. 3k shows a comparison between a single realization of (black) computed from the fast system (2) and the amplitude (blue) computed from the slow-flow system (5) and (6) for 5 cycles, starting at 200. Figure 3m shows an enlarged version of Fig. 3l with between 600 and 800 cycles. The corresponding Bloch modes are visualized in Fig. 3n. Similar aggregating behavior around the synchronized state with is observed for , see Fig. 4.

Fig. 3

Numerical examples performed on the slow-flow system (5) and (6) in the deterministic limit with . The distribution of the projected potential , given by Eq. (21), as a function of and is shown in (a)–(c) for increasing values of the coupling nonlinearity K. The solid black curves parametrize the minima of with respect to 2B over . The dotted lines bound the transition region, separating the domains where the in-phase (In-phase) and the anti-phase synchronized state (Anti-phase) are the global minimum of the potential , given by Eq. (9), respectively. Beyond the dashed black line are linearly unstable rotating waves (LUW). The blue, black and red crosses in Fig. (c) mark the location of the minimum of the projected potential , given by Eq. (21), for , and , respectively. The variation of over 2|B| is shown in (d), (g) and (j). Time traces from random initial conditions of the amplitudes and phase differences are shown in (e), (f) (), (h) and (i) () for oscillators (100 realizations) and in (k) and (l) for oscillators with (200 realizations). The oscillators are indicated by blue, red and cyan color, respectively. In the small inset in (k), the numerical solution (black) of the fast system (2) is compared with the amplitude (blue) computed from the the slow-flow system (5) and (6) for five cycles, starting at 200. The numbers above the trajectories in (i) and (l) indicate how many realizations converge to the solution classes with (CCW and CW rotating waves, respectively). The inset (m) shows an enlarged version of the time traces with in (l) between 600 and 800 cycles. In (n), it is indicated how the 97 realizations in (l) with are distributed among different Bloch modes

Fig. 4

Numerical experiments near and in the transition region, which separates the two domains where the in-phase (In-phase) and the anti-phase synchronized state (Anti-phase) are the global minimum of the potential , given by Eq. (9), respectively, for an oscillator ring with members. a Distribution of the projected potential , given by Eq. (21), as a function of and . The colored crosses indicate the location of the minimum of the projected potential , given by Eq. (21), over B for six evenly spaced values of the linear resistive coupling , respectively. The corresponding distributions of the projected potential over the phase difference 2|B| is shown in (b)–(g). The positions of the minima of are indicated in (e) and (f), where they do not correspond to discrete Bloch modes. For each value of , 500 realizations of the slow-flow system (5) and (6) were computed over 5000 cycles from random initial conditions. In (h)–(m), the phase differences , , are shown in blue, red and cyan, respectively. Time traces of other oscillators appear similar. In (h), (i), (j) and (m), the first 20 realizations converging to each Bloch mode are shown and in (k) and (l) the first 50 of all realizations are shown. The average of over all realizations at 5000 cycles, , is displayed in the same insets. In (h)–(m), the smaller insets in the upper right corners show the relative share of the total realizations of that fall into the bins corresponding to the different Bloch modes. Note that while in (h)–(j) and (m), all realizations converge to discrete Bloch modes, in (k) and (l) (transition region), trajectories generally ( and , respectively) converge to unsteady quasi-limit cycles which are superpositions of CW and CCW rotating waves, see Fig. 5. However, Bloch modes are also observed in the transition region; see, for example, the time trace marked with an asterisk () in (l), which corresponds to

Fig. 5

Example of a transition mode in the steady state, computed from the slow-flow system (5) and (6) for (). Shown is a the variation of the amplitudes and b the phase differences for a single realization over 300 cycles in the steady state, c the distribution of the projected potential over (magnified view of Fig. 4f with a different colormap) and d the phases at 30,000 cycles. The colors of the time traces in (a) and (b) correspond to those of the small rings in (d). Time traces of other transition modes appear similar, with more or less pronounced asymmetry between different oscillators, but with absolute phase differences near the minimum of over 2|B| (see Fig. 4k ,l). The phase distribution along the ring shown in (c) corresponds to a standing wave-type pattern of the acoustic pressure along the ring. Animated examples of transition mode time traces are contained in the Supplementary Material

Figure 3 shows that while a change in from 0.5 to leads to O(1) changes in the phase differences , the corresponding change in the amplitudes in the steady state is small compared to . To better understand this, note that on the one hand, for , the synchronized state is globally attracting, and the system aggregates around low-order Bloch modes , which implies a small influence of and K on the limit cycle amplitude , see Eqs. (15)–(17). On the other hand, for and , we have , i.e., and K tend to cancel each other out. For this reason, and because the aeroacoustic coupling is considered as a small perturbation of the thermoacoustic limit cycles (, ), the coupling-induced changes in the amplitude are in general assumed to be small. For , solutions with reach limit cycle amplitudes exceeding that of the synchronized state: . However, such small variations of the oscillation amplitude are difficult to quantify in experiments on real aero- and thermoacoustic systems, which are typically plagued by significant noise levels. Accordingly, the present work mostly focuses on the phase differences , whose distribution quantifies emergent patters of the acoustic pressure along the ring.

More numerical experiments were performed on the slow-flow system (5) and (6) to explore the transition region. The corresponding results are shown in Fig. 4. Figure 4a shows the projected potential , given by Eq. (21), as a function of and . The colored crosses in 4a mark the location of the minimum of the projected potential , given by Eq. (21), over B for six evenly spaced values of the linear resistive coupling , respectively. For the same values of , the insets in Fig. 4b–g show the distributions of over the phase difference 2|B|. The positions of the minima of are indicated in Fig. 4e, f, where they do not correspond to discrete Bloch modes. Five hundred realizations of the slow-flow system (5) and (6) were computed over 5000 cycles from random initial conditions at each value of . The insets in Fig. 4h–m show the evolution of the phase differences of the oscillators with and 3 in blue, red and cyan, respectively. In Fig. 4h–j, m, the first 20 realizations converging to each Bloch mode are shown, and in Fig. 4k ,l, the first 50 of all realizations are shown. The average of over all realizations at 5000 cycles, , is displayed in the same insets. In Fig. 4h–m, the smaller insets in the upper right corner show the relative share of the total realizations of that fall into the bins corresponding to the different Bloch modes. Note that while in Fig. 4h–j, m, all realizations converge to discrete Bloch modes, in Fig. 4k, l (transition region), trajectories generally ( and , respectively) converge to unsteady quasi-limit cycles which we call “transition modes.” These solutions are superpositions of CW and CCW rotating waves, see Fig. 5. Consistent with the linear stability analysis performed above, in the transition region, Bloch modes are also observed; see, for example, the time trace marked with an asterisk () in Fig. 4l, which corresponds to .

We observe in Fig. 4 that outside the transition region, for or , when the global minimum of the potential , given by Eq. (9), coincides with a fixed point of the slow-flow system (5) and (6), all trajectories converge to discrete Bloch modes. For and , nearly symmetric unimodal and bimodal distributions of the phase difference over and are observed in the numerical experiments, respectively.

Inside the transition region, for , we observed in our numerical experiments that trajectories generally do not converge to the phase-locked fixed points of the slow-flow system (5) and (6), but to unsteady quasi-limit cycles, the transition modes, which feature small-scale, piecewise linear periodic variation of the amplitudes and phase differences ’. In the transition region, the mean over all realizations matches well with the minimum of the projected potential over , see Fig. 4e, f, k and I.

An example of a transition mode in the steady state, computed from the slow-flow system (5) and (6) for (), is shown in Fig. 5. Shown in Fig. 5a, b are the variation of the amplitudes and the phase differences for , respectively, for a single realization over 300 cycles in the steady state. Figure 5c, d shows the distribution of the projected potential over (magnified view of Fig. 4f with a different colormap) and the phases along the ring at the end time, respectively. The colors of the time traces in Fig. 5a, b correspond to those of the small rings in Fig. 5d. The phase differences perform slow, small-scale periodic motion over time. Because the transition modes feature both positive and negative phase differences , these solutions are superpositions of CW and CCW rotating waves which appear as standing wave-type patterns of the acoustic pressure along the ring. After a transient, the amplitudes all remain quasi-steady, performing small-scale, periodic motion around . Time traces of other transition modes appear similar, with more or less pronounced asymmetry between different oscillators, but with absolute phase differences near the minimum of over 2|B|, see Fig. 4k, l.

Animated examples of transition mode time traces (, , and ), computed from the slow-flow system (5) and (6), are contained in the SupplementaryMaterial.

The numerical results reported in Figs. 4 and 5 suggest that in the transition region, the occurrence of quasi-limit cycles near the minimum of the projected potential , given by Eq. (21), is generic. Despite these solutions not being fixed points of the slow-flow system (5) and (6), through this unsteady motion, the system reaches lower values of the potential , given by Eq. (9), and hence the transition modes are typically favored over the discrete Bloch modes in the transition region. Further research is required to explore the dynamic nature of these transition modes. Open topics are discussed in Sect. 5.

To better understand the transition modes’ unsteadiness, we consider the continuous Bloch wavenumber and assume the system reaches a quasi-limit cycle where , , and . Of course this solution can never be exactly observed, because these quasi-limit cycles are not fixed points of the deterministic part of the slow-flow system given by Eqs. (5) and (6). The normalized potential is independent of the signs in k, see Eq. (21). However, the signs do impact the projection of the slow-flow equations (5) and (6) onto the quasi-limit cycle solutions :where the fact that in the steady state, was used and , are obtained by replacing b with in Eqs. (16) and (17), respectively. For a distribution of with varying signs in k along the ring, Eq. (26) implies unsteadiness of the Bloch wavenumber . This unsteadiness carries over to the equation for through and . The ad hoc argument made above to explain the transition modes’ unsteadiness is not an exact description of the nonlinear dynamics in the phase space (see Eq. 19), which are responsible for the occurrence of transition modes and require further investigation.

In this section, the potential landscape of the deterministic system was investigated. To this end, the potential of the amplitude-phase system (5) and (6) was projected onto the phase-locked, uniform-amplitude quasi-limit cycle solutions to obtain the normalized projected potential , which compactly describes the (de-)synchronization transition, dividing the parameter space into two domains where the in-phase and the anti-phase synchronized state are the global minimum of the potential , respectively. These domains are separated by a transition region, whose size is proportional to the coupling nonlinearity K and which vanishes for . Notably, is both non-dimensional and independent of the number of oscillators N.

Steady-state statistics of the noise-driven system

Building on our analysis of the projected potential in the previous section, we begin our investigation of the noise-driven system by studying the behavior of under varying coupling nonlinearity K and noise intensity . For , the fixed points of the deterministic part of the slow-flow system (5) and (6) are still characterized by uniform amplitudes along the ring and constant phase differences between neighboring oscillators. The phase differences are still defined by Eq. (12), while the limit cycle amplitudes are affected by the noise intensity :In the case , the trajectories of the slow-flow system are governed by a gradient system perturbed by white noise forcing, described by Eqs. (7) and (8). Similar to the previous section, it is assumed that in the steady state, each oscillator performs stochastically fluctuating motion near the quasi-limit cycles of the deterministic system with uniform amplitude and absolute phase differences . The normalized projected potential depends now also on the normalized noise intensity and is computed for different values of K and G by substituting the limit cycle amplitude (27) and into the potential , given by Eq. (9), and dividing the result by the factor for normalization. For compactness, the resulting expression is only listed implicitly here:where and are obtained by replacing b with in Eqs. (12) and (27), respectively.

The distribution of the projected potential as a function of and is studied in Fig. 6. The dashed black lines mark the transition region . Note that depends now on both K and G. Above and below the transition region, the in-phase and the anti-phase synchronized state are the global minimum of the potential , given by Eq. (9), respectively. The solid black curves in the transition region parametrize the minima of with respect to 2B, defined by , over . In the bottom right inset in Fig. 6, . We observe that depends monotonously on both K and G, and it is therefore justified for qualitative analyses to study the system for fixed values of K and G as we do in this work. If K or G are increased, for , the critical value of the linear resistive coupling , which is computed numerically in the noise-driven case, increases in magnitude, and hence the size the transition region increases. Furthermore, for increasing K or G, the overall potential landscape becomes increasingly biased toward the synchronized solution. The top row of Fig. 6 should be interpreted with care, because at high noise levels, the amplitudes are in general not all simultaneously close to some amplitude and the projection onto the limit cycle solutions loses its validity. Despite this, the figure illustrates correctly the monotonous trends of , given by Eq. (28), with increasing coupling nonlinearity K and noise intensity G. The curve parametrizing the minima of with respect to 2B over is obtained by numerically computing the zero level set of .

Fig. 6

Distribution of the projected potential , given by Eq. (28), as a function of and for varying noise intensity G and coupling nonlinearity K. The dashed black lines mark the transition region . In the bottom right inset, . Above and below the transition region, the in-phase and the anti-phase synchronized state are the global minimum of the potential , given by Eq. (9), respectively. In the transition region, the solid black curve, defined by , parametrizes the minima of with respect to 2B over

For , in the steady state, the probability P of the fast system (2) being within a domain in the 2N-dimensional phase space of the fast system (2), , is defined aswhere the stationary PDF of the fast variablescan be deduced from the stationary Fokker–Planck equation (see Ref. [95]) of the slow-flow system (5) and (6), given bywhere the newly introduced (slow) variablesare related to (, ) by and , whereThe derivation of is detailed in Appendix B. The final result is, where , , the potential is given by Eq. (9) and is a normalization constant such that . In the figures below, the PDFs are normalized by numerical integration over the shown, finite domains.

In line with our analysis of the deterministic system in the previous section, we now project onto the phase-locked, uniform-amplitude quasi-limit cycle solutions to obtain the projected stationary PDF , whereThis result implies that on the submanifold of phase-locked, uniform-amplitude quasi-limit cycles, all oscillators perform uncorrelated stochastically fluctuating motion around the deterministic limit cycles and that in the steady state, the PDF of each of the oscillators, up to normalization, is given by .

Fig. 7

Geometric analogy between emergent patterns in the oscillator ring, quantified by the phase difference , and polar coordinates. Possible nonnegative discrete Bloch wavenumbers are shown for oscillators. The x-axis () corresponds to in-phase and the y-axis () to anti-phase synchronization. Values of b in between the two extremes correspond to CWW rotating waves. The lower half-plane (not shown) is symmetric, with negative values of b and CW waves

The PDF of the phase differences is obtained by evaluating the limit cycle amplitude in Eq. (27) in terms of using Eqs. (16) and (17) with and substituting the result into the potential in Eq. (36). The result is, up to normalization,To graphically compare the joint PDF to normalized histograms of and obtained from time traces, we propose a geometric analogy between emergent patterns in the oscillator ring, quantified by the phase difference , and polar coordinates. The analogy is shown in Fig. 7. Possible nonnegative discrete Bloch wavenumbers are shown for oscillators. The x-axis () corresponds to in-phase and the y-axis () to anti-phase synchronization. Values of b in between the two extremes correspond to CWW rotating waves. The lower half-plane (not shown) is symmetric, with negative values of b and CW waves.

In Fig. 8, we analyze a time trace of the slow-flow system (5) and (6) with (), so that the considered condition lies outside the transition region. The oscillators are indicated by blue, red and cyan colors, respectively. In Fig. 8a, we observe that the amplitudes fluctuate stochastically near the deterministic limit cycle amplitude , which is indicated in green. The phase difference , shown in Fig. 8b performs similar fluctuations and jumps intermittently between . Time traces of other oscillators appear similar. The phase distribution along the ring after 5000 cycles is shown in Fig. 8c. Figure 8d shows the spectrogram of the reconstructed fast oscillating signal , i.e., the estimated signal power , where is the power spectral density of the signal and is the equivalent noise bandwith [99] of a small frequency increment around f (the shown frequency domain is discretized into 257 such increments), over time t as a function of the normalized frequency . We observe that the signal power is mainly concentrated near the natural eigenfrequency . The average of the signal power over the shown frequency and time domains, denoted by , is also indicated in the inset. Figure 8e shows a comparison of the analytical PDF of the phase differences (red) with a histogram of (blue) obtained from the time trace shown in Fig. 8b. As expected, for the considered conditions (strong amplifying coupling), we find symmetric bimodal distributions of the PDFs, with peaks around the anti-phase synchronized state with . Spectro- and histograms of other oscillators appear similar.

Figure 9a depicts the amplitudes of the time trace shown in Fig. 8a over the last five cycles for . In Fig. 9b, the distribution of the amplitudes along the ring at the end time is shown. We observe that the assumption of slowly varying, uniform amplitudes is approximately satisfied.

Fig. 9

a Amplitudes of the time trace shown in Fig. 8 over the last 5 cycles. b Distribution of the amplitudes along the ring after 5000 cycles

The stochastic averaging method is validated against numerical simulations in Fig. 10 by comparing histograms of and , computed over 10,000 cycles from the fast and the averaged system, given by Eq. (2) and Eqs. (5) and (6), respectively, for varying linear resistive coupling and noise intensity G. Histograms of other oscillators appear similar. The results show good qualitative agreement between the two methods. Both methods describe a symmetric, bimodal distribution for . At high noise intensity, the PDF becomes more spread out and less discernible patterns are visible. Note that the amplitude depends on the noise intensity G and changes from the left to the right column. Sharper shapes of the PDFs emerge for longer times, but the computational cost of simulating the averaged system (5) and (6) became prohibitive for longer time traces than those shown in Fig. 10. Nevertheless, the analytical results derived above in this section enable further validation of the stochastic averaging method.

Fig. 10

Validation of the stochastic averaging method against numerical simulations. Histograms of and , computed over 10,000 cycles from the fast and the averaged system, given by Eq. (2) and Eqs. (5) and (6), respectively, are compared for varying linear resistive coupling and noise intensity G. Histograms of other oscillators appear similar. Note that the amplitude depends on the noise intensity G and changes from the left to the right column

In Fig. 11, the stochastic averaging method is validated against analytical results by comparing histograms of and computed from the fast system (2) (Sim.) to the joint PDF , given by Eq. (36) (Theory), for varying linear resistive coupling and noise intensity G. Histograms of other oscillators appear similar. At low noise, time traces were computed for 10,000 cycles and at high noise for 100,000 cycles. For high noise and amplifying coupling, the maxima of the PDF are indicated by black crosses. The results show good qualitative agreement between the two methods. Note that the amplitude depends on the noise intensity G and changes from the left to the right column.

Fig. 11

Validation of the stochastic averaging method against analytical results. Histograms of and , computed over 10,000 cycles at low noise and 100,000 cycles at high noise from the fast system (2) (Sim.) are compared to the joint PDF defined in Eq. (36) (Theory) for varying linear resistive coupling and noise intensity G. For high noise and amplifying coupling, the maxima of the PDF are indicated by black crosses. Histograms of the other oscillators appear similar. Note that the amplitude changes from the left to the right column

For completeness, we show in Fig. 12 transients of the deterministic time traces contained in the Supplementary Material, showing the color-coded phase differences as a function of time and the oscillator number j. For the stochastic time traces, the figure depicts the evolution over the first 1000 cycles. Parameter values are indicated over the respective insets.

Fig. 12

Transients of deterministic time traces contained in the Supplementary Material, showing the color-coded phase differences as a function of time and the oscillator number j. For the stochastic time traces, the figure depicts the evolution over the first 1000 cycles. Parameter values are indicated over the respective insets

The spectrograms of , shown in Fig. 8d for , can be decomposed into spectrograms of different rotating wave components by defining the modal amplitude of the Bloch mode with wavenumber b using the discrete Fourier transform [25]:In Fig. 13, we plot spectrograms of for . Spectrograms for positive and negative b appear similar. Random, intermittent energy transfer between different modes is observed. The average of the signal power over the shown frequency and time domains, denoted by , is also indicated in the insets. As expected, for the given conditions (, ), the anti-phase synchronized state with dominates the power spectrum.

Fig. 13

Spectrograms of the modal amplitudes for , showing the estimated signal power as a function of time t and the normalized frequency . The modal amplitudes are obtained by decomposing the reconstructed fast pressure signals from the example shown in Fig. 8 (time trace of the slow-flow system (5) and (6) over 5000 cycles for ) into different Bloch mode components using the discrete Fourier transform, see Eq. (38). Spectrograms for positive and negative b appear similar. The average of the signal power over the shown frequency and time domains, denoted by , is also indicated in the insets

In Fig. 14, the stationary PDF of the phase differences , given by Eq. (37), is compared to histograms of obtained from time series simulations of the fast system (2) over 15,000 cycles for varying noise intensity G and coupling nonlinearity K. For high noise and amplifying coupling, the maxima of the PDF are indicated by black crosses. The plots show good agreement between the two methods. Histograms obtained from other oscillators or from time traces of the averaged system (5) and (6) appear similar.

Fig. 14

Comparison of the stationary PDF of the phase differences , given by Eq. (37) (Theory), to histograms of obtained from time series simulations of the fast system (2) over 15,000 cycles (Simulation). For high noise and amplifying coupling, the maxima of the PDF are indicated by black crosses. Histograms obtained from other oscillators or from the averaged system (5) and (6) appear similar

The evolution of the joint PDF , given by Eq. (36), is visualized in Fig. 15 for varying (), with (“Low noise” in Fig. 11). The values of the normalized linear resistive coupling correspond to those used in the numerical experiments reported in Fig. 4. In the transition region, the maxima of the PDF are indicated by black crosses. The figure demonstrates how the above analysis, combined with the geometric analogy introduced in Fig. 7, enables a compact description of the (de-)synchronization transition in dependence of , independent of the number of oscillators N. The joint PDF reproduces the transition from a unimodal to a bimodal distribution as is varied from 0.5 to , which was also observed in the numerical experiments on the deterministic system reported in Figs. 3 and 4.

Fig. 15

Visualization of the (de-)synchronization transition. Shown is the evolution of the joint PDF , given by Eq. (36), for varying with (“Low noise” in Fig. 11). The values of the normalized linear resistive coupling correspond to those used in the numerical experiments reported in Fig. 4. In the transition region, the maxima of the PDF are indicated by black crosses. For the given conditions,

In this section, we analyzed the steady-state statistics of the noise-driven ring of oscillators. We introduced a geometric analogy between the emergent patterns in the ring, quantified by the phase difference , and polar coordinates to enable a graphical description of the (de-)synchronization transition as a function of the linear resistive coupling. The stationary joint PDF and the PDF of the phase differences were obtained by projecting the exact solution of the stationary Fokker–Planck equation onto the phase-locked, uniform-amplitude quasi-limit cycle solutions. The stochastic averaging method was validated by comparing histograms from time series of the fast to those of the averaged system and to the analytical PDFs, showing good agreement over the parameter range considered.

Discussion

In Fig. 13, spectrograms of the modal amplitudes of Bloch modes with different wavenumber b were shown. In Fig. 8 of Ref. [25], similar spectrograms were obtained from acoustic pressure measurements at the can outlets () of a real-world gas turbine with cans. The difference is that in the former, the Bloch modes with dominate the power spectrum, while in the latter, the modes with are dominant. By decreasing to values within the transition region, this clustering of the signal power around non-maximal Bloch wavenumbers b can be reproduced. The corresponding spectrograms are shown in Fig. 16, showing the distribution of the estimated signal power over time t as a function of the normalized frequency . The time traces were computed for () with the fast system (2). Spectrograms for positive and negative b appear similar. The average of the signal power over the shown frequency and time domains, denoted by , is also indicated in the insets. We observe that the signal power is now concentrated around the Bloch modes with . If is further decreased, the signal power shifts toward lower values of b. Note that the inclusion of a coupling nonlinearity K in the model is necessary to observe this effect, as without it, there is no transition region (see Fig. 6). The intermittent energy transfer between different Bloch waves observed in Ref. [25] is qualitatively reproduced by our model. It is an interesting insight that despite the signal power appearing in the spectrograms of the Bloch modes with discrete wavenumber b, the underlying deterministic dynamics are mainly dominated by quasi-steady, non-classical transition modes which are superpositions of CW and CCW rotating waves. This fact cannot be understood from simple observation of the spectrograms, which may suggest that the Bloch modes exclusively dominate the underlying dynamics. The spectrograms reported in Ref. [25] show a shift of the frequency around which the signal power is concentrated with varying Bloch wavenumber b. This feature is not reproduced by our model, which may be because the reactive coupling was neglected. Besides including reactive coupling effects, a topic for future research is to quantify in more detail the time-dependent energy transfer between different Bloch waves observed in both real-world experiments and our model. For such transient phenomena, which are not the focus of this study, differences between the fast system (2) and the averaged system (5) and (6) are expected in the presence of noise, i.e., for .

Fig. 16

Spectrograms of the modal amplitudes for , showing the estimated signal power as a function of time t and the normalized frequency . The modal amplitudes are obtained by decomposing the pressure signals from time series of the fast system (2) with () over 5000 cycles into different Bloch mode components using the discrete Fourier transform, see Eq. (38). Spectrograms for positive and negative b appear similar. The average of the signal power over the shown frequency and time domains, denoted by , is also indicated in the insets

We note that several questions remain open regarding the transition modes observed in Sect. 3 including, but not limited to, the following:The plots in Fig. 14 show that at low noise, the PDF of the phase differences is highly sensitive with respect to changes in . Another topic for future research could be to investigate whether this sensitivity, described by the analytical expression for given in Eq. (37), can be exploited to perform parameter identification of the linear resistive coupling .

Are there different types of transition modes, distinguished by the distribution of along the ring, and can they be systematically classified?

How is the ring size N related to the complexity of the transition modes?

What is their sensitivity with respect to initial conditions? Especially, is it possible to predict, from the initial conditions, if a trajectory will converge to a discrete Bloch mode or to a transition mode?

The following are also topics for future research: Asymmetries have been neglected in this study, but are inherent in real systems. To understand their effects, further investigation is required. Furthermore, in real gas turbines, “true” azimuthal waves may appear in the continuous, annular plenum concurrently with can-annular modes, i.e., apparent rotating waves that emerge due to communicating thermoacoustic modes in individual, but coupled control volumes. The classification of these combined modes is also a topic for future research.

The results reported in Figs. 3 and 4 bear some similarity to the studies found in Refs. [33, 67, 100] where phase-flip bifurcations are studied. The difference is that the referenced works consider transient bifurcations which occur during operation as some parameter is varied, while our investigation is focused on steady-state phenomena with all parameters constant during operation. Extending the current approach to the non-stationary case, and studying transient bifurcations in a similar system could prove a fruitful direction for subsequent studies.

We conclude our discussion with the following remarks: As described in Sect. 3, Bloch modes with positive or negative wavenumber b appear as CW or CCW rotating waves, i.e., waves whose nodal lines spin in CW or CCW direction. A Bloch mode with wavenumber b has exactly |b| nodal lines. However, for finite N, if , there appears to be only a single nodal line which spins in the opposite direction. This can be understood as a spatial analogue of the wagon-wheel effect, where wheels appear to spin in the opposite direction as an effect of finite temporal sampling rates [97]. The effect is illustrated in Fig. 17, where the acoustic pressure along the ring is shown at a fixed time t. The diametric black lines indicate the nodal lines. The radii of the small circles scale with . As the nodal lines of the mode spin in CCW direction, they first hit the oscillators with and then those with etc. Thus, there appears to be only a single nodal line which spins in CW direction, although the underlying Bloch mode is a CCW wave with five nodal lines. Movies of the corresponding time traces are included in the Supplementary Material.

Fig. 17

Spatial analogue of the wagon-wheel effect. The acoustic pressure along the ring is visualized at a fixed time t for and , respectively. The diametric black lines indicate the nodal lines. The radii of the small circles scale with . As the nodal lines of the mode spin in CCW direction, they first hit the oscillators with and then those with etc. Thus, there appears to be only a single nodal line which spins in CW direction, although the underlying Bloch mode is a CCW wave with five nodal lines. Movies of the corresponding time traces are included in the Supplementary Material

Conclusions

In this work, a ring of noise-driven oscillators was analyzed. Deterministic and stochastic averaging was performed to eliminate the fast oscillating terms. Numerical experiments were performed on the noise-free system to motivate the direction of the study. By projecting the potential of the slow-flow variables onto the phase-locked, uniform-amplitude quasi-limit cycle solutions, a compact description of the (de-)synchronization transition in the ring was obtained. These results were adapted to the noise-driven ring of oscillators to derive, in analytical form, the steady-state statistics of the fast system. Special focus was placed on the phase difference , which quantifies emergent patterns along the ring. The stochastic averaging procedure was validated against numerical simulations and analytical results. We have demonstrated that a simple oscillator model with symmetric, purely resistive coupling can reproduce the intermittent energy transfer between different Bloch mode components observed in real-world gas turbines.