Synchronization and chimeras in a network of four ring-coupled thermoacoustic oscillators

Abstract

We take a complex systems approach to investigating experimentally the collective dynamics of a network of four self-excited thermoacoustic oscillators coupled in a ring. Using synchronization metrics, we find a wide variety of emergent multi-scale behaviour, such as (i) a transition from intermittent frequency locking on a $\mathbb {T}^{3}$ quasiperiodic attractor to a breathing chimera, (ii) a two-cluster state of anti-phase synchronization on a periodic limit cycle, and (iii) a weak anti-phase chimera. We then compute the cross-transitivity from recurrence networks to identify the dominant direction of the coupling between the heat-release-rate ($q^{\prime }_{\mathbb {X}}$) and pressure ($p^{\prime }_{\mathbb {X}}$) fluctuations in each individual oscillator, as well as that between the pressure ($p^{\prime }_{\mathbb {X}}$ and $p^{\prime }_{\mathbb {Y}}$) fluctuations in each pair of coupled oscillators. We find that networks of non-identical oscillators exhibit circumferentially biased $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ coupling, leading to mode localization, whereas networks of identical oscillators exhibit globally symmetric $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ coupling. In both types of networks, we find that the $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ coupling can be symmetric or asymmetric, but that the asymmetry is always such that $q^{\prime }_{\mathbb {X}}$ exerts a greater influence on $p^{\prime }_{\mathbb {X}}$ than vice versa. Finally, we show through a cluster analysis that the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ interactions play a more critical role than the $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions in defining the collective dynamics of the system. As well as providing new insight into the interplay between the $p^\prime_{\mathbb{X}}\text{--}p^\prime_{\mathbb{Y}}$ and $p^\prime_{\mathbb{X}}\text{--}q^\prime_{\mathbb{X}}$ coupling, this study shows that even a small network of four ring-coupled thermoacoustic oscillators can exhibit a wide variety of collective dynamics. In particular, we present the first evidence of chimera states in a minimal network of coupled thermoacoustic oscillators, paving the way for the application of oscillation quenching strategies based on chimera control.

1. Introduction

Over the past several decades, tightening emission standards have prompted gas turbine manufacturers to deploy lean-premixed combustion (Lieuwen & Yang 2005). However, this strategy is known to be susceptible to thermoacoustic instability (O'Connor, Hemchandra & Lieuwen 2016). Thermoacoustic instability is typically caused by constructive interactions between the heat-release-rate (HRR) oscillations of an unsteady flame and the pressure oscillations of the combustor (Culick 2006). Such interactions can arise from various mechanisms, such as hydrodynamic instabilities (Poinsot et al. 1987), equivalence ratio fluctuations (Lieuwen & Zinn 1998) and entropy waves (Candel 2002). If the HRR and pressure oscillations are sufficiently in phase, energy can be transferred from the flame to the acoustic field via the mechanism of Rayleigh (1945), resulting in self-excited flow oscillations whose frequencies are often close to those of the natural acoustic modes of the system (Culick 2006). Thermoacoustic instability can also arise from an intrinsic feedback mechanism involving upstream propagating acoustic waves emitted by the flame itself (Hoeijmakers et al. 2014; Emmert, Bomberg & Polifke 2015), with no requirement for acoustic reflection at the combustor boundaries (Silva et al. 2015; Buschmann, Mensah & Moeck 2020). Regardless of the specific feedback mechanism, however, thermoacoustic oscillations can exacerbate thermomechanical stresses and flame blowoff/flashback, limiting the performance and service life of the overall combustion system (Poinsot 2017).

The analysis, prediction and control of thermoacoustic oscillations have been the subject of extensive research (Candel 2002; Dowling & Morgans 2005; Lieuwen & Yang 2005; Culick 2006; Sujith, Juniper & Schmid 2016; Poinsot 2017; Juniper & Sujith 2018; Heckl 2019; Polifke 2020; Schuller, Poinsot & Candel 2020). Most of the existing literature has focused on single combustors because these have simple geometries and well-defined boundary conditions while still being governed by the same physical mechanisms that are responsible for the thermoacoustic feedback loops described above. In many industrial gas turbines, however, there exists not just a single combustor in isolation, but multiple combustors coupled to one another (Luque et al. 2015). Such multi-combustor systems, known as can-annular systems (Bethke et al. 2002), can host a variety of thermoacoustic modes, whose stability and dynamics are determined not only by the flame–acoustic interactions occurring within each individual combustor (i.e. intra-combustor interactions), but also by the bidirectional acoustic interactions occurring between directly/indirectly coupled combustors (i.e. inter-combustor interactions) (Mongia et al. 2003; Kaufmann et al. 2008; Luque et al. 2015; Farisco, Panek & Kok 2017; Bonciolini & Noiray 2019; Ghirardo et al. 2019; von Saldern, Moeck & Orchini 2021a). Although considerable research exists on intra-combustor flame–acoustic interactions (Candel 2002; Dowling & Morgans 2005; Lieuwen & Yang 2005; Culick 2006; Sujith et al. 2016; Poinsot 2017; Juniper & Sujith 2018; Heckl 2019; Polifke 2020; Schuller et al. 2020), less is known about inter-combustor acoustic–acoustic interactions. In this experimental study, we take a complex systems approach to investigating both types of interactions, with a focus on how they give rise to collective dynamics in a self-excited thermoacoustic system consisting of four turbulent lean-premixed combustors coupled in a ring configuration.

1.1. Synchronization and chimeras in coupled thermoacoustic oscillators

The thermoacoustics of multi-combustor systems can be analysed conveniently in the framework of mutual synchronization (Pikovsky, Rosenblum & Kurths 2003; Balanov et al. 2008). This involves treating each combustor as an individual self-excited oscillator (e.g. undergoing limit-cycle motion) and examining how the mutual coupling of those oscillators can change their phase and amplitude dynamics (Sujith & Unni 2020, 2021). This approach has underpinned numerous studies on coupled thermoacoustic oscillators. For example, using concepts from mutual synchronization, Thomas et al. (2018a,b) investigated numerically the effects of dissipative coupling, time-delay coupling, and external noise on the transition to amplitude death in a coupled Rijke-tube model. Using similar concepts, Dange et al. (2019) observed experimentally in-phase/anti-phase synchronization, phase-flip bifurcations and partial amplitude death in two coupled Rijke tubes powered by electric heaters. Similarly, Biwa, Tozuka & Yazaki (2015), Hyodo & Biwa (2018) and Hyodo, Iwasaki & Biwa (2020) observed experimentally in-phase/anti-phase synchronization and amplitude death in two coupled thermoacoustic oscillators powered by electric heaters and laminar Bunsen flames.

The above studies focused only on laminar systems, but practical combustors are almost always turbulent (Lieuwen 2012). Recognizing this, Jegal et al. (2019) and Moon et al. (2019) have investigated experimentally the mutual synchronization of two turbulent lean-premixed combustors coupled via a cross-talk tube. They found a variety of collective dynamics, including in-phase/anti-phase synchronization, desynchronization associated with quasiperiodicity, and complete/partial amplitude death. In a follow-up study, Moon et al. (2020a) investigated how changing the dimensions of the cross-talk tube can affect the dissipative and time-delay coupling between the two combustors as well as their collective dynamics. Recently, these dynamics were modelled phenomenologically by Guan et al. (2021) using two canonical self-excited temporal oscillators (Van der Pol oscillators) coupled via dissipative and time-delay terms, demonstrating the universality of the observed synchronization phenomena. This was inspired by the work of Bonciolini & Noiray (2019), who used two coupled stochastically driven nonlinear oscillators to model the synchronization of thermoacoustic modes in two sequential turbulent combustors.

The use of turbulent combustors in the above studies marks a step closer to practical conditions. However, those studies were limited to only two combustors at a time, whereas heavy-duty land-based gas turbines typically contain a larger number of combustors coupled in a ring configuration. The synchronization of thermoacoustic modes in can-annular and annular combustors has attracted growing interest in recent years. For example, Moeck et al. (2019) used an oscillator model to explore the nonlinear coupling between azimuthal and axisymmetric modes in annular combustors. They found that synchronization between a standing azimuthal mode and an axisymmetric mode can occur if the two modes exhibit similar resonance frequencies and growth rates. In can-annular combustors, the existence of clusters of thermoacoustic modes with similar frequencies was first revealed by Ghirardo et al. (2019) and was subsequently explained by von Saldern, Orchini & Moeck (2021b) using Bloch-wave theory. On the experimental front, Moon et al. (2020b) and Moon, Yoon & Kim (2021) recently investigated the thermoacoustics of four ring-coupled turbulent combustors, but they did not do so in a synchronization framework, leaving open questions about the simultaneous interactions between different thermoacoustic modes and their collective dynamics. In particular, the thermoacoustics of multi-combustor systems is known to depend on both (i) the intra-combustor interactions between HRR and pressure oscillations, and (ii) the inter-combustor interactions between different thermoacoustic modes. Depending on the coupling and system parameters, these interactions can produce a variety of complex multi-scale dynamics, such as quasiperiodicity, chaos, frequency/phase locking, clustering and chimeras (Juniper & Sujith 2018; Sujith & Unni 2020).

Chimeras are spatiotemporal patterns in which regions of synchrony (coherence) and asynchrony (incoherence) coexist (Panaggio & Abrams 2015). They were first discovered by Kuramoto & Battogtokh (2002) in a population of non-locally coupled phase oscillators governed by the complex Ginzburg–Landau equation. It was found that under certain conditions, the spatial domain splits into two parts: one populated by mutually synchronized oscillators evolving at a common frequency, and one populated by desynchronized oscillators evolving at distributed frequencies (Kuramoto & Battogtokh 2002). Abrams & Strogatz (2004) later named this hybrid pattern a chimera – after the fire-breathing creature from Greek mythology endowed with body parts from different animals – in order to highlight the simultaneous coexistence of coherent and incoherent oscillators in the same network. Since their discovery, chimeras have attracted considerable research attention, leading to the theoretical prediction of numerous variants (Parastesh et al. 2021). These variants have been classified in different ways, such as on the basis of the spatiotemporal evolution of the coherent and incoherent domains (e.g. breathing chimeras, alternating chimeras and travelling chimeras), the emergence of multiple coherent clusters within the incoherent domain (e.g. multi-headed chimeras), and the amplitude evolution of the oscillators (e.g. amplitude chimeras and chimera death) (Parastesh et al. 2021). Although chimeras were originally defined for ensembles of identical oscillators, this definition has since evolved to include ensembles of non-identical oscillators as well (Halatek & Frey 2018). In the real world, chimeras have been observed experimentally in various chemical, mechanical, optical and electrical systems (Parastesh et al. 2021), and they are thought to hold the key to a better understanding of the mutual interactions within such systems (Panaggio & Abrams 2015).

In thermoacoustics, chimeras were first reported by Mondal, Unni & Sujith (2017) for a population of local HRR oscillators representing the reactive flow field of a bluff-body stabilized turbulent premixed combustor. Since then, similar chimeras have been reported by Pawar et al. (2019) for a swirl-stabilized combustor, and by Hashimoto et al. (2019) for a model rocket combustor. However, although pioneering, those studies (Mondal et al. 2017; Hashimoto et al. 2019; Pawar et al. 2019) focused exclusively on single-combustor systems where each pixel group in a flame image was taken to be an individual oscillator in a large spatially-extended network. To our knowledge, chimeras have yet to be observed in a multi-combustor system where each combustor is taken to be an individual oscillator in a small network.

Small networks are minimal systems typically containing between three and ten coupled oscillators. For such networks, Ashwin & Burylko (2015) recently proposed a state of weak chimera characterized by the coexistence of (i) two or more oscillators in frequency synchronization and (ii) one or more oscillators evolving at frequencies different from that of the synchronized ensemble. Shortly after the study by Ashwin & Burylko (2015), the first experimental evidence of weak chimeras was reported by Wojewoda et al. (2016) for three coupled pendulums, and by Hart et al. (2016) for four optoelectronic oscillators. Since then, weak chimeras have been investigated theoretically and experimentally in various systems, such as Stuart–Landau oscillators (Kemeth, Haugland & Krischer 2018), pendulum-like nodes (Maistrenko et al. 2017), electrochemical oscillators (Bick, Sebek & Kiss 2017) and candles (Manoj et al. 2019). To date, however, weak chimeras have yet to be observed in a thermoacoustic system, where mutually constructive interactions between HRR sources and their surrounding acoustic field can give rise to destructive self-excited flow oscillations (Lieuwen & Yang 2005). Establishing the existence of chimeras in a thermoacoustic system would open up new opportunities for the application of oscillation quenching strategies based on chimera control (Bick & Martens 2015; Parastesh et al. 2021). In this study, we present the first evidence of chimera states – namely a breathing chimera and a weak anti-phase chimera – in a multi-combustor system undergoing thermoacoustic oscillations. Moreover, we investigate both the intra- and inter-combustor interactions, and how they give rise to collective multi-scale dynamics.

1.2. Complex networks in thermoacoustics

Recent years have seen complex networks emerge as a powerful tool for investigating the connectivity patterns in various systems, such as our climate, the internet and the human brain (Strogatz 2001). Complex networks are based on network theory, which uses graphs to represent the elements of a system as nodes, and the interactions between them as links. The overarching goal is to better understand how a network of interacting elements – each of which may have its own individual dynamics and coupling architecture – will act collectively (Boccaletti et al. 2006). A focal point of current research is to relate the topology of complex networks built from a physical system to the underlying spatiotemporal dynamics of that system, with a view to developing improved methods of prediction and control (Newman 2018).

Complex networks have seen various applications in fluid mechanics (Donges et al. 2009; Gao et al. 2013; Taira, Nair & Brunton 2016; Kobayashi et al. 2019b; Murugesan, Zhu & Li 2019; Hachijo et al. 2020; Iacobello, Ridolfi & Scarsoglio 2020). In the past five years, they have been used increasingly to characterize and control thermoacoustic oscillations in combustion systems. For example, Murugesan & Sujith (2015) used a visibility algorithm to construct complex networks from time traces of the unsteady pressure in a turbulent combustor during its transition from combustion noise to thermoacoustic instability. They found that networks associated with combustion noise have a scale-free structure, which becomes replaced by an ordered topology at the onset of thermoacoustic instability. Gotoda et al. (2017) used both weighted and unweighted $\varepsilon$-recurrence networks, alongside modified visibility graphs, to show that a turbulent combustor operating near the flame blowout limit can exhibit a network structure with small-world features, implying a large clustering coefficient but a small average path length. At around the same time, Godavarthi et al. (2017) used $\varepsilon$-recurrence networks to analyse the dynamical transitions occurring in a turbulent combustor and found that the characteristic path length is a reliable precursor of thermoacoustic instability. Later, Godavarthi et al. (2018) used recurrence networks to investigate the coupling between the HRR and pressure oscillations in a turbulent combustor transiting through different dynamical states, such as high-dimensional deterministic chaos (combustion noise), intermittency, and limit cycles (thermoacoustic instability). More recently, Guan et al. (2019c) and Guan, Gupta & Li (2020) used filtered visibility graphs to distinguish between noise-contaminated limit cycles and deterministic chaos in both laminar and turbulent combustors.

The use of visibility and recurrence algorithms is not the only way to construct complex networks. For example, Unni et al. (2018) and Krishnan et al. (2019a) built spatial networks in which the connectivity between two nodes is determined by the Pearson (1895) correlation coefficient between the flow velocities at those two nodes. They used such networks to identify (i) the flow regions driving thermoacoustic instability and (ii) the optimal location at which to apply passive control via micro-jet injection. To investigate the spatiotemporal dynamics of the acoustic power sources in a turbulent combustor, Krishnan et al. (2019b) used weighted time-varying spatial networks in which two nodes are considered connected if and only if the product of the HRR and pressure fluctuations at both nodes is positive. Murayama & Gotoda (2019) combined a phase synchronization parameter and the determinism (connection strength) of cross recurrence plots to build weighted networks, and then used them to explain how thermoacoustic oscillations can be weakened with secondary air injection. Recently, Krishnan et al. (2021) constructed time-varying weighted spatial turbulent networks based on the Biot–Savart law, and showed that thermoacoustic instability can be suppressed by targeting the primary network hubs (i.e. the fluid elements with high vorticity) with steady air jets. In summary, the topology and scaling properties of complex networks formed from measurements of combustion systems can be analysed to reveal hidden spatiotemporal patterns in the reactive flow field and to guide control strategies, among other applications. In this study, we use complex networks to investigate the directionality of the intra- and inter-combustor interactions occurring in a network of ring-coupled thermoacoustic oscillators.

1.3. Machine learning in thermoacoustics

For over two decades, machine learning has been used in combustion science to extract actionable information and insight from data (Kalogirou 2003). It has recently attracted even greater interest owing to a confluence of advancements in data collection and storage, computational hardware, and machine learning algorithms (Brunton, Noack & Koumoutsakos 2020). Machine learning algorithms can be classified into three broad groups (supervised, semi-supervised and unsupervised) depending on the degree to which the data are labelled. In this study, we focus on unsupervised machine learning because our aim is to discover hidden patterns in unlabelled combustor data (HRR and pressure signals), rather than to build models for prediction and classification based on new data. Two common applications of unsupervised machine learning are dimensionality reduction and cluster analysis. Dimensionality reduction tools, such as proper orthogonal decomposition (Berkooz, Holmes & Lumley 1993) and dynamic mode decomposition (Schmid 2010; Schmid et al. 2011), have been used to identify the dominant flow structures in various thermoacoustic systems (Mariappan, Sujith & Schmid 2015; Noiray & Denisov 2017; Passarelli et al. 2021; Shoji et al. 2021). By contrast, cluster analysis has attracted less attention in thermoacoustics. Nevertheless, clustering algorithms can be used to partition unlabelled data into distinct and meaningful groups, even without expert knowledge. This could facilitate the discovery of hidden patterns and interactions, as well as the development of reduced-order models for improved physical understanding and control (Brunton et al. 2020). It could also open up new possibilities for data-driven forecasting of the onset of thermoacoustic instability in self-excited combustion systems (Sarkar et al. 2016).

1.4. Contributions of the present study

In this experimental study, we take a complex systems approach to investigating the collective dynamics of a self-excited thermoacoustic system consisting of four turbulent lean-premixed combustors coupled in a ring configuration. Using synchronization metrics and recurrence networks, we examine both the intra-combustor flame–acoustic interactions and the inter-combustor acoustic–acoustic interactions. However, unlike previous studies where the latter interactions were investigated as discrete thermoacoustic modes (Moon et al. 2020b, 2021), here we treat each combustor as an individual self-excited thermoacoustic oscillator and explore how multiple such oscillators can interact in a ring-coupled network to form various synchronous and asynchronous patterns, such as a breathing chimera and a weak anti-phase chimera. We then use recurrence network analysis to identify the dominant direction of the bidirectional coupling (i) between the pressure signals in each pair of directly/indirectly coupled oscillators, and (ii) between the HRR and pressure signals in each individual oscillator. Finally, we use a hybrid machine learning algorithm to perform clustering in a three-dimensional feature space defined by global measures extracted from joint recurrence networks. In this way, we are able to determine which of the two types of interactions (intra- or inter-combustor) plays a more critical role in defining the collective dynamics of the overall system.

This study has two main contributions. First, it shows that even a small network of four ring-coupled combustors can exhibit a wide variety of collective dynamics. These dynamics include intermittent synchronization, clustering, a breathing chimera, and a weak anti-phase chimera, encompassing a diverse mix of order and disorder in the spatiotemporal structure. It is important to recall that these dynamics are not unique to our particular system, but are universal to minimal networks of coupled oscillators (Abrams & Strogatz 2004; Maistrenko et al. 2017; Kemeth et al. 2018). This implies that these dynamics can be modelled – and thus understood, predicted and controlled – using only low-order oscillator equations. Second, this study shows that combining cluster analysis and recurrence network analysis can lead to a versatile tool with which to explore the interactions within and between thermoacoustic oscillators. When combined with chimera control techniques (Bick & Martens 2015; Parastesh et al. 2021), this hybrid machine learning framework could help to guide the design of new ring-coupled combustion systems with reduced susceptibility to thermoacoustic instability.

2. Experimental set-up

Figure 1(a–c) shows the experimental set-up, which is identical to that used by Moon et al. (2020b, 2021). It consists of four cylindrical combustors, each containing a lean-premixed CH$_4$–air flame stabilized in a turbulent swirling injector flow. The four combustors are coupled together in a ring configuration via a full-annular cross-talk section (inner diameter 43 mm) mounted perpendicular to the combustor axis. The length of each combustor is adjustable via a movable piston; in this study, we use two different combustor lengths ($1020$ and $1620$ mm), but the axial position of the cross-talk section is always 20 mm upstream of these lengths ($\xi = 1000$ and $1600$ mm, respectively). The reactant mixture contains gaseous premixed CH$_4$ and air, whose flow rates are metered individually using four thermal mass flow controllers (Teledyne Instruments HFM-D-301 for CH$_4$, and Sierra Instruments FlatTrak 780S for air). We measure the pressure fluctuations in each combustor ($p^{\prime }$) using a piezoelectric transducer (PCB 112A22: sensitivity 14.5 mV kPa$^{-1}$) mounted at the injector plane (figure 1b). This measurement location has been shown by Moon et al. (2020b) to be sufficiently far from the pressure nodes of the system to produce a reliable signal-to-noise ratio. We measure the HRR fluctuations of each flame ($q^{\prime }$) using a photomultiplier tube (Hamamatsu H7732-10) viewing through a bandpass optical filter ($309\pm 5$ nm, OH$^{\ast }$ chemiluminescence). We sample simultaneously the $p^{\prime }$ and $q^{\prime }$ signals from all four combustors at 12 kHz for 4 s on a 16-bit analogue-to-digital converter (TEAC LX-110). This study is guided by two main control parameters: the equivalence ratio of the flame ($\phi$), and the axial position of the cross-talk section ($\xi$), which sets the combustor length as $\xi + 20$ mm. The Reynolds number, defined based on the outer diameter of the mixing section, is held constant at approximately 22 000. Further details on this experimental set-up can be found in Moon et al. (2020b, 2021).

Figure 1. Thermoacoustic system consisting of four turbulent lean-premixed combustors coupled in a ring configuration: (a) isometric view, (b) top cross-sectional view, and (c) end cross-sectional view of the annular cross-talk (XT) section. The dimensions shown are in millimetres. Further details can be found in Moon et al. (2020b, 2021). Panel (d) shows the network architecture, which contains two types of inter-combustor interactions: (i) direct coupling between any two adjacent oscillators, as represented by four pairwise links between $p^{\prime }_{\mathbb {X}}$ and $p^{\prime }_{\mathbb {Y}}$ (solid lines: C1–C2, C2–C3, C4–C3 and C4–C1); and (ii) indirect coupling between any two opposite oscillators, as represented by two pairwise links between $p^{\prime }_{\mathbb {X}}$ and $p^{\prime }_{\mathbb {Y}}$ (dash-dotted lines: C1–C3 and C4–C2). Intra-combustor interactions are captured by the coupling between $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ within each individual oscillator (C1, C2, C3 and C4).

As mentioned earlier, we consider the four combustors as a network of four ring-coupled thermoacoustic oscillators. The architecture of this network is illustrated in figure 1(d). The acoustic interactions between any two oscillators can be decomposed into two classes, depending on the coupling type: (i) direct coupling between any two adjacent oscillators, as represented by four pairwise links between $p^{\prime }_{\mathbb {X}}$ and $p^{\prime }_{\mathbb {Y}}$ (figure 1(d), solid lines: C1–C2, C2–C3, C4–C3 and C4–C1); and (ii) indirect coupling between any two opposite oscillators, as represented by two pairwise links between $p^{\prime }_{\mathbb {X}}$ and $p^{\prime }_{\mathbb {Y}}$ (figure 1(d), dash-dotted lines: C1–C3 and C4–C2). The subscripts $\mathbb {X}$ and $\mathbb {Y}$ are used to index into the four oscillators (1, 2, 3 and 4). As for the intra-combustor flame–acoustic interactions, these are captured by the coupling between $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ within each individual oscillator (C1, C2, C3 and C4).

3. Data analysis via complex systems theory

In this section, we give an overview of the complex systems tools used to analyse the $p^{\prime }$ and $q^{\prime }$ data. For a detailed discussion of these tools, the reader is referred to the books by Webber & Zbilut (2005) and Brunton & Kutz (2019), and to the review papers by Donner et al. (2011) and Zou et al. (2019).

3.1. Kuramoto order parameter

A proven way to identify chimera states is to quantify the phase coherence of all the oscillators in a network. Here we do this using the Kuramoto order parameter (Kuramoto 2003)

(3.1)\begin{equation} R_K( t ) = \frac{1}{N}\left| {\sum_{j = 1}^{N} {{{\rm e}^{{\rm i}\,{\theta _j}(t)}}} }\right|, \end{equation}

where $\theta _j$ is the phase of each oscillator (i.e. the $p^{\prime }$ signal), and $N$ is the total number of oscillators. The oscillators are incoherent when $R_K=0$, but are coherent when $R_K=1$.

3.2. $\varepsilon$-recurrence networks

The recurrence plot (RP) was introduced by Eckmann, Kamphorst & Ruelle (1987) as a graphical means of identifying dynamical states and patterns in time series data using the fundamental property of recurrence. It is a two-dimensional binary bitmap whose elements are defined by the recurrence matrix

(3.2) \begin{equation} \boldsymbol{\mathsf{R}} = \Theta (\varepsilon - \| \boldsymbol{X}_i(d) - \boldsymbol{X}_j(d) \|), \end{equation}

where $\boldsymbol {X}_i(d) = ( {{x_i},{x_{i + \tau }},\ldots,{x_{i + \tau ( {d - 1} )}}} )$ is the $i$th state vector of the system, $\| \cdot \|$ is a norm (e.g. $L_1$ norm, $L_2$ norm or $L_{\infty }$ norm), $\tau$ is the embedding time delay, $d$ is the embedding dimension, $\Theta (\cdot )$ is the Heaviside function, and $\varepsilon$ is the recurrence threshold. The distance matrix, $\| \boldsymbol {X}_i(d) - \boldsymbol{X}_j(d) \|$, is binarized using $\Theta (\cdot )$ such that $\boldsymbol{\mathsf{R}}_{ij} = 1$ when $\| \boldsymbol {X}_i(d) - \boldsymbol{X}_j(d) \| < \varepsilon$, but $\boldsymbol{\mathsf{R}}_{ij} = 0$ otherwise (Marwan et al. 2007).

If the input data are bivariate, as they are in this study (e.g. $p^{\prime }_{\mathbb {X}}$ and $p^{\prime }_{\mathbb {Y}}$, or $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$), then two extensions of the RP are possible: the cross recurrence plot (CRP) and the joint recurrence plot (JRP). The CRP is used to compare the phase trajectories of two different time signals in the same phase space, providing information on the internal coupling between them. The CRP is defined by the cross recurrence matrix

(3.3) \begin{equation} \boldsymbol{\mathsf{CR}} = \Theta (\varepsilon - \| \boldsymbol{X}_i(d) - \boldsymbol{Y}_i(d) \|), \end{equation}

where $\boldsymbol {X}_i(d)$ is the $i$th state vector reconstructed from one time series, and $\boldsymbol{Y}_i(d)$ is the $i$th state vector reconstructed from another time series.

The JRP is used to investigate when two phase trajectories, reconstructed from two time signals, recur simultaneously in their own phase spaces. In other words, it is a measure of the joint probability of two systems recurring at the same time. The JRP is defined by the joint recurrence matrix

(3.4) \begin{equation} \boldsymbol{\mathsf{JR}} = \Theta( {\varepsilon - \| {\boldsymbol {X}_i (d) - \boldsymbol{X}_j (d)} \|} )\cdot \Theta ( {\varepsilon - \| {\boldsymbol {Y}_i (d) - \boldsymbol{Y}_j (d)} \|} ), \end{equation}

where $\boldsymbol {X}_i(d)$ is the $i$th state vector reconstructed from one time series, and $\boldsymbol{Y}_i(d)$ is the $i$th state vector reconstructed from another time series.

Quantitative information on the nonlinear dynamics of a system can be extracted from RPs (including CRPs and JRPs) using recurrence quantification analysis (RQA) (Webber & Zbilut 2005). This involves computing statistical measures (e.g. the recurrence rate, determinism and laminarity) based on the geometric patterns (e.g. diagonal or vertical lines) present in RPs. Once computed, such RQA measures can be used in various ways, such as to distinguish between different types of synchronization and to forecast the onset of critical transitions (Marwan et al. 2007). In combustion science, RQA has been used to detect thermoacoustic instability (Nair, Thampi & Sujith 2014; Hernandez-Rivera et al. 2019; Pagliaroli & Troiani 2020; Pavithran, Unni & Sujith 2021), flame blowout (Unni & Sujith 2016; De et al. 2020) and flame flashback (Christodoulou et al. 2016).

Alternatively, the geometric patterns in RPs can be analysed in the framework of complex networks, via the use of $\varepsilon$-recurrence network analysis. Introduced by Marwan et al. (2009), such an analysis involves casting each state vector, $\boldsymbol {X}_i(d)$, into a vertex (node) of a network. Two vertices are considered connected by an edge if $\boldsymbol{\mathsf{R}}_{ij} = 1$. The binary adjacency matrix of an undirected and unweighted recurrence network is defined as

(3.5)\begin{equation}\boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{R}} - \boldsymbol{\mathsf{I}}_{N}, \end{equation}

where $\boldsymbol{\mathsf{I}}_N$, the $N$-dimensional identity matrix, is used to exclude self-loops. Like RPs, CRPs and JRPs can also be reinterpreted in the framework of complex networks. To facilitate this, Feldhoff et al. (2012) proposed the inter-system recurrence matrix

(3.6) \begin{equation} \boldsymbol{\mathsf{IR}} = \left( {\begin{array}{@{}cc@{}} \boldsymbol{\mathsf{R}}_\mathbb{X}(\varepsilon_{{R}}) & \boldsymbol{\mathsf{CR}}_\mathbb{XY}( \varepsilon_{{CR}}) \\ \boldsymbol{\mathsf{CR}}_\mathbb{YX}( \varepsilon_{{CR}}) & \boldsymbol{\mathsf{R}}_\mathbb{Y}( \varepsilon_{{R}}) \end{array}}\right), \end{equation}

where $\boldsymbol{\mathsf{R}}_\mathbb {X}(\varepsilon _{{R}})$ and $\boldsymbol{\mathsf{R}}_\mathbb {Y}(\varepsilon _{{R}})$ are the recurrence matrices of dynamical systems $\mathbb {X}$ and $\mathbb {Y}$, respectively, and $\boldsymbol{\mathsf{CR}}_\mathbb {XY}( \varepsilon _{{CR}})$ and $\boldsymbol{\mathsf{CR}}_\mathbb {YX}( \varepsilon _{{CR}})$ are the corresponding cross recurrence matrices, with $\boldsymbol{\mathsf{CR}}_\mathbb {YX}^{T} = \boldsymbol{\mathsf{CR}}_\mathbb {XY}$. Moreover, $\varepsilon _{{R}}$ and $\varepsilon _{{CR}}$ are the thresholds for the recurrence matrices and cross recurrence matrices, respectively. The binary adjacency matrix of an inter-system recurrence network can then be defined as

(3.7)\begin{equation} \boldsymbol{\mathsf{IA}} = \boldsymbol{\mathsf{IR}} - \boldsymbol{\mathsf{I}}_{N}. \end{equation}

The joint recurrence network is defined analogously to the recurrence network:

(3.8)\begin{equation} \boldsymbol{\mathsf{JA}} = \boldsymbol{\mathsf{JR}} - \boldsymbol{\mathsf{I}}_{N}. \end{equation}

When constructing $\varepsilon$-recurrence networks, we do not use the full time series of $p^{\prime }$ and $q^{\prime }$ because this would cause the adjacency matrix to be excessively large ($48\,000 \times 48\,000$ elements). Instead, we split the time series into shorter segments using a sliding window of 0.4 s, which is long enough to capture at least $30$ and $100$ cycles of the lowest and highest frequency components, respectively. For comparison, previous analyses involving $\varepsilon$-recurrence networks used a non-overlapping sliding window spanning around 43 cycles of the dominant mode (Godavarthi et al. 2018). To reduce noise, we use a window overlap ratio of 0.5, resulting in 19 shorter time segments. We use a recurrence threshold equal to the fixed recurrence rate ($RR$), following recommendations by Feldhoff et al. (2012) that $RR_{\mathbb {X}\mathbb {Y}} < RR_{\mathbb {X}}, RR_{\mathbb {Y}}$ and $RR_{\mathbb {X}}, RR_{\mathbb {Y}} \leq 0.05$. For the inter-system recurrence networks, we use a threshold of 0.04 for the cross recurrence matrices and 0.05 for the recurrence matrices. For the joint recurrence networks, we use a threshold of 0.05 for the recurrence matrices. A sensitivity analysis (Appendix A) reveals that the results are insensitive to the exact threshold values used, so long as they are within the range suggested by Feldhoff et al. (2012).

We use four measures to quantify the network structure: the cross-transitivity, the global edge density, the global average path length and the global clustering coefficient. We use these specific measures because they have been proven to be able to characterize the mutual synchronization of various coupled systems (Zou et al. 2019), including those in thermoacoustics (Sujith & Unni 2020). Watts & Strogatz (1998) originally proposed the clustering coefficient as a measure of the mean fraction of triangles formed by three different vertices of a network. The clustering coefficient is defined as

(3.9) \begin{equation} \mathcal{C} = \frac{1}{N}\sum_{i = 1}^{N} {\frac{{\sum_{j,k = 1}^{N} {{\boldsymbol{\mathsf{A}}_{ij}}{\boldsymbol{\mathsf{A}}_{{jk}}}{\boldsymbol{\mathsf{A}}_{{ki}}}} }}{{{\boldsymbol{k}_i}( {{\boldsymbol{k}_i} - 1} )}}}, \end{equation}

where $\boldsymbol{\mathsf{A}}_{ij}$ is an element (edge) in the adjacency matrix $\boldsymbol{\mathsf{A}}$, and $\boldsymbol{k}_i = {\sum _{j = 1}^{N} {{\boldsymbol{\mathsf{A}}_{{jk}}}} }$ is the degree centrality, a measure of the number of edges associated with a given vertex $i$. This measure, however, is known to underestimate the actual fraction of triangles in some networks, particularly those dominated by vertices of a small degree, e.g. scale-free networks. To overcome this problem, Barrat & Weigt (2000) proposed the transitivity as an alternative way to quantify the clustering of a network. The transitivity is defined as

(3.10) \begin{equation} {\mathcal{T}} = \frac{\displaystyle{\sum_{i,j,k = 1}^{N} {{\boldsymbol{\mathsf{A}}_{ij}}{\boldsymbol{\mathsf{A}}_{{jk}}}{\boldsymbol{\mathsf{A}}_{{ki}}}} }}{{\displaystyle \sum_{i,j,k = 1}^{N} {{\boldsymbol{\mathsf{A}}_{ij}}{\boldsymbol{\mathsf{A}}_{{ki}}}} }}. \end{equation}

More recently, Feldhoff et al. (2012) proposed the cross-transitivity

(3.11) \begin{equation} {\mathcal{T}^{\mathbb{X}\mathbb{Y}}} = \frac{\displaystyle{\sum_{{i \in \boldsymbol{V}_{\mathbb{X}};\ j,k \in {\boldsymbol{V}_\mathbb{Y}}}} {{\boldsymbol{\mathsf{IA}}_{ij}}{\boldsymbol{\mathsf{IA}}_{{jk}}}{\boldsymbol{\mathsf{IA}}_{{ki}}}} }}{{\displaystyle\sum_{i \in {\boldsymbol{V}_{\mathbb{X}}; j,k \in {\boldsymbol{V}_\mathbb{Y}}}} {{\boldsymbol{\mathsf{IA}}_{ij}}{\boldsymbol{\mathsf{IA}}_{{ki}}}}}}, \end{equation}

where $\boldsymbol{\mathsf{IA}}_{ij}$ is an element (edge) in $\boldsymbol{\mathsf{IA}}$, and $\boldsymbol{V}_\mathbb{X}$ and $\boldsymbol{V}_\mathbb {Y}$ are two disjunct subsets of the whole vertex set $\boldsymbol{V}$ such that $\boldsymbol{V}_{\mathbb {X}} \cap \boldsymbol{V}_{\mathbb {Y}} = \emptyset$ and $\boldsymbol{V}_{\mathbb {X}} \cup \boldsymbol{V}_{\mathbb {Y}} = \boldsymbol{V}$. It is worth noting that ${\mathcal {T}^{\mathbb {X}\mathbb {Y}}}$ and ${\mathcal {T}^{\mathbb {Y}\mathbb {X}}}$ are analogous to $\mathcal {T}$, counting the number of cross-triangles over the number of cross-triples. The cross-transitivity has been used by Godavarthi et al. (2018) to investigate the mutual interactions between $p^{\prime }$ and $q^{\prime }$ in a turbulent premixed combustor undergoing limit-cycle oscillations. This network measure was found to be able to identify the dominant direction of the bidirectional coupling between the pressure and HRR fields, providing physical insight for the development of new detection and control strategies. In § 4.2, we use the cross-transitivity to identify the dominant direction of the bidirectional coupling within and between thermoacoustic oscillators.

The second network measure that we use is the global edge density

(3.12)\begin{equation} {\rho _g} = \frac{1}{{N( {N - 1} )}}\sum_{i = 1}^{N} {{\boldsymbol{k}_i}}, \end{equation}

where $\boldsymbol{k}_i$ is the degree centrality at vertex $i$, and ${\rho _g} \approx RR_{JR}$, i.e. the recurrence rate of the JRP.

The third network measure that we use is the global average path length

(3.13)\begin{equation} \mathcal{L}_g = \frac{1}{{N(N - 1)}}\sum_{i,j = 1}^{N} {{\boldsymbol{\mathsf{l}}_{ij}}}, \end{equation}

where $\boldsymbol{\mathsf{l}}_{ij}$ is the shortest path length between vertices $i$ and $j$. Murugesan & Sujith (2016) have shown that $\mathcal {L}_g$ can capture the changes occurring in a turbulent combustor as it transitions from combustion noise to thermoacoustic instability. During combustion noise, many network nodes exist with many links to their neighbours. During thermoacoustic instability, the network becomes more organized, with nodes having links to only their adjacent neighbours, resulting in a longer average path between nodes. Thus $\mathcal {L}_g$ can be used to detect changes in the system dynamics based on changes in the network connectivity.

The fourth network measure that we use is the global clustering coefficient (see (3.9)), which is denoted by $\mathcal {C}_g$ in this study. Unlike $\mathcal {L}_g$, $\mathcal {C}_g$ is large for a highly connected network (e.g. combustion noise) but is small for a poorly connected network (e.g. thermoacoustic instability). This suggests that $\mathcal {C}_g$ can be used to detect different forms and degrees of mutual synchronization. Further details on these network measures can be found in the review papers by Donner et al. (2011) and Zou et al. (2019).

3.3. Cluster analysis of recurrence network measures

Recently, Kobayashi et al. (2019a) proposed a multi-step hybrid framework combining complex networks and machine learning. First, they constructed ordinal partition transition networks from time traces of the pressure and HRR fluctuations in a swirl-stabilized turbulent combustor. Second, they used principal component analysis to build a two-dimensional feature space containing the first and second components estimated from the probability distribution of the transition patterns. Third, they used a support vector machine with $k$-means clustering to classify the principal components into three groups, each corresponding to a distinct state of the combustor: combustion noise, transition, and thermoacoustic instability. Finally, they monitored the transition state percentage $R_t$, which is defined as the duration ratio between the states of transition and combustion noise. They found that a threshold of $R_t = 0.50$ provides sufficiently early warning for thermoacoustic instability to be suppressed preemptively via secondary air injection at the flame base.

In this study, we propose a similar hybrid framework combining complex networks and machine learning. However, rather than using the framework for early detection of thermoacoustic instability, we use it to determine which of the two types of interactions ($p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ or $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$) contributes more to the collective dynamics of the system. We use a standard Gaussian mixture model (GMM) to perform clustering in a three-dimensional feature space defined by three global measures extracted from joint recurrence networks: $\mathcal {C}_g$, $\rho _g$ and $\mathcal {L}_g$. We examine all three network measures for both $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ and $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions.

Although both the GMM and $k$-means algorithms rely on the use of cluster centres, only the former can account for data covariance, implying that it is more flexible in discovering clusters of different shapes (Press et al. 2007). This flexibility has led to GMM clustering being applied to various observation types, such as human skin tones (Yang & Ahuja 1998), seismic events (Kuyuk et al. 2012) and pulsars (Lee et al. 2012). As § 4.3 will show, that flexibility will also prove to be useful here as the data clusters in our system tend to be of different shapes. In general, the use of GMM clustering involves a model of the form

(3.14)\begin{equation} f\left( {{x_j},\Theta } \right) = \sum_{p = 1}^{k} {{\alpha _p}\,{\mathcal{N}_p} ({{x_j},{\mu _p},{\sigma _p}})}, \end{equation}

where $\Theta$ is a vector of unknown parameters (with mean $\mu _p$ and variance $\sigma _p$), $f(\cdot )$ is the measured probability density function (PDF), and $\mathcal {N}_p(\cdot )$ is the PDF of the sample $x_j$ and is normally distributed. Each $\mathcal {N}_p(\cdot )$ is weighted by $\alpha _p$ such that $\alpha _1+\alpha _2+\cdots +\alpha _k=1$, and $k$ is the number of mixtures (Reynolds 2009). We use an expectation-maximization algorithm to converge iteratively to the maximum likelihood parameters. We standardize the data features by removing their mean and by scaling them to unit variance before clustering.

4. Results and discussion

4.1. Collective dynamics: synchronization and chimeras

We start by examining the collective dynamics of the system under three exemplary operating conditions. These are represented by three networks that differ in their spatial distributions of $\phi$ and $\xi$: network I (§ 4.1.1) and network II (§ 4.1.2) are each populated by identical oscillators, whereas network III (§ 4.1.3) is populated by two different types of oscillators.

4.1.1. Network I: intermittent frequency locking and a breathing chimera

Figure 2 shows the collective dynamics of network I, whose oscillators are identical in equivalence ratio ($\phi _{1,2,3,4} = 0.61$) and cross-talk position ($\xi _{1,2,3,4} = 1600$ mm). We first consider the early stages of the experiment, $0 \leq t \leq 2$ s. In this interval (figure 2(a1), where the time span $0.98 \leq t \leq 1.02$ s is magnified), the time traces of $p^{\prime }$ for all four oscillators show temporally synchronized peaks (red bands) and troughs (blue bands), indicating that the entire network is in a state of global in-phase synchronization (Pikovsky et al. 2003; Balanov et al. 2008). This synchronicity extends to the HRR fluctuations as well (figure 2(a2), $0.98 \leq t \leq 1.02$ s): strong temporal alignment can be observed between the peaks of $p^{\prime }$ and $q^{\prime }$ (red and yellow bands) and between their troughs (blue and green bands). The fact that $p^{\prime }$ and $q^{\prime }$ are evolving in phase (i.e. with a phase difference of less than ${\rm \pi} /2$) implies a positive Rayleigh (1945) integral ($\int p^{\prime } q^{\prime }\,\mathrm {d}t > 0$), providing a physical mechanism by which energy is transferred from the flames to the acoustic modes to sustain the observed self-excited thermoacoustic oscillations (Lieuwen & Yang 2005; Nicoud & Poinsot 2005; Magri, Juniper & Moeck 2020).

Figure 2. Collective dynamics of network I, which exhibits a transition from intermittent frequency locking on a $\mathbb {T}^{3}$ quasiperiodic attractor to a breathing chimera. Shown at the top are time traces of (a1) $p^{\prime }$ and (a2) $p^{\prime }$ and $q^{\prime }$ for each of the four oscillators in the network (C1, C2, C3 and C4); both $p^{\prime }$ and $q^{\prime }$ have been normalized by their respective maximum values from the entire network. Also shown are the spectrograms and PSDs of (b1–4) $p^{\prime }$ and (c1–4) $q^{\prime }$, with the PSDs computed via the algorithm of Welch (1967). Panels (b1, c1), (b2, c2), (b3, c3) and (b4, c4) correspond to oscillators C1, C2, C3 and C4, respectively. The figure also shows (d1, e1) the temporal variation of $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ and $\Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$, alongside (d2, e2) their probability distributions, $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ and $\zeta _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$, where the values of $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ and $\Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ in (d2, e2) are wrapped around the interval $[-{\rm \pi}, {\rm \pi}]$. In the inset of (e1), each curve has been shifted by even integer multiples of ${\rm \pi}$ for clearer visualization. The ( f1) time trace, ( f2) spectrogram and PSD of the Kuramoto order parameter $R_K$ are shown in order to evaluate the phase coherence of the network. In (d1–2, e1–2), the dark- and light-grey regions denote in-phase and anti-phase dynamics, respectively. In (d1, f1–2), the yellow regions denote epochs of global in-phase synchronization.

Although $p^{\prime }$ and $q^{\prime }$ evolve in phase, their dynamics are not simply periodic: the spectrogram and power spectral density (PSD) of both $p^{\prime }$ (figures 2b(1–4)) and $q^{\prime }$ (figures 2c(1–4)) reveal the coexistence of three discrete modes, whose frequencies are incommensurable ($f_1 = 167$, $f_2 = 186$, $f_3 = 201$ Hz). This implies the presence of quasiperiodicity on an ergodic three-dimensional torus attractor (3-torus), otherwise known as $\mathbb {T}^{3}$ quasiperiodicity (Hilborn 2000). Similar $\mathbb {T}^{3}$ quasiperiodic states have been observed by Guan et al. (2019a,b) in a laminar thermoacoustic system consisting of a single Bunsen-flame combustor subjected to external periodic forcing. Other systems in which $\mathbb {T}^{3}$ quasiperiodicity has been observed include Rayleigh–Bénard convection cells (Gollub & Benson 1980), barium–sodium niobate crystals (Martin, Leber & Martienssen 1984) and electrical circuits (Borkowski et al. 2015). In nonlinear dynamical systems, quasiperiodic states on $\mathbb {T}^{n}$ with $n \geq 3$ are known to be unstable to arbitrarily small disturbances (Newhouse, Ruelle & Takens 1978). As a result, they tend to collapse into disorganized states, such as chaos, via a sequence of folding and stretching operations (Hilborn 2000). In our system, we find that the $\mathbb {T}^{3}$ quasiperiodic state is indeed unstable, existing only transiently between asynchronous epochs, before eventually giving way to a breathing chimera in the late stages of the experiment ($t \geq 2$ s), as will be discussed below.

To verify that all four oscillators in the network are in the same $\mathbb {T}^{3}$ quasiperiodic state, we use the Hilbert transform to compute the instantaneous phase difference, $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$, between the pressure signals in each pair of directly/indirectly coupled oscillators ($\mathbb {X}$ and $\mathbb {Y}$). Focusing still on the early stages ($0 \leq t \leq 2$ s), we find that all six oscillator-pair combinations exhibit long epochs in which $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ oscillates within $\pm {\rm \pi}/2$ of even integer multiples of ${\rm \pi}$, indicating in-phase synchronization (figure 2(d1), yellow regions and, in the left inset, dark-grey regions). During these long synchronous epochs, the time-averaged slope of $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$, or $\langle \dot {\Delta \psi }_{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}\rangle$, is zero for all combinations, indicating that all four oscillators are evolving at the same time-averaged frequency, a phenomenon known as frequency locking (Pikovsky et al. 2003). The fact that the instantaneous values of $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ do not remain steady in time, however, implies the absence of phase locking (Pikovsky et al. 2003). The presence of frequency locking without phase locking is a characteristic state in synchronization theory, known as phase trapping. Phase trapping is not unique to thermoacoustic systems (Balusamy et al. 2015; Kashinath, Li & Juniper 2018; Guan et al. 2019a) or even to fluid mechanical systems (Li & Juniper 2013b,a), but can be found in various physical systems, such as solid-state lasers (Thévenin et al. 2011). Its detection here in a ring-coupled thermoacoustic system thus provides further evidence of its universality (Pikovsky et al. 2003; Balanov et al. 2008).

Interspersed between those long epochs of in-phase frequency locking (figure 2(d1), yellow regions) are short epochs of asynchrony in which $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ jumps by even integer multiples of ${\rm \pi}$. These jumps, known as phase slips, are a classic feature of forced or coupled self-excited systems exposed to strong or unbounded noise (Pikovsky et al. 2003). In effect, such noise drives the system from one stable equilibrium point (potential well) to an adjacent one. In our thermoacoustic system, the source of such noise is thought to be turbulence in the underlying reactive flow field. Put together, these observations indicate that during $0 \leq t \leq 2$ s, the system switches intermittently between two distinct regimes: (i) in-phase frequency locking on a $\mathbb {T}^{3}$ quasiperiodic attractor; and (ii) desynchronization due to noise-induced phase slipping. In this study, we refer to this alternating state as intermittent frequency locking. Similar states have been observed previously in other forced or coupled self-excited systems, such as the human brain exposed to external visual stimuli (Pisarchik, Chholak & Hramov 2019). However, it is important to recognize that the intermittent frequency locking observed here differs from the intermittent phase synchronization observed by Pawar et al. (2017): the former involves switching between frequency-locked $\mathbb {T}^{3}$ quasiperiodicity and desynchronization in a small network of four ring-coupled combustors represented by their pressure signals, whereas the latter involves switching between phase-locked periodicity and desynchronization in a single isolated combustor represented by its pressure and HRR signals. Crucially, in our ring-coupled network, all four oscillators are either globally synchronized or globally desynchronized at any given instant (within $0 \leq t \leq 2$ s). In other words, the inter-combustor $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ interactions of the entire network enter the frequency-locked epochs at the same time, and then they switch to the asynchronous epochs also at the same time (figure 2(d1)). Later, however, we will see that such simultaneous global switching does not continue indefinitely, with the network eventually transitioning to a breathing chimera.

Next we examine the probability distribution of $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$, denoted as $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$, where $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ is wrapped around the interval $[-{\rm \pi}, {\rm \pi}]$ (figure 2(d2)). In a continuously desynchronized system (i.e. one without any synchronous epochs), $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ would drift unboundedly in time, causing $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ to be distributed uniformly across all possible values of $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$. By contrast, in a continuously synchronized system (i.e. one without any asynchronous epochs), $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ would be locked to certain values, causing $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ to exhibit dominant peaks. In the special case of intermittent frequency locking, the tail properties of those $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ peaks would depend on the amplitude of phase trapping and on the frequency and duration of phase slipping. In figure 2(d2) (top sub-panel, $0 \leq t \leq 2$ s), we find that nearly all the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ interactions lead to a unimodal $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ distribution with a peak at $-{\rm \pi} /2 < \Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}} < {\rm \pi}/2$. This is consistent with the observed intermittent switching between long epochs of in-phase frequency locking and short epochs of desynchronization associated with phase slipping (figure 2(d1)). However, one specific pair of indirectly coupled oscillators (C1–C3) exhibits a bimodal $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ distribution centred on $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}} \approx 0$. This is caused by the square-like waveform of $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ for C1–C3 (figure 2(d1), see left inset), which shows that C1 leads C3 roughly half of the time, and vice versa the other half. Meanwhile, the other pair of indirectly coupled oscillators (C2–C4) exhibits a unimodal $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ distribution rather than a bimodal one. This indicates an asymmetry in the evolution of $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$, despite the overall system being composed of identical oscillators. We attribute this asymmetry to subtle differences in the operating conditions (e.g. the temperature, velocity and $\phi$ of the reactants) and in the dimensions of the combustors and cross-talk tubes. Indeed, Moon et al. (2019) have shown that two nominally identical uncoupled combustors operated at nominally identical conditions can exhibit different limit-cycle amplitudes in their pressure oscillations. Moreover, Ghirardo et al. (2019) have shown that even slight asymmetries in can-annular combustors can have a noticeable effect on their modal dynamics. Nevertheless, for all six oscillator-pair combinations, every peak in $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ resides well within the in-phase limits of $-{\rm \pi} /2 < \Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}} < {\rm \pi}/2$ (figure 2(d2), top sub-panel, dark-grey region), indicating that the intermittent epochs of phase slipping are neither frequent nor long enough to prevent the four oscillators from evolving in phase with one another on a time-averaged basis.

Moving on to the intra-combustor flame–acoustic interactions, we again use the Hilbert transform to compute the instantaneous phase difference, $\Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$, but this time between the pressure and HRR signals in each individual oscillator ($\mathbb {X}$). We find that $\Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ evolves differently across all four oscillators (figure 2(e1), $0 \leq t \leq 2$ s), despite them being nominally identical in geometry and operating conditions. The cause of this asymmetry is believed to be related to that observed in $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ between C1–C3 and C2–C4. All four oscillators exhibit phase slipping between $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ but to different degrees, with C3 showing the largest and most frequent phase slips. This may be attributed to the relatively strong low-frequency HRR components in C3 (figure 2(c3), see inset). These components are at linear combinations of $f_1$, $f_2$ and $f_3$, implying that they may be due to beating, possibly promoted by the slow temporal variations observed in the mode amplitudes of the HRR signal. These low-frequency components are more prominent in $q^{\prime }_{\mathbb {X}}$ than in $p^{\prime }_{\mathbb {X}}$ (figure 2(c3) versus figure 2(b3)), which would suggest that they are caused by vortical interactions rather than by acoustic interactions. Similar low-frequency components in the HRR signal, arising from an intrinsic hydrodynamic mode in the flame, have been identified by Guan et al. (2019c) as a potential cause of $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ decoupling in a liquid-fuelled combustor with a turbulent diffusion flame. The fact that the oscillator with the lowest degree of phase slipping between $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ (C1) exhibits the weakest low-frequency HRR components lends support to this hypothesis.

Between the epochs of phase slipping, phase locking occurs, as evidenced by $\Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ remaining relatively constant in time (figure 2(e1), see inset, where the data have been shifted by even integer multiples of ${\rm \pi}$ for clearer visualization). This indicates the presence of intermittent phase locking between $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ in each oscillator. The absence of continuous phase locking between $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ is not a violation of the Rayleigh criterion. This is because although the instantaneous values of $\Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ do not always remain at a fixed value, the probability distribution of $\Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$, denoted as $\zeta _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$, still shows a statistical preference for in-phase $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ dynamics ($-{\rm \pi} /2 < \Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}} < {\rm \pi}/2$) in each oscillator (figure 2(e2), top sub-panel). Finally, it is worth recalling that the evolution of $\Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ differs among the four oscillators, with phase slipping in some oscillators often coinciding with phase locking in other oscillators (figure 2(e1)). This stands in stark contrast to the evolution of $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ (figure 2(d1)), demonstrating that global synchronization of $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ across the entire network does not necessarily imply simultaneous phase locking of $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ in each oscillator.

We now consider the late stages of the experiment, $2 \leq t \leq 4$ s. We find that a key difference relative to the early stages ($0 \leq t \leq 2$ s) is that synchronization of $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ no longer occurs simultaneously in every pair of oscillators. Although intermittent frequency locking continues to occur across the entire network, its degree and timing vary locally. For example, in a sample time window $2.3 \leq t \leq 2.34$ s (figure 2(d1), green region), three oscillator-pair combinations (C1–C2, C2–C3, C1–C3) exhibit frequency locking, as evidenced by their $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ oscillating boundedly with $\langle \dot {\Delta \psi }_{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}\rangle = 0$. In the same time window, however, the other three combinations (C3–C4, C4–C1, C2–C4) exhibit desynchronization, as evidenced by their $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ slipping in time. Crucially, in the epochs of frequency locking, $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ shows none of the statistical preference for even integer multiples of ${\rm \pi}$ (in-phase synchronization) seen in the early stages. Taking C2–C3 as an example, we find that some of its frequency-locked epochs show temporal switching between anti-phase and in-phase synchronization (figure 2(d1), top middle inset), some epochs show only in-phase synchronization (figure 2(d1), bottom right inset), and other epochs show only anti-phase synchronization (figure 2(d1), top right inset). Such switching between anti-phase and in-phase synchronization can be seen directly in the pressure traces (figure 2(a1)), where the peaks of C2 are initially aligned with the troughs of C3 ($t \approx 2.3$ s) but then become aligned with the peaks of C3 shortly afterwards ($t \approx 2.34$ s). As expected, this mix of in-phase and anti-phase synchronization causes $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ to be uniformly distributed with no clear peak (figure 2(d2), bottom sub-panel, $2 \leq t \leq 4$ s). From these observations, we can conclude that the system has transitioned from a state in which intermittent frequency locking occurs simultaneously across the entire network ($0 \leq t \leq 2$ s) to one in which intermittent frequency locking occurs at different times and to different degrees in different parts of the network ($2 \leq t \leq 4$ s). As noted in § 1.1, such a hybrid state containing both synchronous and asynchronous spatial domains is called a chimera (Abrams & Strogatz 2004).

To identify the specific type of chimera present in this system, we compute $R_K$ using the instantaneous phase of $p^{\prime }$ in all four oscillators (Kuramoto 2003). On examining the time trace, spectrogram and PSD of $R_K$ (figures 2( f1–2)), we find two different types of behaviour. Early on ($0 \leq t \leq 2$ s), we observe long epochs in which $R_K$ is relatively steady and close to 1 (figures 2( f1–2), yellow regions), indicating that the phasors of all four oscillators are well aligned with one another, consistent with frequency locking occurring in the entire network. Interspersed between those long epochs of $R_K \approx 1$ are short epochs in which $R_K$ fluctuates between 0.5 and 1; these $R_K$ fluctuations arise from rotating phasors associated with phase slipping. Later on ($2 \leq t \leq 4$ s), we find that $R_K$ oscillates continuously with a large amplitude, indicating that the system has transitioned to a non-stationary state in which the positions of the coherent (synchronous) and incoherent (asynchronous) spatial domains vary in time. Such a state is known as a breathing chimera and was first discovered by Abrams et al. (2008) in a network of phase oscillators. An examination of the $R_K$ spectra shows that the breathing frequencies are $f_3-f_1$, $f_3-f_2$ and $f_2-f_1$ (figure 2( f2)), which are linear combinations of the three incommensurable modes identified earlier in the $\mathbb {T}^{3}$ quasiperiodic attractor. It is worth mentioning that a breathing chimera was observed previously by Mondal et al. (2017) in a turbulent combustor during its transition to intermittency. In that study, however, the instantaneous phase was extracted from the HRR field in a single combustor, without any inter-combustor coupling. In the present study, we show that a breathing chimera can also emerge in a small network of four ring-coupled thermoacoustic oscillators, strengthening the universality of this state.

As for the intra-combustor flame–acoustic interactions in the late stages ($2 \leq t \leq 4$ s), we find that the trends identified in the early stages ($0 \leq t \leq 2$ s) are still present. In particular, intermittent phase locking between $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ continues to occur to varying degrees in all four oscillators (figure 2(e1)): C1 exhibits almost continuous phase locking, but C2, C3 and C4 exhibit increasingly pronounced epochs of phase slipping amidst intermittent phase locking. Nevertheless, as in the early stages, all four oscillators spend considerable time with their $\Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ at even integer multiples of ${\rm \pi}$, causing $\zeta _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ to exhibit a dominant peak within the in-phase limits of $-{\rm \pi} /2 < \Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}} < {\rm \pi}/2$ (figure 2(e2), bottom sub-panel, $2 \leq t \leq 4$ s). This implies that across the entire network, $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ are in phase on a time-averaged basis, even though they are not always so on an instantaneous basis. The in-phase relationship between $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ in each oscillator (figure 2(e2), bottom sub-panel) stands in stark contrast to the absence of any statistically preferred phase relationship between $p^{\prime }_{\mathbb {X}}$ and $p^{\prime }_{\mathbb {Y}}$ in the entire network (figure 2(d2), bottom sub-panel). This contrasting behaviour provides further evidence that the phase dynamics of the intra-combustor flame–acoustic interactions does not necessarily dictate that of the inter-combustor acoustic–acoustic interactions.

In summary, we have shown that under certain conditions, a small network of four ring-coupled thermoacoustic oscillators can transition from (i) intermittent frequency locking on a $\mathbb {T}^{3}$ quasiperiodic attractor to (ii) a breathing chimera in which the positions of the synchronous and asynchronous spatial domains vary in time. Furthermore, we have provided evidence showing that the phase dynamics of the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ interactions between coupled oscillators does not necessarily match the phase dynamics of the $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions within those oscillators.

4.1.2. Network II: a two-cluster state of anti-phase synchronization

Figure 3 shows the collective dynamics of network II, whose oscillators are identical in equivalence ratio ($\phi _{1,2,3,4} = 0.65$) and cross-talk position ($\xi _{1,2,3,4} = 1600$ mm). In terms of operating conditions, network II differs from network I (§ 4.1.1) in that it is at a slightly higher equivalence ratio (0.65 versus 0.61). The layout of figure 3 (network II) is identical to that of figure 2 (network I). On examining the time traces of $p^{\prime }$ (figure 3(a1)), we find that unlike network I, network II exhibits no dynamical transitions during the entire sampling interval. Instead, it remains in a stationary state where anti-phase synchronization occurs between any two adjacent oscillators, i.e. between any two directly coupled oscillators (C1–C2, C2–C3, C3–C4, C4–C1). In the literature on can-annular combustors, this is known as a push–pull mode, which Ghirardo, Moeck & Bothien (2020) and von Saldern et al. (2021b), among others, have recently studied via low-order modelling. Owing to the ring-coupled architecture of our system, the fact that anti-phase synchronization occurs between any two directly coupled oscillators (C1–C2, C2–C3, C3–C4, C4–C1) implies that in-phase synchronization occurs between any two indirectly coupled oscillators (C1–C3, C2–C4). The result is a globally synchronous state in which all four oscillators evolve at the same limit-cycle frequency ($f_1 = 210$ Hz, as per the PSDs of both $p^{\prime }$ and $q^{\prime }$; figures 3(b1–4) and 3(c1–4)) but are split into two clusters based on their instantaneous phases: cluster {C1,C3} is in anti-phase synchronization with cluster {C2, C4}. Similar states of clustering have been observed previously in networks of Stuart–Landau oscillators (Manrubia, Mikhailov & Zanette 2004; Premalatha et al. 2018), chemical oscillators (Wang, Kiss & Hudson 2000; Kiss, Zhai & Hudson 2005) and candle-flame oscillators (Manoj et al. 2019; Manoj, Pawar & Sujith 2021). The phase difference $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ for each of the six oscillator-pair combinations remains largely constant in time (figure 3(d1)), indicating that the pressure oscillations in the entire network are not only frequency locked but also phase locked (Pikovsky et al. 2003), with no sign of the phase trapping or slipping seen in network I (figure 2(d1)). As noted earlier, the directly coupled oscillators (C1–C2, C2–C3, C3–C4, C4–C1) are anti-phase synchronized, implying that their $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ is locked to odd integer multiples of ${\rm \pi}$ (figure 3(d1), light-grey regions). By contrast, the indirectly coupled oscillators (C1–C3, C2–C4) are in-phase synchronized, implying that their $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ is locked to even integer multiples of ${\rm \pi}$ (figure 3(d1), dark-grey regions). This mix of anti-phase and in-phase synchronization leads to a clear separation of peaks in the $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ distribution (figure 3(d2)), which is a characteristic feature of clustering based on the instantaneous phase (Manrubia et al. 2004).

Figure 3. The same as in figure 2 but for network II, which exhibits a two-cluster state of anti-phase synchronization on a periodic limit cycle.

In contrast to the mix of anti-phase and in-phase dynamics observed in the inter-combustor $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ interactions, the intra-combustor $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions show only in-phase dynamics. This is apparent in the time traces of $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ (figure 3(a2)) and in the evolution of $\Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ (figures 3(e1–2)). Although all the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ interactions are continuously phase locked (figure 3(d1)), the $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions are continuously phase locked for only C1, and are intermittently phase locked otherwise (figure 3(e1)). This behaviour is similar to that observed in network I (figure 2(e1)), demonstrating that continuous phase locking of $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ across the entire network does not necessarily imply continuous phase locking of $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ in each oscillator. Nevertheless, an examination of $\zeta _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ shows that in all four oscillators, $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ remain in phase on a time-averaged basis, with C1 showing the tallest $\zeta _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ peak (figure 3(e2)), in line with its $\Delta \psi _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ being continuously phase locked (figure 3(e1)).

Figure 3( f1) shows that $R_K$ remains close to zero for the entire sampling interval, which would typically suggest a disordered state of randomly aligned phasors. Here, however, the low values of $R_K$ are caused not by randomly aligned phasors, but by a systematic cancellation of phasors between the two halves of the network: the phasors of cluster {C1,C3} cancel those of cluster {C2,C4}. A similar cancellation of phasors, arising also from anti-phase synchronized clustering, has been reported previously in a large population of strongly coupled relaxation oscillators (Călugăru et al. 2020). The spectrum of $R_K$ is weak and flat, with no noticeable peaks (figure 3(f2)). Recent analysis of a network of coupled Stuart–Landau oscillators by Joseph & Pakrashi (2020) has shown that the conditions favourable to anti-phase synchronization include a small number of oscillators (typically less than 20), low connectivity, weak coupling over distance, and strong symmetry in the network topology. Our thermoacoustic system satisfies all of these conditions: it has only four oscillators (four combustor nodes), each oscillator is directly coupled to only its two adjacent neighbours (i.e. nearest neighbour coupling), the coupling is achieved with a single annular cross-talk section, and the network topology is a symmetric ring. The discovery of anti-phase synchronization in such a network of four ring-coupled thermoacoustic oscillators is thus consistent with the low-order network analysis of Joseph & Pakrashi (2020).

4.1.3. Network III: a weak anti-phase chimera

Figure 4 shows the collective dynamics of network III, whose oscillators are identical in cross-talk position ($\xi _{1,2,3,4} = 1000$ mm) but not in equivalence ratio ($\phi _{1,3} = 0.61$, $\phi _{2,4} = 0.57$). On examining the time traces of $p^{\prime }$ (figure 4(a1)), we find that the network is stationary but divided into two halves: each of the two pairs of indirectly coupled identical oscillators (C1–C3, C2–C4) is in anti-phase synchronization, producing push–pull modes (Moon et al. 2021), but the two pairs are not synchronized with each other. In other words, half of the network (C1–C3) evolves at one frequency ($f_1 = 262$ Hz), while the other half (C2–C4) evolves at a different frequency ($f_2 = 77$ Hz), as can be seen in the PSDs of both $p^{\prime }$ (figures 4(b1–4)) and $q^{\prime }$ (figures 4(c1–4)). In recent experiments, Moon et al. (2020b) observed a similar frequency-based partitioning of the spatial domain and attributed it to mode localization induced by a loss of rotational symmetry in the network.

Figure 4. The same as in figure 2 but for network III, which exhibits a weak anti-phase chimera.

Examining the inter-combustor $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ interactions, we find that $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ for both pairs of indirectly coupled identical oscillators (C1–C3, C2–C4) is approximately constant in time, with values hovering near odd integer multiples of ${\rm \pi}$, interrupted by occasional phase slips (figure 4(d1)). In the $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ curve (figure 4(d2)), this gives rise to peaks at $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}} = \pm {\rm \pi}$, indicating that both pairs of indirectly coupled identical oscillators (C1–C3, C2–C4) are undergoing intermittent phase locking in an anti-phase manner (Pikovsky et al. 2003). By contrast, in all four pairs of directly coupled non-identical oscillators (C1–C2, C2–C3, C3–C4, C4–C1), $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ drifts unboundedly in time (figure 4(d1)), indicating desynchronization associated with phase drifting. This is caused by the oscillator frequency alternating between $f_1$ and $f_2$ around the network. Thus $\zeta _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ is uniformly distributed across all possible values of $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ (figure 4(d2)). Collectively, these observations indicate that network III is in a stationary state of weak anti-phase chimera. As alluded to in § 1.1, a weak anti-phase chimera is a special type of chimera in which two or more oscillators in a network evolve in anti-phase synchronization, while at least one other oscillator evolves at a frequency different from that of the synchronized ensemble (Ashwin & Burylko 2015). Weak anti-phase chimeras have been predicted theoretically by Maistrenko et al. (2017) in a minimal network of three identical pendulum-like nodes. Here, we present experimental evidence showing that a weak anti-phase chimera can emerge in a thermoacoustic system, namely a small network of four non-identical thermoacoustic oscillators coupled in a ring configuration.

Turning now to the intra-combustor $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions, we find that intermittent phase locking occurs in all four oscillators, but to different degrees (figure 4(e1)): the two oscillators at $\phi _{1,3} = 0.61$ (C1 and C3) exhibit only occasional phase slips, whereas the two oscillators at $\phi _{2,4} = 0.57$ (C2 and C4) exhibit frequent phase slips. In the former two oscillators (C1 and C3), the peak in $\zeta _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ sits near the boundary between in-phase and anti-phase dynamics (figure 4(e2)), indicating that the Rayleigh criterion is barely satisfied in C1 and C3. By contrast, in the latter two oscillators (C2 and C4), the peak in $\zeta _{p^{\prime }_{\mathbb {X}}q^{\prime }_{\mathbb {X}}}$ sits well within the in-phase limits (figure 4(e2)), indicating that the Rayleigh criterion is satisfied in C2 and C4, with appreciable energy being transferred from the flames to the acoustic field. This increase in energy transfer is believed to be the physical cause of the higher pressure amplitudes observed in C2 and C4 relative to C1 and C3 (figure 4(a1)).

Figure 4( f1) shows that $R_K$ for network III is higher than that for network II (figure 3(f1)), indicating greater phase coherence. Although phasor cancellation occurs in both pairs of anti-phase synchronized oscillators (C1–C3, C2–C4), the process is not perfect owing to phase slips and minor fluctuations in $\Delta \psi _{p^{\prime }_{\mathbb {X}}p^{\prime }_{\mathbb {Y}}}$ (figure 4(d1)). Given that the phasors in C1–C3 rotate at $f_1$ while those in C2–C4 rotate at $f_2$, a new spectral component emerges at their linear combination ($f_1 - f_2$), as can be seen in the spectrogram and PSD of $R_K$ (figure 4(f2)).

4.2. Identifying the dominant coupling direction via recurrence networks

Using recurrence network analysis (§ 3.2), we compute the cross-transitivity to identify the dominant direction of the bidirectional coupling (i) between the pressure signals in each pair of directly/indirectly coupled oscillators, $\mathcal {T}_{p^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {Y}}}$ and $\mathcal {T}_{p^{\prime }_{\mathbb {Y}} p^{\prime }_{\mathbb {X}}}$ (figure 5a), and (ii) between the pressure and HRR signals in each individual oscillator, $\mathcal {T}_{p^{\prime }_{\mathbb {X}} q^{\prime }_{\mathbb {X}}}$ and $\mathcal {T}_{q^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {X}}}$ (figure 5b). Identifying the degree of asymmetry in these interactions can reveal physical insight into the network coupling architecture, thus aiding the development of new control strategies.

Figure 5. Recurrence network analysis of the three networks in § 4.1: (a) cross-transitivity of $p^{\prime }_{\mathbb {X}}$ and $p^{\prime }_{\mathbb {Y}}$ for each pair of directly/indirectly coupled oscillators; (b) cross-transitivity of $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ for each individual oscillator. The vertical marker bars represent the standard deviation. Also shown are network diagrams illustrating the coupling architecture in (c) network I, (d) network II, and (e) network III.

For any two coupled oscillators $\mathbb {X}$ and $\mathbb {Y}$, if $\mathcal {T}_{p^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {Y}}} < \mathcal {T}_{p^{\prime }_{\mathbb {Y}} p^{\prime }_{\mathbb {X}}}$, then the phase trajectory of $\mathbb {Y}$ is pulled towards that of $\mathbb {X}$, implying that their bidirectional coupling is biased in the direction $\mathbb {X} \rightarrow \mathbb {Y}$ (Feldhoff et al. 2012). However, it should be emphasized that this does not necessarily imply that oscillator $\mathbb {X}$ has unidirectional authority over oscillator $\mathbb {Y}$. Instead, it simply means that the influence of oscillator $\mathbb {X}$ on oscillator $\mathbb {Y}$ is greater than the influence of oscillator $\mathbb {Y}$ on oscillator $\mathbb {X}$. In other words, the two oscillators are still interacting bidirectionally, albeit through asymmetric coupling.

Before discussing the results, we note that although network I is non-stationary (§ 4.1.1), its cross-transitivity data remain statistically stationary over the entire sampling interval (not shown, for brevity). On this basis, we time-average the cross-transitivity data for network I in the same way as we do for networks II and III, which are inherently stationary. Even so, figures 5(a,b) reveal variable data scatter in some cases (the vertical marker bars represent the standard deviation), making it non-trivial to compare different datasets. Here, we consider $\mathcal {T}_{p^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {Y}}}$ to be statistically different from $\mathcal {T}_{p^{\prime }_{\mathbb {Y}} p^{\prime }_{\mathbb {X}}}$ if $|\mathcal {T}_{p^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {Y}}} - \mathcal {T}_{p^{\prime }_{\mathbb {Y}} p^{\prime }_{\mathbb {X}}}| / \max (\mathcal {T}_{p^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {Y}}}, \mathcal {T}_{p^{\prime }_{\mathbb {Y}} p^{\prime }_{\mathbb {X}}}) \geq 0.1$. We use an analogous criterion to compare $\mathcal {T}_{p^{\prime }_{\mathbb {X}} q^{\prime }_{\mathbb {X}}}$ and $\mathcal {T}_{q^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {X}}}$.

In networks I and II, we find that $\mathcal {T}_{p^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {Y}}} \approx \mathcal {T}_{p^{\prime }_{\mathbb {Y}} p^{\prime }_{\mathbb {X}}}$ for every pair of directly/indirectly coupled oscillators (figure 5a). This shows that the inter-combustor acoustic–acoustic coupling in both networks is globally symmetric, with no dominant direction of authority between any two nodes, regardless of whether they are directly or indirectly coupled. Such globally symmetric coupling is due to each network being populated by identical oscillators. As for the intra-combustor flame–acoustic coupling, we find that, somewhat counterintuitively, $\mathcal {T}_{p^{\prime }_{\mathbb {X}} q^{\prime }_{\mathbb {X}}} \approx \mathcal {T}_{q^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {X}}}$ for only C1 and C2 in network I, and for only C1 and C4 in network II (figure 5b), indicating that symmetric coupling between the HRR and pressure fields exists in only half of each network, despite identical oscillators all around. Furthermore, we find that $\mathcal {T}_{p^{\prime }_{\mathbb {X}} q^{\prime }_{\mathbb {X}}} > \mathcal {T}_{q^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {X}}}$ for the other half of each network, indicating asymmetric coupling in which the HRR field exerts a greater influence on the pressure field than vice versa. This particular form of asymmetric $q^{\prime }_{\mathbb {X}} \rightarrow p^{\prime }_{\mathbb {X}}$ coupling has also been observed by Godavarthi et al. (2018) in a turbulent premixed combustor without any inter-combustor coupling. The inter- and intra-combustor interactions identified in networks I and II are summarized graphically in figures 5(c) and 5(d), respectively.

In network III, we find that $\mathcal {T}_{p^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {Y}}} < \mathcal {T}_{p^{\prime }_{\mathbb {Y}} p^{\prime }_{\mathbb {X}}}$ for every pair of directly coupled oscillators (figure 5(a), C1–C2, C2–C3, C3–C4, C4–C1). This indicates asymmetric coupling, with a dominant direction that forms a clockwise loop around the perimeter of the network (figure 5(e), $\text {C1}\rightarrow \text {C2}\rightarrow \text {C3}\rightarrow \text {C4}\rightarrow \text {C1}$). As for the two pairs of indirectly coupled oscillators (C1–C3, C2–C4), we find that $\mathcal {T}_{p^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {Y}}} \approx \mathcal {T}_{p^{\prime }_{\mathbb {Y}} p^{\prime }_{\mathbb {X}}}$, indicating symmetric coupling. This is due to the coexistence of two opposing coupling paths. Between C1 and C3, these paths are $\text {C1}\rightarrow \text {C2}\rightarrow \text {C3}$ and $\text {C3}\rightarrow \text {C4}\rightarrow \text {C1}$; analogous paths exist between C2 and C4. Ghirardo et al. (2019, 2020) showed that mode localization in coupled can-annular combustors can be caused by asymmetric local perturbations, such as those arising from variations in the can geometry or the flame response. We speculate that mode localization in network III (§ 4.1.3) arises from the asymmetric spatial distribution of equivalence ratio ($\phi _{1,3} = 0.61$, $\phi _{2,4} = 0.57$), which splits the network in half with non-identical oscillators and thus biases the inter-combustor coupling in the circumferential direction. By contrast, because both network I (§ 4.1.1) and network II (§ 4.1.2) contain identical oscillators, they exhibit globally symmetric coupling with no mode localization. Moving on to the intra-combustor flame–acoustic coupling (figure 5b), we find that the half of the network at $\phi _{1,3} = 0.61$ (C1 and C3) has symmetric coupling, whereas the other half, at $\phi _{2,4} = 0.57$ (C2 and C4), has asymmetric coupling. As is the case for networks I and II, the asymmetric flame–acoustic coupling is biased such that the HRR field exerts a greater influence on the pressure field than vice versa (figure 5e). A practical implication of this asymmetric $q^{\prime }_{\mathbb {X}} \rightarrow p^{\prime }_{\mathbb {X}}$ coupling is that if the goal is to control thermoacoustic oscillations by disrupting the intra-combustor $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions, then modifying the flame response may be more effective than modifying the acoustic field. However, as the next subsection will show, if the inter-combustor $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ interactions could be modified as well, then doing so may provide a more effective means of control.

In all three networks, we find that the coupling between any two signals ($p^{\prime }$ or $q^{\prime }$) becomes more symmetric as their degree of phase locking increases. For example, in network I (§ 4.1.1), $p^{\prime }_{\mathbb {X}}$ and $q^{\prime }_{\mathbb {X}}$ experience stronger phase locking in C1 than in C4. As a result, the difference between $\mathcal {T}_{p^{\prime }_{\mathbb {X}} q^{\prime }_{\mathbb {X}}}$ and $\mathcal {T}_{q^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {X}}}$ is smaller in C1 than in C4 (figure 5b). This trend holds not just for the intra-combustor flame–acoustic interactions, but also for the inter-combustor acoustic–acoustic interactions. For example, in network III (§ 4.1.3), the pressure signals are more synchronized in the indirectly coupled oscillators (C1–C3, C2–C4) than they are in the directly coupled oscillators (C1–C2, C2–C3, C3–C4, C4–C1). As a result, the difference between $\mathcal {T}_{p^{\prime }_{\mathbb {X}} p^{\prime }_{\mathbb {Y}}}$ and $\mathcal {T}_{p^{\prime }_{\mathbb {Y}} p^{\prime }_{\mathbb {X}}}$ is smaller in the former group than in the latter group (figure 5a). Unravelling the relationships between the synchronization behaviour and coupling architecture of such networks will be crucial for understanding, predicting and controlling the dynamics of ring-coupled thermoacoustic systems.

4.3. Cluster analysis of recurrence network measures

We use the hybrid machine learning framework proposed in § 3.3 to discover patterns in the network structure. The aim is to determine whether the collective dynamics is dominated by the inter-combustor $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ interactions or by the intra-combustor $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions. We perform a GMM cluster analysis in a three-dimensional feature space defined by three global measures extracted from joint recurrence networks: $\mathcal {C}_g$, $\rho _g$ and $\mathcal {L}_g$. The results of this analysis are shown in figure 6 for both $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ and $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ objects. For both types of objects, we find that the optimal number of clusters is $n_c = 4$ (Appendix B).

Figure 6. GMM cluster analysis in a feature space defined by three global measures extracted from joint recurrence networks: the global clustering coefficient $\mathcal {C}_g$, the global edge density $\rho _g$, and the global average path length $\mathcal {L}_g$. Panel (a) shows $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ objects, and (b) shows $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ objects. The bottom row shows a graphical representation of the cluster distribution for each of the three networks from § 4.1.

Starting with the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ data (figure 6a), we find that each of the four clusters ($\mathcal {G}1$, $\mathcal {G}2$, $\mathcal {G}3$, $\mathcal {G}4$) is homogeneous, containing objects from only a single network. For example, $\mathcal {G}1$ and $\mathcal {G}2$ contain objects from networks I and II, respectively; these two clusters are aligned parallel to the main diagonal in the $\rho _g$–$\mathcal {C}_g$ plane (figure 6(a), see inset). Similarly, both $\mathcal {G}3$ and $\mathcal {G}4$ contain objects from network III. Specifically, $\mathcal {G}3$ captures only the directly coupled oscillators (C1–C2, C2–C3, C3–C4, C4–C1); these objects are scattered in the $\mathcal {L}_g$–$\mathcal {C}_g$ plane but are concentrated at a particular $\rho _g$ slice. By contrast, $\mathcal {G}4$ captures only the indirectly coupled oscillators (C1–C3, C2–C4); these objects are scattered relatively evenly across the entire feature space. It is worth noting that although both the indirectly coupled oscillators in network III (figure 6(a), $\mathcal {G}4$) and the directly coupled oscillators in network II (figure 6(a), $\mathcal {G}2$) exhibit anti-phase synchronization on a periodic limit cycle, only the former undergoes intermittent phase locking. As a result, the cluster structures of $\mathcal {G}4$ and $\mathcal {G}2$ are markedly different from each other. This shows that even subtle differences in the synchronization dynamics can be identified via changes in the recurrence network measures.

Unlike the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ data (figure 6a), the $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ data (figure 6b) are grouped into some heterogeneous clusters, i.e. clusters containing objects from more than one network. For example, $\mathcal {G}2$ contains objects from both networks I and III, while $\mathcal {G}3$ contains objects from both networks II and III. When viewed from above (figure 6(b), see inset), $\mathcal {G}1$, $\mathcal {G}2$ and $\mathcal {G}3$ are aligned along the cross-diagonals of the $\rho _g$–$\mathcal {C}_g$ plane, while $\mathcal {G}4$ is concentrated at the core.

In summary, by performing a GMM cluster analysis of three global measures extracted from joint recurrence networks, we have shown that the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ objects form only homogeneous clusters. By contrast, the $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ objects form both homogeneous and heterogeneous clusters. This suggests that the network features arising from the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ interactions are more distinctive than those arising from the $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions. From this, we can conclude that the inter-combustor acoustic–acoustic interactions are more important than the intra-combustor flame–acoustic interactions in defining the collective dynamics of the system.

5. Conclusions

In this experimental study, we have taken a complex systems approach to investigating the collective dynamics of four turbulent lean-premixed combustors coupled in a ring configuration. We treated each combustor as an individual self-excited thermoacoustic oscillator and explored how multiple such oscillators can interact in a network to form various synchronous and asynchronous patterns, many of which have not been observed previously in ring-coupled combustors. Specifically, we considered the intra-combustor flame–acoustic interactions as a process of mutual coupling between the HRR fluctuations ($q^{\prime }_{\mathbb {X}}$) and pressure fluctuations ($p^{\prime }_{\mathbb {X}}$) in each individual combustor. Similarly, we considered the inter-combustor acoustic–acoustic interactions as a process of mutual synchronization between the pressure fluctuations ($p^{\prime }_{\mathbb {X}}$ and $p^{\prime }_{\mathbb {Y}}$) in directly/indirectly coupled combustors.

Using synchronization metrics derived from the Hilbert transform and the Kuramoto order parameter, we found a wide range of complex multi-scale dynamics, depending on the spatial distribution of the equivalence ratio and cross-talk position. These dynamics include (network I) a transition from intermittent frequency locking on a $\mathbb {T}^{3}$ quasiperiodic attractor to a breathing chimera in which synchronous and asynchronous spatial domains move around in time, (network II) a two-cluster state of anti-phase synchronization on a periodic limit cycle, and (network III) a weak anti-phase chimera. In both the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ and $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions, we found evidence of phase slipping or drifting occurring between epochs of phase locking or trapping. However, we found that the phase dynamics of $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ does not always match that of $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$, with phase or frequency locking in the former often coinciding with phase slipping or drifting in the latter. Finally, we found that identical oscillators in a network (networks I and II) evolve globally at the same time-averaged frequencies, but that non-identical oscillators in a network (network III) exhibit mode localization, with different parts of the network evolving at different frequencies. This mix of frequencies is characteristic of a weak chimera. Along with the breathing chimera found in network I, this constitutes the first evidence of chimera states in a minimal network of coupled thermoacoustic oscillators.

We then used recurrence network analysis, namely the cross-transitivity, to identify the dominant direction of the bidirectional coupling in $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ and $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$. We found that the mode localization observed in network III arises from its non-identical oscillators biasing the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ coupling in the circumferential direction. By contrast, because both networks I and II contain only identical oscillators, they exhibit globally symmetric $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ coupling with no mode localization. Counterintuitively, even with identical oscillators, both networks I and II exhibit a combination of symmetric and asymmetric $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ coupling. The asymmetric $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ coupling is always biased such that the HRR field exerts a greater influence on the pressure field than vice versa ($q^{\prime }_{\mathbb {X}} \rightarrow p^{\prime }_{\mathbb {X}}$). This suggests that modifying the flame response may be more effective than modifying the acoustic field, if the aim is to control thermoacoustic oscillations by disrupting the intra-combustor $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions. However, if the inter-combustor $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ interactions can be modified as well, then doing so may offer a more effective means of control, as our cluster analysis shows. Thus identifying the degree of asymmetry in the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ and $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions can provide physical insight into the network coupling architecture and help to guide the design of new control strategies.

Finally, we proposed a hybrid framework, combining unsupervised machine learning and complex networks, to discover hidden patterns in the network structure. We performed a GMM cluster analysis in a three-dimensional feature space defined by three global measures extracted from joint recurrence networks: the global clustering coefficient ($\mathcal {C}_g$), the global edge density ($\rho _g$), and the global average path length ($\mathcal {L}_g$). We found that the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ objects form only homogeneous clusters, whereas the $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ objects form both homogeneous and heterogeneous clusters. This indicates that, compared with the intra-combustor flame–acoustic interactions, the inter-combustor acoustic–acoustic interactions are more distinctive and thus play a more critical role in defining the collective dynamics of the system.

The implications of this study are twofold. First, we have shown that even a small network of four ring-coupled combustors can exhibit a wide variety of collective dynamics, such as intermittent frequency locking, in-phase/anti-phase synchronization, clustering, $\mathbb {T}^{3}$ quasiperiodicity, a breathing chimera and a weak anti-phase chimera. These dynamics encompass a broad mix of order and disorder, producing spatiotemporal patterns where coherent and incoherent domains coexist. Despite their complexity, however, these dynamics are known to be universal to minimal networks of coupled oscillators (Abrams & Strogatz 2004; Maistrenko et al. 2017; Kemeth et al. 2018). Our findings thus lend support to the use of canonical low-order models (e.g. the Van der Pol and/or Stuart–Landau oscillators) to understand, predict and control the thermoacoustics of coupled combustion systems. Second, we have shown that cluster analysis can be combined with recurrence network analysis to create a powerful tool with which to explore the $p^{\prime }_{\mathbb {X}}$–$p^{\prime }_{\mathbb {Y}}$ and $p^{\prime }_{\mathbb {X}}$–$q^{\prime }_{\mathbb {X}}$ interactions occurring between and within thermoacoustic oscillators, respectively. When combined with chimera control techniques (Bick & Martens 2015; Parastesh et al. 2021), this hybrid machine learning framework could offer valuable clues as to how certain chimera states (e.g. chimera death, Zakharova, Kapeller & Schöll 2014) can be reached. In turn, this would open up new pathways to controlling self-excited thermoacoustic oscillations in coupled combustion systems, potentially improving the performance and service life of gas turbines.