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.


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 leanpremixed 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][12][13][14] and [15][16][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 62% efficiency, and which can be supplied with sustainably produced H 2 (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][19][20][21], General Electric [22][23][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 canannular gas turbines. The theoretical studies in Refs. [27,28] use periodic (Bloch-type) boundary conditions to simplify the analysis of reduced-order models of canannular combustors. The effect of resistive (damperlike) 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][31][32][33].

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 U tot : 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 f W/U tot , 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].
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 U tot , 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][47][48][49][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 syn-chronized 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][64][65][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 C j , the acoustic pressure at the can outlet η j , the aperture width W , cross-section area A a 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 q j , which characterize the flame response to acoustic perturbations. Figure 1b shows an abstraction of the same system for N = 12 cans from an axial perspective, the amplitude A j and phase ϕ j of the acoustic pressure at the can outlet η j , 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 2π/N along the ring, and that j = 1 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 U tot 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 A a /A is smaller than the ratio of effective aperture depth to can length (d + √ A a )/L, respectively: A a /A < (d + √ A a )/L < 1. 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  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 C j , the acoustic pressure at the can outlet η j , the aperture width W , cross-section area A a 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 q j , which characterize the flame response to acoustic perturbations. b Abstraction of the system illustrated in (a), viewed from an axial perspective, for N = 12 cans. Shown are the amplitude A j and phase ϕ j of the acoustic pressure at the can outlet η j , respec-tively, 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 U tot 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 "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 ψ 0 of the can with corresponding eigenfrequency ω 0 . In a first approximation, we assume that the mode shape of ψ 0 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 ω ≈ ω 0 , where ω = 2π f and f is the frequency of oscillation, the acoustic pres-sure field p j in the jth can is approximated by the unimodal projection where η j 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 ψ 0 (L) = 1 (see Fig. 1a). In this case, η j is equivalent to the acoustic pressure p j 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 ν > 0, 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 smallamplitude 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][74][75][76], chemistry [77][78][79] and biology [80][81][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][84][85] and Kuromoto's classic model [86], which has become a paradigmatic example of phase synchronization in coupled oscillators [87][88][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][69][70]. Such loworder 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 ξ j of equal intensity Γ with the autocorrelation functions ξ j ξ j,τ R = R ξ j (t)ξ j (t + τ )dt = Γ δ(τ), where · R 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 canannular combustor: The sign before ϑ is positive because the coupling is amplifying if λ < 0 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 λ = 0 with nonzero saturation κ = 0 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.
We define the coordinate change to amplitude-phase variables {A j , ϕ j }, j = 1, . . . , N , as follows: where η j = A j cos (ϕ j + ω 0 t),η j = −ω 0 A j sin (ϕ j + ω 0 t) and arctan(y, x) is the two-argument arctangent function. Assuming that (A) the variation of A j and ϕ j 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: is the phase difference and ζ j and χ j are 2N uncorrelated white noise sources with the same noise intensity Γ /2ω 2 0 , respectively. The above definitions of ϕ j and Δ j 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 ϕ j ≡ mod(ϕ j , 2π) and Δ j ≡ mod(Δ j + π, 2π) − π . These operations, which leave Eqs. (5) and (6) unchanged, are equivalent to adding integer multiples of ±2π to ϕ j and Δ j so that the sums lie in the domains [0, 2π ] 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 j, m = 1, . . . , N and the potential V is defined as follows: Equations (7) and (8) imply that trajectories (A j , ϕ j ), j = 1, . . . , N , are stochastically perturbed and attracted toward lower values of the potential V 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) P(η m ,η m /ω 0 , t), m = 1, . . . , N , 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 t → ∞.
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 [0, t end ], where 1/ν is the relaxation time (∼ time is takes to reach the steady state) and τ a is the correlation time of the signal "a". The correlation time is defined as τ a = ∞ 0 + ξ j ξ j,τ R dτ , where (·) + denotes the one-sided limit from above (see Ref. [93], p. 65, Eq. (3.61) with τ = t 2 − t 1 ). 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 ξ j , ζ j and χ j will generally have nonzero correlation times. However, because we model ξ j , ζ j and χ j as white noise signals with (theoretically) vanishing correlation times, we assume throughout this work that 1/ν max(τ ξ j , τ A j , τ ϕ j ), 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, Γ = 0. Under the averaging assumptions (A) and (B), the limit cycles of the deterministic part of the fast system (Eq. 2 with ξ j ≡ 0) are given by fixed points {Ȧ j = 0,φ j = 0} of the slowflow system (5) and (6) with ζ j ≡ χ j ≡ 0. We find that these fixed points describe solutions with uniform amplitude along the ring and equal phase difference between neighboring oscillators: The periodicity condition ϕ j ≡ ϕ j+N 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 Solutions with |b| > 0 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, , 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 (Δ * b > 0) 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 b = 0 is the in-phase synchronized solution, for which all oscillators reach maximum and minimum pressure simultaneously. For even N , we call the solution with b = N /2 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 mod(N , |b|) = 0, 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. Fig. 2 Illustration of the limit cycle solutions corresponding to b = 0 and b = ±1 for N = 3 oscillators. Shown are the values of the phase difference Δ * b and the phase distribution ϕ j 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 b = 0 is the synchronized solution, for which all oscillators are in phase For vanishing coupling λ = ϑ = 0, the system (2) corresponds to N decoupled van der Pol oscillators, whose weak-nonlinearity limit cycle solution [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)- where A j and Δ j 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]) yieldṡ This 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 λ < ν/2 sin(π b/N ) 2 , the linear stability of the corresponding limit cycles, demonstrated by Eq. (13), shows that all Bloch modes are local minima of the potential V 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 N = 3 oscillators, the limit cycle solutions corresponding to b = 0 and b = ±1 are shown in Fig. 2. Shown are the values of the phase difference Δ * b and the phase distribution ϕ j along the ring.
Substituting Eqs. (10)-(12) into the slow-flow system (5) and (6) yields the limit cycle amplitudes where we have defined Λ = 2λ/ν is the normalized linear resistive coupling and K = 16ϑ/κ the normalized coupling nonlinearity. Limit cycles for a given value of the Bloch wavenumber b can only exist if otherwise the amplitude A * b 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 P(Ω) of the 2Ndimensional phase space Ω, where defined by P(Ω) : 1}. These quasilimit cycle solutions include the case of true limit cycles, which are fixed points of the slow-flow system (5) and (6), for k( j) ≡ 1. However, we keep k general because, as we show below, quasi-steady solutions with uniform amplitudes A j ≡ A * β ≈ const. and uniform absolute phase differences |Δ j | ≡ |Δ * β | = |πβ/N | ≈ const. play a crucial role in the nonlinear dynamics of the oscillator ring. The projection of the potential V defined by Eq. (9) onto the submanifold P(Ω) reads Normalizing this result, and using the fact that for a thermoacoustic instability, sign(ν) = 1, yields where U P = κV P /2N |ν| 2 is the normalized projection of the potential V onto the phase-locked quasi-limit cycle solutions. Note U P is independent of the signs in k and is symmetric with respect to positive and negative B. We show below that U P 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 U P , given by Eq. (21), around B = 0 (synchronized state) is and the expansion around B = ±π/2 is For Λ > 0 (dissipative coupling) and B ≈ 0, we consider Eq. (22), which shows that the global minimum of the potential V is at B = 0, i.e., the synchronized state. Conversely, for Λ < Λ c = −K (strong amplifying coupling) and B ≈ ±π/2, we consider Eq. (23) which shows the global minima of V are at 2B = Δ * ±N /2 = ±π (or Δ * ±floor(N /2) for odd N ), corresponding to anti-phase synchronization.
For 0 > Λ > Λ c , which we call the "transition region" in the following, the minima of the projected potential U P , given by Eq. (21), lie between B = 0 and B = ±π/2. This parameter range is further investigated below. By computing the zero level set of the slope of U P with respect to B, we obtain the curve parameterizing the minima of U P in the transition region as a function of Λ: Numerical examples were performed and analyzed to explore the projected potential U P , 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 In  Fig. 3a-c, the distribution of the projected potential U P as a function of Λ and 2|B| = |Δ * b | is shown. These insets illustrate how a nonzero coupling nonlinearity K > 0 leads to a transition region for 0 > Λ > Λ c = −K , bounded by the dotted lines, which separates the domains where the inphase (In-phase) and the anti-phase synchronized state (Anti-phase) are the global minimum of the potential V, given by Eq. (9), respectively. The solid black curve indicates the location of the minimum of U P 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 U P , given by Eq. (21), for Λ = 0.5, −0.5 and −0.75, 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 U P as a function of 2|B|. In Fig. 3e, f, (Λ = 0.5), h, i (Λ = −0.5), it is shown how, in a ring of N = 3 oscillators with dissipative coupling, all trajectories converge to the synchronized state and for strong amplifying coupling with Λ < Λ c = −K , all trajectories converge to rotating wave solutions, with a symmetric distribution between CW (Δ j < 0) and CCW (Δ j > 0) rotating waves. The number of realizations that converge to either solution class is indicated above the Δ j -curves. As shown in Fig. 3k Fig. 3k shows a comparison between a single realization of η 1 (black) computed from the fast system (2) and the amplitude A 1 (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 Δ j > 0 between 600 and 800 cycles. The corresponding Bloch modes are visualized in Fig. 3n. Similar aggregating behavior around the synchronized state with b = 0 is observed for Λ > 0, see Fig. 4. Figure 3 shows that while a change in Λ from 0.5 to −0.5 leads to O(1) changes in the phase differences Δ j , the corresponding change in the amplitudes A j in the steady state is small compared to A * 0 . To better understand this, note that on the one hand, for Λ > 0, the synchronized state is globally attracting, and the system aggregates around low-order Bloch modes b ≈ 0, which implies a small influence of Λ and K on the limit cycle amplitude A * b , see Eqs. (15)- (17). On the other hand, for Λ < 0 and 0 < K 1, 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 (|Λ| 1, K 1), the coupling-induced changes in the amplitude are in general assumed to be small. For Λ < Λ c = −K , solutions with |b| > 0 reach limit cycle amplitudes A * b 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 Δ j , whose distribu-tion quantifies emergent patters of the acoustic pressure η j 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 U P , given by Eq. (21), as a function of Λ and 2|B| = |Δ * β |. The colored crosses in 4a mark the location of the minimum of the projected potential U P , 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 U P over the phase difference 2|B|. The positions of the minima of U P 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 Δ j of the oscillators with j = 1, 2 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 |Δ 1 | over all realizations at 5000 cycles, |Δ 1 (t end )| R , is displayed in the same insets. In Fig. 4h-m, the smaller insets in the upper right corner show the relative share R(Δ 1 )/R tot of the total realizations R tot of Δ 1 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 (100% and 95%, 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 b = 4.
We observe in Fig. 4 that outside the transition region, for Λ > 0 or Λ < Λ c , when the global minimum of the potential V, 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 Λ > 0 and Λ < 0, nearly symmetric unimodal and bimodal distributions of the phase difference Δ 1 over Δ * b and Δ * β are observed in the numerical experiments, respectively.
Inside the transition region, for 0 > Λ > Λ c , we observed in our numerical experiments that trajec- In-phase Anti-phase 1. 16 2.50 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 V, given by Eq. (9), respectively, for an oscillator ring with N = 12 members. a Distribution of the projected potential U P , given by Eq. (21), as a function of Λ and 2|B| = |Δ * β |. The colored crosses indicate the location of the minimum of the projected potential U P , given by Eq. (21), over B for six evenly spaced values of the linear resistive coupling Λ, respectively. The corresponding distributions of the projected potential U P over the phase difference 2|B| is shown in (b)-(g). The positions of the minima of U P 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 Δ j , j = 1, 2, 3, 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 |Δ 1 | over all realizations at 5000 cycles, |Δ 1 (t end )| R , is displayed in the same insets. In (h)-(m), the smaller insets in the upper right corners show the relative share R(Δ 1 )/R tot of the total realizations R tot of Δ 1 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 (100% and 95%, 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 b = 4 tories 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 A j and phase differences Δ j '. In the transition region, the mean |Δ 1 (t end )| R over all realizations matches well with the minimum Shown is a the variation of the amplitudes A j and b the phase differences Δ j ≡ mod(ϕ j − ϕ j−1 + π, 2π)− π for a single realization over 300 cycles in the steady state, c the distribution of the projected potential U P over 2|B| = |Δ * β | (magnified view of Fig. 4f with a different colormap) and d the phases ϕ j 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 |Δ j | near the minimum of U P over 2|B| (see Fig. 4k ,l). The phase distribution ϕ j along the ring shown in (c) corresponds to a standing wave-type pattern of the acoustic pressure  Fig. 5a, b are the variation of the amplitudes A j and the phase differences Δ j = mod(ϕ j − ϕ j−1 + π, 2π) − π for j = 1, . . . , N , respectively, for a single realization over 300 cycles in the steady state. Figure 5c, d shows the distribution of the projected potential U P over 2|B| = |Δ * β | (magnified view of Fig. 4f with a different colormap) and the phases ϕ j 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 Δ j perform slow, small-scale periodic motion over time. Because the transition modes feature both positive and negative phase differences Δ j , these solutions are superpositions of CW and CCW rotating waves which appear as standing wave-type patterns of the acoustic pressure η j along the ring. After a transient, the amplitudes A j all remain quasi-steady, performing small-scale, periodic motion around A * 0 . Time traces of other transition modes appear similar, with more or less pronounced asymmetry between different oscillators, but with absolute phase differences |Δ j | near the minimum of U P over 2|B|, see Fig. 4k, l. Animated examples of transition mode time traces (A j , ϕ j , Δ j and η j ), computed from the slow-flow system (5) and (6), are contained in the Supplementary Material.
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 U P , 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 V, 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 Of course this solution can never be exactly observed, because these quasi-limit cycles are not fixed points of the deterministic part of the slowflow system given by Eqs. (5) and (6). The normalized potential U P 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 {A j ≡ A * β , Δ j ≡ k( j)Δ * β }: where the fact that in the steady state,φ j ≈ 0 was used and ν β , κ β are obtained by replacing b with β in Eqs. (16) and (17), respectively. For a distribution of Δ j with varying signs in k along the ring, Eq. (26) implies unsteadiness of the Bloch wavenumber β. This unsteadiness carries over to the equation foṙ A * β 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 V of the amplitude-phase system (5) and (6) was projected onto the phase-locked, uniformamplitude quasi-limit cycle solutions to obtain the normalized projected potential U P , which compactly describes the (de-)synchronization transition, dividing the parameter space into two domains where the inphase and the anti-phase synchronized state are the global minimum of the potential V, respectively. These domains are separated by a transition region, whose size is proportional to the coupling nonlinearity K and which vanishes for K = 0. Notably, U P is both nondimensional and independent of the number of oscillators N .

Steady-state statistics of the noise-driven system
Building on our analysis of the projected potential U P in the previous section, we begin our investigation of the noise-driven system by studying the behavior of U P under varying coupling nonlinearity K and noise intensity Γ . For Γ = 0, 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 Γ = 0, 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 A j ≡ A * β and absolute phase differences |Δ j | ≡ |Δ * β |. The normalized projected potential U P depends now also on the normalized noise intensity G = κΓ /8|ν| 2 ω 2 0 and is computed for different values of K and G by substituting the limit cycle amplitude (27) and Δ j ≡ Δ * β into the potential V, given by Eq. (9), and dividing the result by the factor 2N |ν| 2 /κ for normalization. For compactness, the resulting expression is only listed implicitly here: where A * β and Δ * β are obtained by replacing b with β in Eqs. (12) and (27), respectively.
The distribution of the projected potential U P as a function of Λ and 2|B| = |Δ * β | is studied in Fig. 6. The dashed black lines mark the transition region 0 > Λ > Λ c . Note that Λ c 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 V, given by Eq. (9), respectively. The solid black curves in the transition region parametrize the minima of U P with respect to 2B, defined by ∂U P /∂ B = 0, over Λ. In the bottom right inset in Fig. 6, Λ c = −0.356. We observe that U P 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 K > 0, the critical value of the linear resistive coupling Λ c , 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 A j are in general not all simultaneously close to some amplitude A * β and the projection onto the limit cycle solutions loses its validity. Despite this, the figure illustrates correctly the monotonous trends of U P , given by Eq. (28), with increasing coupling nonlinearity K and noise intensity G. The curve parametrizing the minima of U P with respect to 2B over Λ is obtained by numerically computing the zero level set of ∂U P /∂ B.
For t → ∞, in the steady state, the probability P of the fast system (2) being within a domain D in the 2N -dimensional phase space of the fast system (2), Σ = {η 1 ,η 1 , . . . , η N ,η N }, is defined as η 1 ,η 1 , . . . , η N ,η N can be deduced from the stationary Fokker-Planck equation (see Ref. [95]) of the slow-flow system (5) and (6), given by where the newly introduced (slow) variables The derivation of P ∞ is detailed in Appendix B. The final result is the potential V is given by Eq. (9) and N is a normalization constant such that P(D → Σ, t) = 1. 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 P ∞ onto the phase-locked, uniform-amplitude quasi-limit cycle This result implies that on the submanifold of phaselocked, 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 P ∞ P . The PDF of the phase differences P P,Δ is obtained by evaluating the limit cycle amplitude A * β in Eq. (27) in terms of Δ * β using Eqs. (16) and (17) with b → β and substituting the result into the potential V in Eq.

Anti-phase
In-phase CCW waves  (36). The result is, up to normalization, To graphically compare the joint PDF P ∞ P (A * β , Δ * β ) to normalized histograms of A j and Δ j obtained from time traces, we propose a geometric analogy between emergent patterns in the oscillator ring, quantified by the phase difference Δ j , and polar coordinates. The analogy is shown in Fig. 7. Possible nonnegative discrete Bloch wavenumbers b ≥ 0 are shown for N = 12 oscillators. The x-axis (b = 0) corresponds to in-phase and the y-axis (b = N /2 = 6) 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 Λ = −0.5 (Λ c = −0.356), so that the considered condition lies outside the transition region. The oscillators j = 1, 2, 3 are indicated by blue, red and cyan colors, respectively. In Fig. 8a, we denoted by dB f,t , is also indicated in the inset. e Comparison of the analytical PDF of the phase differences P ∞ P,Δ (red) with a histogram of Δ 1 (blue) from the time trace shown in (b). f Histograms of A 1 and Δ 1 , represented using the geometric analogy described in Fig. 7. Spectro-and histograms of other oscillators appear similar observe that the amplitudes A j fluctuate stochastically near the deterministic limit cycle amplitude A * 6 , which is indicated in green. The phase difference Δ 1 , shown in Fig. 8b performs similar fluctuations and jumps intermittently between ±π . Time traces of other oscillators appear similar. The phase distribution ϕ j along the ring after 5000 cycles is shown in Fig. 8c. Figure 8d shows the spectrogram of the reconstructed fast oscillating signal η 1 (t), i.e., the estimated signal power L η 1 η 1 ( f ) f , where L η 1 η 1 is the power spectral density of the signal η 1 and Δf 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 ω/ω 0 . We observe that the signal power is mainly concentrated near the natural eigenfrequency ω 0 . The average of the signal power over the shown frequency and time domains, denoted by dB f,t , is also indicated in the inset. Figure 8e shows a comparison of the analytical PDF of the phase differences P ∞ P,Δ (red) with a histogram of Δ 1 (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 Δ * ±6 = ±π . 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 j = 1, . . . , N . In Fig. 9b, the distribution of the amplitudes A j along the ring at the end time is shown. We observe that the assumption of slowly varying, uniform amplitudes is approximately satisfied.
The stochastic averaging method is validated against numerical simulations in Fig. 10 by comparing histograms of A 1 and Δ 1 , 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 Λ < 0. At high noise intensity, the PDF becomes more spread out and less discernible patterns are visible. Note that the amplitude A * 0 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 Cyles A j 5000 0 0 Fig. 9 a Amplitudes A j of the time trace shown in Fig. 8 over the last 5 cycles. b Distribution of the amplitudes A j along the ring after 5000 cycles

Normalized linear coupling
Normalized noise intensity Av.

Low noise High noise
Amplifying Dissipative Fig. 10 Validation of the stochastic averaging method against numerical simulations. Histograms of A 1 and Δ 1 , 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 A * 0 depends on the noise intensity G and changes from the left to the right column 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. In Fig. 11, the stochastic averaging method is validated against analytical results by comparing histograms of A 1 and Δ 1 computed from the fast system (2) (Sim.) to the joint PDF P ∞ P , 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 The spectrograms of η j , shown in Fig. 8d for j = 1, can be decomposed into spectrograms of different rotating wave components by defining the modal amplitude η b of the Bloch mode with wavenumber b using the discrete Fourier transform [25]: In Fig. 13, we plot spectrograms of η b for b = 0, . . . , 6. Spectrograms for positive and negative b appear similar. Random, intermittent energy transfer between different modes is observed. The average of the signal  Fig. 14 Comparison of the stationary PDF of the phase differences P ∞ P,Δ (Δ * b ), given by Eq. (37) (Theory), to histograms of Δ 1 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 power over the shown frequency and time domains, denoted by dB f,t , is also indicated in the insets. As expected, for the given conditions (Λ = −0.5, Λ c = −0.356), the anti-phase synchronized state with b = ±6 dominates the power spectrum.
In Fig. 14, the stationary PDF of the phase differences P ∞ P,Δ , given by Eq. (37), is compared to histograms of Δ 1 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.
The evolution of the joint PDF P ∞ P , given by Eq. (36), is visualized in Fig. 15 for varying Λ (Λ c = −0.356), with G = 2.88 × 10 −1 ("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 P ∞ P reproduces the transition from a unimodal to a bimodal distribution as Λ is varied from 0.5 to −0.5, which was also observed in the numerical experiments on the deterministic system reported in Figs. 3 and 4. 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 Δ j , 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 P ∞ P and the PDF of the phase differences P ∞ P,Δ were obtained by projecting the exact solution of the stationary Fokker-Planck equation onto the phaselocked, 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 η b 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 (∼ η j ) of a real-world gas turbine with N = 12 cans. The difference is that in the former, the Bloch modes with b = 5, 6 dominate the power spectrum, while in the latter, the modes with b = 3, 4 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 ω/ω 0 . The time traces were computed for Λ = −0.3 (Λ c = −0.356) 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 dB f,t , is also indicated in the insets. We observe that the signal power is now concentrated around the Bloch modes with b = 4, 5. 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 Γ = 0.
We note that several questions remain open regarding the transition modes observed in Sect. 3 including, but not limited to, the following: -Are there different types of transition modes, distinguished by the distribution of Δ j 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 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 P ∞ P,Δ given in Eq. (37), can be exploited to perform parameter identification of the linear resistive coupling λ.
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  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 steadystate phenomena with all parameters constant during operation. Extending the current approach to the nonstationary 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 mod(N , |b|) = 0, 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 η j 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 |η j |. As the nodal lines of the b = 5 mode spin in CCW direction, they first hit the oscillators with j = 3, 9 and then those with j = 4, 10 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 phaselocked, 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 Δ j , 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.
Funding Open access funding provided by Swiss Federal Institute of Technology Zurich.

Data availability statement
The datasets used for generating the plots and results in the present study can be directly obtained from the numerical simulation of the related mathematical equations in the manuscript.

Conflict of interest
The authors declare that they have no conflict of interest.

Supplementary material Animated examples of time traces of
A j , Δ j , ϕ j and η j , computed from the noise-free (G = 0) and noise-driven (G = 2.88 × 10 −1 ) slow-flow system (5) and (6) over five acoustic cycles in the steady state for different values of Λ, are included in the supplementary material. For visualization purposes, the color bar of the time traces of η j in the noise-driven case is cut off at ±1.2 A * 0 . As in Fig. 17, in the plots showing the distribution of the acoustic pressure η j along the ring, the radii of the small circles scale with |η j |.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/. (48) for k = j − 1, j, j + 1. The expressions in Eqs. (45) and (47) are valid for any nonlinear function f . We now restrict ourselves the case where the function f is given by Eq. (42). Averaging is applied to keep only contributions that are relevant on the slow time scale associated with the modulation of the amplitudes A j and phases ψ j . First, following the classic treatment of Krylov and Bogoliubov [91], deterministic averaging is performed on f 1 and f 3 . To facilitate the averaging method, we rewrite these functions by expanding the trigonometric terms. This yields and where Δ j = φ j − φ j−1 = ϕ j − ϕ j−1 and "f.o.t." stands for "fast oscillating terms." Note that the quantities φ j + φ j±1 = ϕ j + ϕ j±1 + 2ω 0 t as well as 2φ j = 2ϕ j + 2ω 0 t define arguments of fast oscillating terms. We assume that the period of an acoustic cycle T = 2π/ω 0 is small compared to the slow time scale 1/ν of the modulation of the amplitude and phase signals. By taking the averages of Eqs. (49) and (50), f 1 (A k , ϕ k ) T and f 3 (A k , ϕ k ) T , respectively, over a time interval {t 0 , t 0 + T } of length T with arbitrary start time t 0 , one removes the fast oscillating terms.
Stochastic averaging of the terms f 2 and f 4 can be performed in different ways, two of which are presented in Refs. [93] (Vol. 2, pp. 105-113) and [94], respectively. The former method, due to Stratonovich, is adopted here. In our derivation, we follow the presentation of a similar result (stochastic averaging of a single thermoacoustic van der Pol oscillator) given in Appendix A of Ref. [70].
For compactness, we use the notation f τ = f (t +τ ) and the dependence of φ j , ϕ j and A j on time is suppressed below. We introduce the random variables σ j,1 = −ξ j sin φ j = f 2 ω 0 and σ j,2 = −ξ j cos φ j = f 4 A j ω 0 . We begin by estimating the expected values of σ j,1 and σ j,2 . We consider a time interval of length θ satisfying 1/ν θ max(τ ξ j , τ A j , τ ϕ j ) ∀ j, where τ a denotes the correlation time of the signal "a". We note that although the signals ξ j , ζ j and χ j are modeled as white noise sources with (theoretically) vanishing correlation times, in reality and in numerical simulations, these signal will generally have nonzero, albeit small correlation times. As a first step, we approximate σ j,1 and σ j,2 as follows: where ϕ j,−θ =ϕ j (t − θ) and ϕ j = ϕ j − ϕ j,−θ . Using the fact that ϕ j,−θ is not correlated with ξ j since θ max(τ ξ j , τ ϕ j ), we write Considering that ξ j is a stationary process and that 1/ν θ , i.e., the variation of the amplitudes A j and phases ϕ j over the considered time interval is negligible, we obtain (56) and σ j,2 R = − 0 −θ ξ j ξ j,τ R sin (ω 0 t + ϕ j ) cos (ω 0 (τ + t) + ϕ j )dτ ω 0 A j .
Expanding the products of cosines in Eqs. (57) and (58) and neglecting the fast oscillating terms gives where L ξ j ξ j is the power spectral density of the noise signal ξ j . We proceed by estimating the correlation function of the zero mean processes σ j,1 = σ j,1 − σ j,1 R and σ j,2 = σ j,2 − σ j,2 R = σ j,2 . In Eqs. (51) and (52), the mean values are determined by the second terms on the RHS and the fluctuating components result from the first terms: σ j,2 ≈ −ξ j cos (ω 0 t + ϕ j,−θ ).
Since θ max(τ ξ j , τ ϕ ), it is assumed that these processes are delta-correlated. Considering first σ j,1 , one can write where the integral in Eq. (63) is expressed as Because the correlation time τ ξ j is much shorter (theoretically zero) than the oscillation period 1/ω 0 , we can neglect the fast oscillating terms and write Following the same derivation for σ j,2 yields the correlation function σ j,2 σ j,2,τ R = (Γ /2)δ(τ ). With the above results, we can write for j = 1, . . . , N , where ζ j and χ j are 2N uncorrelated white noise sources with the noise intensity Γ /2ω 2 0 . Note that the mean of the averaged process f 2 (A, φ) R is not zero and that the stochastic averaging procedure gives rise to a deterministic term which is inversely proportional to A j .
Combining Eqs. (49), (50), (65) and (66), together with the assumption that the variation of A j and ϕ j over an acoustic period is negligible, yields the slowflow amplitude and phase dynamics of the fast system (2), given by the following set of equations: anḋ We assume that the Jacobian J f = ∂ f /∂ y is invertible and can be written as

J f ( y) = Q( y)h( y),
where Q = Q −T is an orthogonal rotation matrix with |det( Q)| = 1 and h is a diagonal matrix. This implies h is related to the metric tensor g = J f T J f by h T ( y)h( y) = g( y).
Now consider a stochastically perturbed gradient system, given by the following Langevin equation in rectilinear coordinates x: where n contains additive white noise components in the directions of the state variables x, which are all assumed to have the same noise intensity Γ /2ω 2 0 . For a system of the type of Eq. (72), the stationary Fokker-Planck equation [95] is given by Eq. (31), which, for P ∞ = 0, has the following normalizable solution: where N ∈ R is a normalization constant. In ycoordinates, the probability of the stationary state of the system being within some domain f −1 (D) is given by where the fact that |det( J f )| = |det(h)| was used. With this, the stationary PDF reads, in y-coordinates, Under the mapping (69), the Langevin equation (72) becomeṡ whereñ = Q T n contains the noise components in the direction of the y-coordinates. The solution of the stationary Fokker-Planck equation P ∞ ( y) associated with the transformed Langevin equation (76) is given by Eq. (75).
The matrix Q = J f h −1 is a block-diagonal matrix of orthogonal matrices, and is thus itself orthogonal. Hence, Eq. (70) is satisfied for the present system.
To derive the solution of the stationary Fokker-Planck equation (31), we note that the slow-flow system (5) and (6) is equivalent to a Langevin equation of the form given in Eq. (76), where y = (A 1 , ϕ 1 , . . . , A N , ϕ N ), and the potential V is defined by Eq. (9). We find from Eq. (82) that With this result, the stationary solution of the Fokker-Planck equation associated with the slow-flow system (5) and (6) can be directly computed from Eq. (75): We now consider the mapping z = d(x) from x to the fast z coordinates given by d(x) = [A 1 cos(φ 1 ), −ω 0 A 1 sin(φ 1 ), . . . , where φ j = ϕ j + ω 0 t and the amplitude-phase coordinates (A j , ϕ j ) are expressed in terms of U j and V j using Eqs. (87) and (88) where and the dependence of φ j on U j and V j is suppressed for compactness. The modulus of the determinant of J d is the positive real constant | det( J d )| = ω N 0 , which is absorbed into the normalization factor N . Since and by the definition of P ∞ (x) in Eq. (86), the stationary PDF of the fast system (2) is given by the potential V is given by Eq. (9) and N is a normalization constant.