{'fraction_non_alphanumeric': 0.03257743565895268, 'fraction_numerical': 0.01090401772055854, 'mean_word_length': 4.877859300820026, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multilayer network analysis to study complex inter-subsystem interactions in a turbulent thermoacoustic system

Thermoacoustic systems are complex systems where the interactions between the hydrodynamic, acoustic and heat release rate fluctuations lead to diverse dynamics such as chaos, intermittency, and ordered dynamics. Such complex interactions cause catastrophically high-amplitude acoustic pressure oscillations and the emergence of order in the spatio-temporal dynamics, referred to as thermoacoustic instability. In this work, we use multilayer networks to study the topology of inter-subsystem interactions between the hydrodynamic, acoustic and combustion subsystems in a bluff-body stabilized turbulent dump combustor during diverse dynamical states. We construct a two-layered network where the layers represent the thermoacoustic power and vorticity fields. The inter-layer links are determined using cross-variable short-window correlations between vorticity and thermoacoustic power fluctuations. We distinguish the topology of inter-layer networks during different dynamical states using network properties such as assortativity and link-rank distributions. During chaotic dynamics, interactions between subsystems are non-localized and spread throughout the flow field of the combustor. During the state of thermoacoustic instability (order), the inter-layer network topology reveals dense inter-layer connections between regions of coherent thermoacoustic power generation and regions where large vortices are shed periodically. Furthermore, we show that such dense inter-layer connections emerge in localized pockets in the recirculation zone in the dump plane of the combustor much prior to the onset of order, that is, during the state of intermittency. We identify these regions of intense inter-subsystem interactions as hubs of the multilayer network, and conjecture that targeted attack at such hubs disrupts inter-subsystem interactions and thus helps mitigate thermoacoustic instability.


Introduction
Thermoacoustic systems, such as those in combustion chambers of gas-turbine engines or a rocket combustors, can exhibit diverse dynamical states owing to nonlinear interactions between the hydrodynamic field, the acoustic pressure oscillations and the heat released due to combustion (Lieuwen & Yang 2005;Sujith & Unni 2020;Sujith & Pawar 2021). Thus, thermoacoustic systems are essentially complex systems. A system is considered to be complex if its dynamics is determined by the interaction of its constituent subsystems. Dynamics of complex systems cannot be fully comprehended by studying any of the subsystems in isolation from the other subsystems as is done in conventional reductionist approaches. Many problems in fluid mechanics pose similar challenges, such as turbulence which has vortical interactions at multiple timescales, multi-phase flows containing interactions between the different phases (such as fluid and particles), and the climate system where interactions between wind flow, cloud cover, solar radiation and the earth's topology together determine the dynamics. In such systems, it is essential to analyse the spatio-temporal interactions both within each subsystem and between subsystems.
The pattern of intra and inter-subsystem interactions and the resulting dynamics are sensitive to the control parameters of the system. In a thermoacoustic system, the structure of nonlinear interactions between the acoustic, hydrodynamic and combustion subsystems varies with the inlet mass flow rates of fuel and air (or Reynolds number), as well as the ratio of these fuel and air mass flow rates (i.e. equivalence ratio). The stable state of operation of a combustor is often referred to as combustion noise. This state is characterized by chaotic spatio-temporal dynamics such as incoherent heat release and irregular shedding of small vortices (George et al. 2018) and low-amplitude aperiodic fluctuations in the acoustic pressure signal (Gotoda et al. 2011;Nair et al. 2014). Moreover, the acoustic pressure time series obtained during the state of combustion noise delineates high-dimensional chaos (Tony et al. 2015). With change in the control parameters, the complex feedback interactions between the subsystems of a thermoacoustic system cause more acoustic driving than acoustic damping in the system. As a result, the system exhibits thermoacoustic instability which is characterized by ordered dynamics including large-amplitude periodic oscillations in the acoustic pressure signal (Lieuwen 2002), coherent heat release and periodic shedding of large coherent structures (vortices) (George et al. 2018;Poinsot et al. 1987). The transition from combustion noise (chaos) to thermoacoustic instability (order) occurs via the route of intermittency which is characterized by epochs of periodicity interspersed by epochs of chaotic dynamics (Nair et al. 2014;Unni & Sujith 2017;Kheirkhah et al. 2017;Ebi et al. 2018).
The occurrence of thermoacoustic instability is catastrophic and can lead to humongous loss of resources and revenue. Several studies have been invested on understanding the phenomenological causes and devising effective means of control of thermoacoustic instability (Candel 2002;Lieuwen & Yang 2005;Sujith & Pawar 2021). For example, researchers use different approaches such as dynamical systems and bifurcation theory (Ananthkrishnan et al. 1998;Lieuwen 2002;Gotoda et al. 2014;Laera et al. 2017), framework of synchronization (Pawar et al. 2017), measuring flame response to the acoustic field perturbations using impulse response functions (Polifke 2020) or flame transfer/ describing functions (Dowling 1997;Noiray et al. 2008;Hemchandra et al. 2011;Schuller et al. 2020), and complex networks (Murugesan & Sujith 2015;Gotoda et al. 2017;Sujith & Unni 2020). Recently, the transition from combustion noise to thermoacoustic instability has been viewed as a transition from chaos to order in the turbulent reacting flow field of a thermoacoustic system (Sujith & Pawar 2021). However, to fully comprehend the emergence of order in the spatiotemporal dynamics of a thermoacoustic system, we must investigate the pattern of nonlinear interactions between multiple constituent subsystems during various dynamical states. For example, we must study how the emergence of periodicity in the acoustic field may be related to the emergence of large vortices in the combustor; or how the spatio-temporal distribution of heat release rate fluctuations are dependent on the acoustic and the hydrodynamic field.
Early seminal works by Poinsot et al. (1987), Sterling & Zukoski (1987) and Schadow & Gutmark (1992) discuss the role of coherent structures in exciting periodic dynamics in the combustion chamber, thus leading to thermoacoustic instability. In these early attempts, techniques such as spark-Schlieren imaging and C 2 radiation maps were used to investigate the simultaneous evolution of large coherent structures and the motion of regions of intense heat release, respectively. These studies describe the complete cycle of formation, growth and decay of large coherent structures (vortices) and the burning of fuel-air mixture carried by these vortices. They show that the reaction zone trails behind the vortices.
Vortices cause mixing of the jet of fresh mixture and the burnt hot products in the flow field. Further, these large coherent structures promote bulk mixing, while the braids of such vortices (regions of high strain rate) ensure fine-scale mixing between hot products and the unburnt reactants (Poinsot et al. 1987;Schadow & Gutmark 1992). Also, vortices may impinge onto the walls of the combustor and break down causing fine-scale mixing of fuel and air. When sufficient mixing occurs, all the fuel-air mixture contained in the large vortex undergoes combustion at once, thus causing large and sudden heat release in the combustion chamber. Since these large vortices are shed periodically, we obtain periodic spikes in the heat release rate fluctuations. In locations where the heat release rate fluctuations are in phase with the acoustic pressure, acoustic driving becomes dominant eventually leading to thermoacoustic instability (Sterling & Zukoski 1987).
Thus, these studies theorized that there are inter-dependencies between the emergence of acoustic modes, coherent structures and heat release rate fluctuations during the state of thermoacoustic instability. Subsequently, the role of acoustic modes on the evolution of coherent structures during the state of thermoacoustic instability has been investigated using different techniques such as phase averaging and proper orthogonal decomposition (Lacarelle et al. 2010;Oberleithner et al. 2011;Tammisola & Juniper 2016), and dynamic mode decomposition (Schmid et al. 2011;Premchand et al. 2019). Moreover, active control techniques have been developed by modifying vortex shedding frequencies and suppressing the coherence of large vortices that aid in exciting thermoacoustic instability (McManus et al. 1993;Paschereit et al. 1999).
A quantitative analysis of such interactions between the vorticity field and an acoustically driven field of heat release rate fluctuations remains an interesting area of research for the thermoacoustic community. Moreover, many interesting questions about such intersubsystem interactions remain unanswered. For example, how does the pattern of intersubsystem interactions differ during various dynamical states such as chaos (combustion noise), intermittency and order (thermoacoustic instability)? Are there localised spatial pockets where intense inter-subsystem interactions occur during or prior to the onset of thermoacoustic instability? Answering these questions can aid in understanding the emergence of order amidst chaos in the spatio-temporal dynamics of a thermoacoustic system.
Complex systems, such as thermoacoustic systems, are best analysed using complex networks. A complex network is a set of nodes connected by links, where these links are defined by the relation between the nodes. Various network construction techniques can be used to construct networks from temporal or spatio-temporal data to infer the dynamics of that system (Gao et al. 2017;Iacobello et al. 2021). Networks have been used extensively to characterize climate systems (Tsonis et al. 2006), fluid mechanics problems involving turbulence (Iacobello et al. 2018;Taira et al. 2016), and thermoacoustic systems (Sujith & Unni 2020).
In thermoacoustic literature, different types of complex networks have been used to investigate the transition between diverse dynamical states in a combustor. For example, complex networks derived from time series of acoustic pressure signal have been used to analyse the emergence of periodicity in the acoustic field of a thermoacoustic system. Complex networks based on visibility algorithm (Murugesan & Sujith 2015) revealed that the scale-free nature of networks derived from the acoustic pressure dynamics of a thermoacoustic system is lost as order emerges amidst chaos in the dynamics of the system. Further, recurrence analysis of acoustic pressure oscillations captures transitions of the dynamical states of a turbulent combustor and offers early warning signals for thermoacoustic instability (Godavarthi et al. 2017). Also, joint recurrence networks were used to analyse the interdependence between the time series of acoustic pressure and heat release rate fluctuations during different dynamical states in a turbulent combustor (Godavarthi et al. 2018). Variation of measures derived from visibility network and recurrence networks were used to detect the onset of thermoacoustic instability as well as flame-blowout (Gotoda et al. 2017).
Another network construction method based on correlation between phase space cycles derived from acoustic pressure time series highlighted the psuedo-periodicity and highdimensional nature of acoustic pressure dynamics during the state of thermoacoustic instability Okuno et al. (2015). Further, Tandon & Sujith (2021) used cycle networks to show that as order emerges from chaos in the acoustic pressure dynamics of a turbulent combustor, the phase space topology morphs monotonically from a set of several highly unstable periodic orbits during the occurrence of chaotic dynamics (combustion noise), to a combination of moderately stable and unstable periodic orbits during intermittency, and finally to a stable limit cycle attractor during ordered dynamics (thermoacoustic instability). Also, Kobayashi et al. (2019) have used ordinal partition transition networks based on the synchronization of acoustic pressure and heat release rate fluctuations to identify the onset of thermoacoustic instability. Moreover, Aoki et al. (2020) used the network entropy of ordinal partition transition networks built from acoustic pressure signal to identify the transition between chaotic and periodic dynamics during the occurrence of intermittency.
Further, transitions in spatio-temporal dynamics of thermoacoustic systems have also been investigated using complex networks. Networks are constructed by treating spatial points in the flow field as nodes. The links between these nodes are determined via some relation between the variation of a certain flow variable, such as velocity or heat release rate fluctuations, at the location of these nodes. Here, the relation itself may be defined via linear or nonlinear measures such as correlation, mutual information or Biot-Savart law for induced velocity. For instance, Unni et al. (2018) used correlation between velocity fluctuations at different locations to construct unweighted networks. Using these timeaveraged spatial networks, the authors identify the different regions where such correlations are most significant during distinct dynamical states. Similarly, Krishnan et al. (2019a) constructed weighted networks from velocity correlations and suggested passive control strategies to mitigate thermoacoustic instability by targeting regions of high node strengths in the velocity-based networks. Recently, Hashimoto et al. (2019) built turbulence networks where the links are determined by the velocity induced at a certain location due to vorticity at another location using the Biot-Savart law, as proposed by Taira et al. (2016). Eventually, Kurosaka et al. (2021) and Krishnan et al. (2021) showed that attacking the spatial locations of primary hubs identified by the turbulence network using localized secondary jet injections can mitigate the occurrence of thermoacoustic instability. Also, complex networks were constructed from the spatial field of thermoacoustic power generated in a turbulent combustor to analyse the feedback between the heat released due to combustion and the acoustic field (Krishnan et al. 2019b;Shima et al. 2021).
Note that, in the literature hitherto, networks are constructed using the information of a single dynamic variable derived from the turbulent reacting flow field of a combustor. Thus, such a network construction represents the interactions within one subsystem of the combustor on what is called a 'single-layer' network. For example, a network built from the velocity fluctuations represents the spatio-temporal interactions within the hydrodynamic field alone and can be called the hydrodynamic layer network. Such single-subsystem (or single-layer) networks do not encode inter-subsystem interactions. However, the diverse dynamics observed in a thermoacoustic system are clearly emergent phenomena arising from the interactions between the hydrodynamic, acoustic and combustion fields (Sujith & Unni 2020). Further, the evolution of each subsystem also depends on the evolution of the other subsystems. Thermoacoustic systems have multiple levels of complexity, that is, complex interactions within subsystems as well as complex interactions across subsystems. One natural solution to further decipher the higher-level complexity of such systems is to investigate the dynamics through a network of single-layer networks, that is, a multilayer network.
A multilayer network consists of layers representing different subsystems of a complex system (Bianconi 2018). Links established between nodes within a layer represent interactions within that subsystem, while links between nodes in different layers represent the interactions between different subsystems. Using multilayer network analysis, one can identify nodes which have strong influence in both layers or in facilitating inter-layer interactions and so on. Multilayer networks have been discussed extensively and well-developed mathematically in recent years (Boccaletti et al. 2014;Kivelä et al. 2014;De Domenico et al. 2013 The application of multilayer networks to fluid mechanics can offer great insight and help reveal deeper complexities in fluid mechanics problems such as those related with multiple phases, interaction of numerous vortical structures, and multiple timescales. Specific to thermoacoustic systems, the approach involving multilayer networks is the key to answering the questions we have raised previously. In this work, we identify the patterns of interaction between the hydrodynamic, acoustic and heat release rate fields during various dynamical states in a turbulent bluff-body stabilized dump combustor. We construct a two-layered network where one layer represents the hydrodynamic subsystem and the other layer represents the dynamics of the acoustically coupled heat release rate field in the combustor. In order to decipher the structure of inter-subsystem interactions during various dynamical states, we analyse how the links are distributed between different spatial locations across the two layers of the combustor. To do so, we introduce the concept of inter-layer network assortativity and inter-layer link-rank distribution. Our results corroborate the theory put forward by Poinsot et al. (1987) and Schadow & Gutmark (1992). That is, we find dense inter-layer connections between regions of coherent vortex shedding in the hydrodynamic layer and the regions of high heat release rate fluctuations in the other layer. Thus, we show that the interaction between coherent structures and the acoustic and combustion subsystems is indeed significant during the state of thermoacoustic instability. Moreover, our analysis reveals that significant localized intersubsystem interactions emerge during the state of intermittency much prior to the onset of thermoacoustic instability. By examining the topology of the inter-layer connections, we identify regions where the interaction between the hydrodynamic field and the acoustically coupled heat release rate fluctuations is predominant. In this manner, we attempt to quantify the role of coherent structures in the emergence of order in the spatio-temporal dynamics of a turbulent thermoacoustic system.
The details of the experimental setup and data processing are discussed in Section 2, while the steps involved in the inter-layer network construction are presented in Section 3. Network properties such as assortativity and link-rank distribution for multilayer networks are discussed in Section 4. The results of multilayer network analysis on different dynamical states in a turbulent bluff-body stabilized dump combustor are discussed in Section 5. We conclude the study and discuss the scope for future work in Section 6.

Experiments in a turbulent thermoacoustic system
Experiments were performed in a turbulent combustor with a backward-facing step and a bluff-body to stabilize the flame for combustion. A schematic of the experimental setup is shown in figure 1. The cross-section of the combustion chamber is 90 mm × 90 mm while its length is 1100 mm. The bluff-body used as a fame holding device is a circular disk of 10 mm thickness and 47 mm diameter located at a distance of 4.5 mm from the inlet. The bluff-body is mounted on a central shaft having a diameter of 16 mm. The fuel is provided through a 16 mm diameter central shaft and ejected through four circumferential holes of 1 mm diameter. Fuel (liquid petroleum gas) and air are mixed in the burner before being inlet into the combustion chamber and it is ignited using a spark plug. This fuel-air mixture is injected at a distance of 110 mm upstream of the dump plane. The variation of fuel and air flow rates is measured in standard litres per minute (SLPM) and is monitored using mass flow controllers (Alicat Scientific, MCR Series). In the experiments, the mass flow rate of fuel is maintained at = 30 ± 0.44 SLPM while the mass flow rate of air ( ) is varied quasi-statically between 480 ± 7.84 to 780 ± 10.24 SLPM. As a result, the Reynolds number ( ) and the equivalence ratio ( ) of the inlet fuel-air mixture also vary with subsequent uncertainties of ±6% and ±0.02, respectively. The Reynolds number and the equivalence ratio ( ) are calculated based on the procedure given in Nair et al. (2014), taking into account the mixture viscosity (Wilke 1950). The range of Reynolds numbers and equivalence ratio reported in the experiment are 18000 − 24000 and 0.99 − 0.60.
The temporal variation of acoustic pressure is measured using a piezoelectric pressure transducer (PCB 103B02 having a sensitivity of 223.4 mV/kPa and measurement uncertainty is ±0.15 Pa) The time series of global heat release rate is acquired using a photomultiplier tube (Hamamatsu H10722-01) fitted with an OH * filter (narrow bandwidth filter centred at 308 nm and 12 nm full width at half maximum). Both the time series are acquired at a sampling rate of 10 kHz for 3 s. Simultaneous high-speed particle image velocimetry (PIV) and OH * chemiluminescence imaging are performed to acquire velocity field measurements and the spatio-temporal variation of unsteady heat release rate, respectively. A high-speed CMOS camera (Phantom-v 12.1) outfitted with an OH * filter (308 nm, 12 nm FWHM) is used to acquire the OH * chemiluminescence images of the flame at a frame rate of 2000 Hz for 1.258 s for each equivalence ratio. The physical region captured has dimensions 67 mm × 44 mm which is mapped to 800 × 600 pixels resolution.
For PIV, the region of interest is illuminated by a double pulsed Nd:YLF laser (527 nm, 25 mJ at 1 kHz frequency, Photonics laser). A high-speed CMOS camera (Photron FASTCAM SA4) equipped with a ZEISS 100 mm lens along with a short bandpass optical filter centred at 527 nm (12 nm FWHM) is used to capture the light scattered (Mie scattering) by the seeding particles (1 TiO 2 particles). We use a fluidized bed seeder for an adequate supply and uniform distribution of TiO2 in the flow. The detailed procedure for optimizing the seeding control is described in Raghunathan et al. (2020). The initial resolution of the Mie scattering images acquired using the Photron camera is 1024 × 1024 pixels at a frequency of 2000 Hz covering a physical region of dimensions 58.5 mm × 58.5 mm. To obtain the velocity and vorticity field data, a cross-correlation algorithm with multiple pass grid refining technique is used on the Mie scattering images. The size of the interrogation window and the maximum overlap is chosen such that, for all flow rates, the resolution of the velocity field remains the same. After numerous post-processing steps discussed by George et al. (2018), the size of the region of interest becomes 41 × 41 and 63 × 63 spatial points for data acquired at inlet air mass flow rates greater than and lesser than 600 SLPM, respectively. The velocity field data acquired after post-processing of the PIV data is obtained at a frequency of 1000 Hz for 1.36 s. For further details of the analysis tools used to extract and process the experimental data, the reader may refer to George et al. (2018).
In our experiments, at low inlet mass flow rates and high equivalence ratio, the combustor operates in stable mode. This state is called combustion noise and is characterized by lowamplitude aperiodic acoustic field fluctuations. As we increase the inlet mass flow rate of air (and subsequently decrease the equivalence ratio), the combustor exhibits the state of intermittency characterized by short burst of high-amplitude periodic oscillations amidst low-amplitude aperiodic fluctuations in the acoustic field. With further increase in the inlet mass flow rate of air (decrease in the equivalence ratio), the combustor operates in unstable mode, referred to as thermoacoustic instability characterized by high-amplitude periodic oscillations in the acoustic field.

Constructing the inter-layer network from experimental data
We construct a two-layered network with inter-layer connections to decipher the spatial pattern of interactions between the hydrodynamic, acoustic and heat release rate fields during different dynamical states in the combustion chamber. Here, a layer of the network represents a subsystem and each layer consists of the same set of nodes representing the same spatial points in the flow field. We consider two different physical quantities derived from the experimental data, which are vorticity ( ) and ′ ′ fluctuations. Here, ′ ′ is the temporal product of acoustic pressure ( ′ ) and the heat release rate fluctuations ( ′ ) at any specific location. In one layer, the nodes represent the vorticity fluctuations at different spatial locations in the combustor and this layer is called the vorticity layer or the -layer. In the other layer, a node represents the time series of ′ ′ at the location corresponding to that node. We refer to this layer as the thermoacoustic power layer or the ′ ′ -layer. Next, an inter-layer link is established between a node in the ′ ′ -layer and a node in the -layer if the time series of ′ ′ and at the respective locations are correlated over short epochs (as explained later). The resulting network is hereafter referred to as inter-layer network.
In order to construct such a network, we make the spatial and temporal resolution of the and ′ ′ fields equal. We compare the analysis windows of the cameras used for chemiluminescence imaging and PIV, and then select the common region with dimensions 58.5 mm × 58.5 mm with a common spatial resolution of 39 × 39 pixels. Data derived from chemiluminescence imaging and the processed PIV data for inlet air flow rates less than 600 SLPM have higher resolution than that for air flow rates more than 600 SLPM, which is reduced to 39 × 39 pixels using bicubic interpolation method (Zhang et al. 2011). As a result, we have the same number of nodes at precisely the same spatial locations in both the layers. In multilayer network terminology, nodes representing the same entity (here, spatial location) in different layers are referred to as replica nodes (Bianconi 2018). We extract the acoustic pressure, vorticity and heat release rate fluctuations at a common temporal resolution of 1000 Hz frequency for 1.25 . After excluding the spatial points that lie in the region of the bluff-body location, we obtain = 1037 nodes containing non-trivial data in each layer.
To calculate the correlation between two time series we use an adjusted correlation function ℛ , which is a variant of the conventional Pearson correlation measure. Note that, if two time series do not have similar probability distributions, the range of their correlation may be distorted from the theoretically expected range of [−1, 1] and may be smaller than this ideal range (Ratner 2009). To correct for the range, we find the maximum and minimum possible correlation between the two time series. Here, the data points of both the input time series are rearranged in ascending order of their values and a correlation is found, called the positive rematch which is essentially the maximum possible correlation. Similarly, if the data of one of the input time series is rearranged in ascending and that of the other in descending order while finding correlation, we obtain the minimum possible correlation called the negative rematch. The actual range of correlation of the two input time series is between the negative and the positive rematch.
Note that, we need to compare the correlation values obtained from correlating the time series of different physical variables for various pairs of locations. To make a fair comparison these correlation values must lie in the same range that is [-1,1]. In order to correct the range of correlation and bring it to the ideal range, we must normalise the Pearson correlation value, as suggested by Ratner (2009). If the conventional Pearson correlation between the two time series is positive (or negative), we normalise it with the positive (or negative) rematch; and the correlation value thus obtained is called the 'adjusted' correlation.
Also, the timescale at which the dynamics ensues in the hydrodynamic field is not necessarily the same as that in the acoustic and the combustion subsystems, especially during the states of combustion noise and intermittency. Due to distinct timescales involved, the calculation of correlation between ′ ′ and fluctuations is not straightforward. The correlations within the hydrodynamic field may be sustained only for short epochs (a few cycles of the dominant acoustic mode) owing to turbulent fluctuations (Lieuwen 2002). Turbulence also introduces variations in the amplitude envelope of the acoustic pressure signal (Lieuwen 2002), and such cycle-to-cycle variability reduces the time window for strong correlations. Thus, high values of cross-variable correlations will be sustained only for short epochs.
Hence, we find cross-variable adjusted correlation for short epochs between the time series of ′ ′ and for different pairs of locations to determine the inter-layer connections. We call this correlation measure as short-window adjusted ( ) correlation. Figure 2 summarizes the procedure of determining -correlation measure. Consider two time series, and obtained at nodes and in layers 1 and 2 of the inter-layer network, respectively. We divide each time series into short non-overlapping windows, where each window has a size of seconds. This window is equivalent to integral number of cycles of the dominant mode of acoustic oscillations ( ′ ) that occur at = 140 Hz frequency. Thus, cycles are equivalent to = / s. The same window size is used across all dynamical states. The data in the th window of the time series is then written as = [( − 1) : ]. Note that, the two time series and correspond to different variables representing dynamics of different subsystems. Such variables not only ensue at different time scales, but also influence each other after a certain time delay, probably, owing to the shedding and convection of coherent structures in the flow field (Lieuwen & Yang 2005;George et al. 2018). To incorporate such time delay, we calculate the correlation between the ℎ window of and a range of windows from − to + of the time series , where, determines the maximum time-lag allowed for correlating two variables. Similarly, correlation is calculated between all windows of the two time series and the maximum (maximized over and ) of such correlation values is defined as the -correlation. The definition of short-window adjusted correlation ( ) is given by equation 3.1.
where, is the number of windows a time series is divided into, and and are the window indices. Here, the optimum window size is chosen by finding the number of cycles ( ) in a window for which the mean of -correlations of time series of all pairs of nodes in distinct layers is maximum. We vary between three to ten and find that such optimum window size occurs when = 4. Here, = 10 corresponding to a maximum time-lag of 0.2 s is chosen to allow us to capture delays associated with hydrodynamic convection timescales. Note that, the time delay associated with correlation between time series at any two spatial locations is variable and depends on the distance between the spatial points and their location with respect to the combustor walls or the bluff-body. The use of -correlation helps extract the maximum of the short-window correlation between the two time series corresponding to such unknown variable time-lags for any pair of locations.
In summary, we use short-window adjusted correlation to determine the weights of inter-layer edges between nodes and in layers 1 and 2 respectively, thus accounting for time-lagged correlation between different spatial locations. In order to track only the strong correlations, we set a threshold of ℎ = 0.75 for the -correlations. The resulting thresholded -correlation matrix is the adjacency matrix of the inter-layer network. Thus, Note that, the adjacency matrix is not symmetric. Particularly, in our system, the interlayer link with weight gives the correlation between the ′ ′ time series at node in the ′ ′ -layer ( 1 ) and the vorticity time series at node in the -layer ( 2 ). On the other hand, gives the correlation between the vorticity time series at node in the -layer and the ′ ′ time series at node in the ′ ′ -layer. Thus, ≠ , and hence is an asymmetric matrix.

Tools for multilayer network analysis
In this section, we discuss the basic properties of multilayer networks used in the current work to analyse the network topology, namely, inter-layer node strengths, inter-layer network assortativity and inter-layer link-rank distribution. We discuss the utility of these tools in identifying the topology of the inter-layer network, and hence, the patterns of inter-subsystem interactions during different dynamical states in a thermoacoustic system.

Inter-layer node strength
The inter-layer node strength of a node in a particular layer is defined as the sum of the weights of its inter-layer links (Bianconi 2018). In a two-layered network, a node in layer 1 has a node strength defined by equation 4.1.
where is the total number of nodes in each layer and 1 , 2 is the weight of the link between node in layer 1 and node in layer 2 . Clearly, the node strength of node in layer 2 (defined by equation 4.2) is different than that of node in layer 1 , since the adjacency matrix is not symmetric (refer Section 3).
The inter-layer node strength of a node signifies the tendency of that node to form inter-layer connections with the nodes in the other layer. Thus, the inter-layer node strength of a node in the ′ ′ -layer signifies if the thermoacoustic power generated at that location is strongly correlated with the hydrodynamic activity at several other locations in the combustion chamber. Similarly, the inter-layer node strength of a node in the -layer signifies if the vorticity dynamics at that location is correlated to the thermoacoustic power generated at different locations in the combustion chamber. Also, nodes with very high node strengths are referred to as hubs. The spatial distribution of such node strengths in the combustor geometry can thus help us identify hubs of the inter-layer network which are essentially spatial pockets having significant inter-subsystem influence.
Also, along with such spatial distribution of node strengths, we must decipher the entire topology of inter-layer connections. For example, we are interested to know if there are any significant hub-to-hub connections across layers; or are the inter-layer connections of a hub in one layer spread over a large spatial region in the other layer. To do so, we borrow the concept of assortativity from single-layer networks and extend it to inter-layer networks.
Inter-layer network assortativity Assortativity measures the tendency of a node to connect to other nodes with similar characteristics in the network (Barabási 2013). Particularly, degree-assortativity examines if high-degree nodes tend to connect to similar high-degree nodes or not (Newman 2002;Barabási 2013). Note that, degree is the total number of connections that a node has in a network. In single-layer networks, a standard measure of assortativity is the degree correlation function ( ) which captures the relation between the degree of connected nodes (Pastor-Satorras et al. 2001;Vázquez et al. 2002;Barabási 2013). The degree correlation function is defined by equation 4.3.
where ( ′ | ) is the conditional probability of finding a node of degree ′ in the neighborhood of a node with degree . Thus, ( ) quantifies the average degree of the neighbors of all those nodes which have degree (Barabási 2013). The variation of ( ) with helps decipher the assortativity of a network. In an assortative network, similar-degree nodes have high tendency to connect. In such a network, higher the degree of a node, higher will be the average degree of its neighbors. Thus, the degree correlation function ( ) is expected to increase with . On the other hand, the degree correlation function decreases with for disassortative networks.
Degree correlations can be used to study the network topology of multilayer networks as well. Recently, various definitions of multilayer network assortativity have been proposed based on degree correlations. Reis et al. (2014) discuss the importance of degree correlations of the inter-layer connections on the structure of directed interconnected networks. Also, de Arruda et al. (2016) proposed the extension of degree correlations based on tensor notation of multilayer networks and analysed the effect of assortativity on epidemic spreading in multilayer networks. Here, we introduce a similar extension of the definition of degree correlations for a two-layered inter-layer network consisting of inter-layer connections alone.
A multilayer network can be represented by what is called a supra-adjacency matrix which is a block matrix built from the adjacency matrices of the intra-layer and inter-layer networks (Bianconi 2018). The diagonal blocks of the supra-adjacency matrix represent the intra-layer connections, which are void in our case. The off-diagonal blocks represent the inter-layer connections in the multilayer network. We construct an unweighted supra-adjacency matrix from the inter-layer network by operating the Heaviside function Θ on the adjacency matrix , as defined in equation 4.4.
The supra-adjacency matrix therefore represents a network of 2 nodes with nodes in each layer. In our system, we note that the replica nodes in each layer essentially represent the same spatial point in the combustor. However, we treat them as distinct nodes in this formalism as they possess distinct information arising from different subsystems. The degree of each node is calculated as the number of inter-layer links based on this unweighted interlayer network. Then, we calculate the degree correlation function for the matrix using the conventional definition as given by equation 4.3. Here, is referred to as the inter-layer degree correlation function. We study the variation of with the inter-layer degree which is derived as the sum of number of inter-layer links in the network represented by .
If ( ) increases with , the inter-layer connectivity is assortative, i.e., a node with high degree in one layer tends to connect only to nodes with high degrees in the other layer. Similarly, if ( ) decreases with , the inter-layer network is disassortative and a hub in one layer tends to have many connections to low-degree nodes in the other layer. Further, such behavior of the degree correlation function can be approximated as a power law, as in equation 4.5 (Barabási 2013).

( ) ∼
( 4.5) where is the correlation exponent. Real-world networks can exhibit different values of for different range of node degrees . Such mixed assortativity has been observed earlier, for example, in citation networks (Barabási 2013). The variation of the degree correlation function with degree provides a global picture of the tendency of assortative connections in a network. However, the local tendency of nodes to connect (or not) to similar-degree nodes is not evident (Noldus & Van Mieghem 2015). The degree-correlation method, while useful, has certain drawbacks and cannot help picture the complete topology of connections in a network (Noldus & Van Mieghem 2015; Litvak & Van Der Hofstad 2013). Further, Xu et al. (2009) explain that the presence of a few hubs or 'super-rich' nodes in a network can induce 'false' disassortativity and distort the degree correlations significantly. By its very definition, a hub is connected to several nodes in the network. Then, many low-degree nodes are obviously connected to the hub, thus inducing negative degree correlations. Moreover, there are fewer hubs than low-degree nodes in a typical network. So, even if a hub tends to connect to other hubs in the network, such hub-to-hub connections are lesser in number compared to connections between a hub and low-degree nodes. In other words, even though a hub may have few local assortative connections (hub-to-hub), a hub-node induces disassortativity in the degree correlation function. Therefore, it is essential to specifically examine if there are hub-to-hub connections in a network irrespective of its global assortativity.
Note that, an assortative inter-layer network derived from a thermoacoustic system (as described in Section 3) would imply that inter-subsystem interactions occur predominantly between the regions with similar node strengths in both the layers. In contrast, a disassortative inter-layer network would imply that the influence of inter-subsystem activity in the hubregions of one layer extends to large spatial regime comprising the low-degree nodes in the other layer. However, we must also examine the connections of hubs specifically. A hub in one layer may have several links to nodes with low degrees in the other layer causing disassortativity, but this hub can also have significant connections to hubs in the other layer. In other words, despite disassortative degree correlations, an inter-layer network can have significant hub-to-hub links between layers. It is therefore necessary to characterize the distribution of links between hub-regions and low-degree regions across the two layers to infer the extent of inter-subsystem activity in our combustor. Thus, we propose the use of link-rank distributions between nodes in different layers as described in the following sub-section.

Inter-layer link-rank distribution
The hubs in a network largely determine the topology of the network and the subsequent interpretations. These large degree nodes in a network are said to be 'rich nodes'. If the rich nodes are interconnected, they form a group referred to as the rich-club ( In an inter-layer network, a rich club would signify the presence of dense hub-to-hub interlayer connections between layers. Particularly, the existence of a rich-club in disassortative networks would signify that there are many hub-to-hub connections in the network as well as significant connections between hubs and low-degree nodes. Thus, we must characterize the distribution of inter-layer links between high as well as low degree nodes in both the layers. To do so, we construct an inter-layer link-rank distribution as described in figure 3. First, we rank the nodes in each layer according to their inter-layer degrees. The highest degree node in a layer has rank one and so on. Note that, the ranking of nodes in each layer is independent of the nodes in the other layer. We thus assign ranks ′ ′ to nodes in the ′ ′ -layer and ranks to nodes in the -layer. Next, we form bins of ranks in each layer where each bin classifies a range of ranks. For example, bin one constitutes nodes having the first 10% ranks, the second bin constitutes nodes with the next 10% ranks and so on. Thus, nodes in each layer are distributed into ten bins according to their ranks. Next, we examine what fraction of the total inter-layer links connect nodes in the ℎ bin of the ′ ′ -layer and nodes in the ℎ bin of the -layer (where and vary from one to ten). This gives us a distribution of inter-layer links classified according to the ranks (degrees) of the nodes connected by these links. We call such a distribution as an inter-layer link-rank distribution. Note that, this inter-layer link-rank distribution proposed here is inspired from the node-node link distribution suggested by Zhou & Mondragón (2004) for single-layer networks.
We can categorize the inter-layer links according to the ranks of the nodes they connect, as shown in figure 3. A link that connects a node with rank ′ ′ in the ′ ′ -layer and a node with rank in the -layer is referred as ′ ′ − to − link. If the link connects low ranking nodes (hub) in the ′ ′ -layer to high ranking nodes (low-degree) in the -layer, then such a link will be called a low-to-high type link, as represented by a green box in the schematic inter-layer link-rank distribution in figure 3. Similarly, there can be low-to-low (white box in figure 3), high-to-low (blue box in figure 3) and high-to-high (black box in figure 3) type of links in an inter-layer network.
The link-rank distribution is very helpful in deciphering several features of the inter-layer network topology. Firstly, if the percentage of low-to-low type of links is high in an interlayer network, we can infer that there are dense hub-to-hub inter-layer connections in the network. Also, if the percentage of low-to-high links is very significant, we infer that hubs in the first layer connect to regions of low node strengths in the other layer. Secondly, we can also compare the different types of links to infer the density of inter-layer connections between different regions of the combustor. Using inter-layer degree correlation function and inter-layer link-rank distribution, we examine the topology of the inter-layer networks derived during various dynamical states in a turbulent combustor and the results are presented in Section 5.

Results
We use two-layered inter-layer networks to analyse the inter-subsystem interactions during diverse dynamical states in a turbulent bluff-body stabilized combustor. To assess the raw flow field, we plot the normalised time series of acoustic pressure and heat release rate fluctuations (row-I) along with the spatio-temporal evolution of the ′ ′ field (row-II) and vorticity (row-III) in figures 4, 7, 8, 11 during the occurrence of combustion noise (chaos), aperiodic and periodic epochs of intermittency and thermoacoustic instability (order), respectively. The ′ and ′ time series are normalized with the corresponding maximum value during each dynamical state to facilitate comparison across different dynamical states. Next, using the tools defined in Section 4, we analyse how the topology of the inter-layer connections changes as the dynamics of the combustor transitions from chaos to order via the route of intermittency. To aid in understanding the topology of the inter-layer network, we also plot a multilayer visualization of the inter-layer connections for specific nodes in both the layers. We consider one sample node in the recirculation zone in the dump plane upstream of the bluff-body and another in the wake downstream of the bluff-body. The two-layered visualization is derived from the inter-layer network constructed from the experimental data during distinct dynamical states.

Combustion noise
During the state of combustion noise, the combustor exhibits low-amplitude chaotic fluctuations in the acoustic pressure ( ′ ) and global heat release rate ( ′ ) signals, as shown in figure 4-(I). Further, the value of ′ ′ is small in magnitude and we observe incoherent spatial patterns in the ′ ′ field observed throughout the turbulent reacting flow field of the combustion chamber, as evident in figure 4-(II) (George et al. 2018). Also, patches of slightly high values of ′ ′ occur in the wake of the bluff-body erratically as is evident from figure 4-(II)B.
In the vorticity field, we observe from experiments that small-sized vortices are shed in the recirculation zone (in the dump plane) and in the wake of the bluff-body. The vorticity (figure 4-(III)) is found to be high and positive in the outer shear layer owing to the anticlockwise vortices shed in this region. The vortices shed in the wake are formed due to the flow turning at the tip of the bluff-body from where we see the origin of the inner shear layer. Vorticity is high and negative in the inner shear layer indicating that clockwise vortices are shed in the wake. We note that, while the streaks of high vorticity in the inner and outer shear layer seem to be consistent in time, the shedding of the small vortices is rather erratic (George et al.

2018).
Results of the analysis of the inter-layer network constructed from the vorticity and the thermoacoustic power field during the occurrence of combustion noise are shown in figure 5. The node strength of a node in the ′ ′ -layer is determined from the inter-layer network by summing up the weights of all its inter-layer connections to nodes in the -layer, as discussed in Section 4. Thus, high node strength of a node in the ′ ′ -layer indicates that the thermoacoustic power generation at that location (node) is well correlated over short epochs with the hydrodynamic activity (vorticity) of the other locations to which it is connected via inter-layer links. Similarly, one can interpret the node strength of nodes in the vorticity layer as well.
We observe that the inter-layer node strengths in the ′ ′ -layer are high in the shear layer region and moderately high in the wake region, and low in the recirculation zone in the dump plane ( figure 5(a)) during the state of combustion noise. In contrast, in the -layer ( figure 5(b)), high inter-layer node strength is obtained in the recirculation zone in the dump plane and in the wake region, and very low node strength in the shear layers. Note that, we do not obtain any identifiable spatial patch of hub-nodes in figure 5(a,b). Such interlayer node strength distribution indicates that the thermoacoustic power generated in the shear layer region upstream of the bluff-body is correlated to the vortex shedding occurring in the recirculation zone in the dump plane and the wake of the bluff-body during the state of combustion noise. We note that the node strengths in the -layer are small, almost negligible, in the shear layer despite larger vorticity in these regions as observed in figure 4-(III). Clearly, thermoacoustic power generation in the combustor is less correlated to the vorticity fluctuations in the shear layer owing to the incoherence induced in the vorticity field by erratic vortex shedding during the state of combustion noise. Node strengths in the -layer are comparatively higher in the recirculation zone in the dump plane and the bluff-body wake, since the combustion of fuel-air mixture in the small vortices shed in these recirculation regions causes heat release.
It is interesting that, while thermoacoustic power generation is maximum in the wake ) with inter-layer degree ( ). exhibits a negative power law variation with where for lowdegree nodes the power law exponent is = −0.04, indicating weak disassortativity, while for high-degree nodes it is = −0.22 indicating disassortativity. (d) The inter-layer link-rank distribution of links between nodes ranked into bins in the ′ ′ -layer and the -layer. High density of inter-layer links exist between low ranking nodes in both layers. Also, significant links exist between low and high ranking nodes in both layers leading to disassortativity in the network.
downstream of the bluff-body (figure 4-(II)), inter-layer network analysis identifies a region of high node strength in the shear layer of the ′ ′ -layer, that is upstream of the bluff-body. Thermoacoustic power is generated in, both, the bluff-body wake and the recirculation zone in the dump plane due to the combustion of the unburnt fuel-air mixture entrained by vortices in these recirculation regions. However, the thermoacoustic power generated in the outer shear layer (as compared to that in the wake of the bluff-body) appears to be more strongly correlated with vortex shedding activity in the combustor.
Next, we study the topology of the inter-layer network obtained during the state of combustion noise using inter-layer network assortativity (as explained in Section 4) in figure 5(c). The plot of the inter-layer degree correlation function ( ) with the interlayer degree ( ) exhibits distinct power-law behaviors for low and high inter-layer degree range, i.e., the inter-layer network exhibits mixed assortativity. In the low-degree range, the degree correlation power-law exponent is negative but close to zero ( = −0.04) indicating very weak disassortativity. In the high-degree range, we observe strong disassortativity with Figure 6: Visualization of inter-layer links derived from the inter-layer network during the state of combustion noise for sample nodes. The green surface represents the -layer, while the red surface represents the ′ ′ -layer. The cut out region in each layer indicates the location of the bluff-body. Sample nodes are selected (a) in the dump plane in the -layer, (b) in the wake region in the -layer, (c) in the dump plane in the ′ ′ -layer, and (d) a node in the wake region in the ′ ′ -layer. The inter-layer connections of most nodes are spread over wide spatial locations indicating non-localized inter-subsystem interactions in the combustion chamber. = −0.22. Note that, such strong disassortativity may be a manifestation of the numerous links made by 'rich' nodes (nodes with very high node strengths) as discussed earlier (false disassortativity discussed in Section 4). Thus, we must examine the distribution of links between nodes with high or low node strengths in the two layers of the inter-layer network. Figure 5(d) shows the inter-layer link-rank distribution derived from inter-layer network during the state of combustion noise, as explained in Section 4. We obtain highest density of links between the lowest ranked bins (low-to-low ′ ′ − to − type), while significant inter-layer links also exist between nodes in the two layers with moderate and high ranking as well. That is, the density of inter-layer connections of all types (low-to-low, high-to-low, low-to-high and high-to-high) are significant and comparable. Also, the link-rank distribution shows that connections between high-degree nodes in each layer are not negligible, counterintuitive to the strong disassortative degree correlations obtained for the range of high-degree nodes ( figure 5(c)).
In summary, we infer that the inter-subsystem interactions are spread throughout the flow field and are not localised to any specific region in the combustion chamber. This interpretation is supported by the visualization of inter-layer links of sample points as shown Thermoacoustic power source appears erratically while thermoacoustic power sink exists predominantly in the shear layer region. Positive incoherent vorticity field exists in the dump plane and the outer shear layer, while a negative vorticity streak is observed in the inner shear layer. in figure 6. Figure 6(a) shows that the inter-layer links of a node selected in the dump plane of the -layer (green-colored layer) connect it to most nodes in the ′ ′ -layer (red-colored layer) throughout the flow field. Figure 6(b) shows that a node in the wake region in the -layer connects predominantly to nodes in the ′ ′ -layer in the bluff-body wake and in the shear layer regions upstream of the bluff-body (regions of high node strengths in the ′ ′ -layer, see figure 6(a)). Similarly, the inter-layer links of a node in the ′ ′ -layer selected in the recirculation zone (6(c)) are spread across in the -layer in the recirculation zone and the wake, but not in the shear layer regions as expected from the spatial distribution of node strengths in figure 5(b). Also, a node in the wake in the ′ ′ -layer (see figure 6(b)) connects predominantly to nodes in the wake in the -layer, with few links spread into the recirculation zone. Clearly, the hydrodynamic activity and thermoacoustic power generated at nodes in the wake region are densely inter-connected. Also, inter-subsystem interactions occur between spatial location upstream and downstream of the bluff-body. Note that, inter-layer links of most nodes in both layers are spread over wide spatial locations and are clearly non-localized.

Intermittency
During the state of intermittency, we observe bursts of periodic oscillations amidst aperiodic fluctuations in the acoustic pressure and the heat release rate signals (Nair et al. 2014). In figure 7 and figure 8, we show the variations in the ′ ′ (row-II) and vorticity fields (row-III) in the combustion chamber during epochs of aperiodic and periodic acoustic pressure dynamics, respectively.
During aperiodic dynamics in intermittency, we observe sinks of thermoacoustic power generation (negative values of ′ ′ are predominant) as evident from figure 7-(II). Note that, the ′ ′ field obtained for aperiodic dynamics during the state of intermittency is very different from that during combustion noise, since we do not obtain incoherent patterns of ′ ′ . Instead, we observe coherent sinks of ′ ′ with few erratic occurrences of ′ ′ source Thermoacoustic power generated during periodic epochs is an order of magnitude higher than during the aperiodic epochs. Source of ′ ′ exists predominantly in the shear layer and also in the wake close to the bluff-body. Large patches of high positive vorticity appear and disappear in the dump plane and the outer shear layer regions due to the shedding of large coherent structures in these regions.
predominantly in the shear layer region such as in figure 7-(II)A,D. On the other hand, during the periodic epochs of intermittency (see figure 8-(II)), we observe coherent regions that are source (or sinks) of thermoacoustic power generation with values which are one order of magnitude higher than that during the aperiodic epochs. These regions of coherent ′ ′ generation occur in the inner shear layer and in the wake close to the bluff-body as observed in figure 8-(II). Also, such coherent source and sinks of ′ ′ occur periodically. Further, from experiments we note that small vortices are shed during aperiodic dynamics in the state of intermittency. From figure 7-(III), we observe incoherence in the vorticity field fluctuations specifically in the recirculation zone. Unlike the vorticity field during the state of combustion noise, a consistent streak of vorticity does not exist in the outer shear layer. On the other hand, during periodic epochs of intermittency (see figure 8-(III)), we observe that the vortices formed in the recirculation zone are larger than those formed during the state of combustion noise or aperiodic epochs of intermittency. These large vortices are shed almost periodically for short epochs in the recirculation zone, after which they traverse to the wake via the shear layer and may break down into smaller vortices or impinge on the walls of the combustor and the bluff-body. Thus, we see distinct patches of large positive vorticity appear and disappear in the recirculation zone (figure 8-(III)B) and then in the outer shear layer (figure 8-(III)C). In contrast, the vortices shed at the tip of the bluff-body are small. So, we obtain an almost consistent streak of negative vorticity in the inner shear layer (figure 8-(III)A-D).
Next, we construct an inter-layer network considering, both, the epochs of periodic and aperiodic dynamics during the occurrence of intermittency. Figure 9 shows the results from the inter-layer network analysis of the spatio-temporal dynamics during the state of intermittency. The node strength distribution of both the layers ( figure 9(a,b)) is peculiar showing a distinct patch of very high node strengths (about 0.9) in the recirculation zone in (d) The inter-layer link-rank distribution of links between nodes ranked into bins in the ′ ′ -layer and the -layer. The colorbar indicates the percentage of links. The density of links connecting hubs in both layers (low-to-low type links) is very high. Also, low-to-high and high-to-low type of links are significant indicating that hubs in one region connect also to regions of moderate or low node strengths in the other layer.
the dump plane downstream of the backward-facing step (refer schematic figure 1(c)); while in other regions the node strengths are moderately high (about 0.6). These spatial patches of hub-nodes occur in similar regions in the dump plane in both the layers. We infer that localized regions of inter-subsystem interactions emerge during the state of intermittency. However, note that, we cannot infer the spatial extent of the inter-subsystem influence of these hub nodes using the spatial node strength distributions alone. Thus, we use the interlayer network assortativity and link-rank distributions to study the topology of the inter-layer network, and thus infer the pattern of inter-subsystem interactions in the combustor.
The inter-layer degree correlation function ( ) derived from the inter-layer network decreases with increase in the inter-layer degree ( ) and exhibits a power-law with the exponent = −0.39, implying strong disassortativity in the topology of the network (figure 9(c)). That is, hubs in one layer have many links to nodes with low node strengths in the other layer. In other words, assortativity analysis reveals that the nodes in the dump plane (hubs) in each layer have significant connections to the shear layer regions and the wake of the bluff-body (regions of low node strengths) in the other layer. However, strong network disassortativity does not rule out the possibility of hub-to-hub connections (as discussed in Section 4) and there are also connections between hubs in both the layers. As a result, we understand that the thermoacoustic power generated in the recirculation zone in the dump plane is strongly correlated with the vorticity dynamics in all regions of the combustor during the state of intermittency. Also, the shedding of vortices in the recirculation zone of the dump plane is correlated over short epochs with the thermoacoustic power generated in all regions of the combustor.
An important question remains unanswered; while the hubs in one layer connect to almost all regions in the other layer, are there significant links between the hub regions in both the layers? To answer this, we examine the inter-layer link-rank distribution. We note that low-tolow type of links have the highest density as evident from figure 9(d) (refer to the schematic in figure 3). Clearly, the density of inter-layer hub-to-hub connections is significant and a 'rich club' is formed in the inter-layer network, counter-intuitive to the strong disassortativity obtained in figure 9(c). We therefore infer that intense inter-subsystem interactions occur between the hydrodynamic, combustion and acoustic subsystems specifically in the localized hub-region in the dump plane during the state of intermittency.
Also, note that, there are significant density of links of the type high-to-low and lowto-high, where as, very few connections are of the high-to-high type (figure 9(d)). Such a link-rank distribution justifies the strong disassortativity (as in figure 9(c)) despite dense hub-to-hub inter-layer connections, since hubs in one layer connect to nodes of all ranks in the other layer. Also, since both high-to-low and low-to-high type of links are equally significant, we understand that the hubs in both the layers are equally influential. In other words, both, the thermoacoustic power generated and the vortices shed in the dump plane are equally crucial in fostering inter-subsystem interactions in the combustor.
We conjecture that, during the state of intermittency, the fuel-air mixture carried by the large vortices shed in the recirculation zone undergoes combustion in the recirculation zone itself. Remember, that large vortices are shed in the recirculation zone in the dump plane, while small vortices are shed in the wake of the bluff-body (figure 8-(III)). Higher vorticity and higher local strain rates induced in the flow by the large vortices promote better mixing of the fuel-air mixture, thus allowing combustion of fuel-air mixture in these vortices upstream of the bluff-body. As a result, intense and coherent heat release occurs in the dump plane and the subsequent downstream locations, that is, in the outer and inner shear layers (as evident in figure 8-(II)A,D). As a result, the vorticity dynamics of the recirculation zone is well correlated with the thermoacoustic power generation in the recirculation zone as well as in the regions downstream.
A visualization of the inter-layer links of sample nodes in the two layers is shown in figure  10. A node in the recirculation zone in the hub-region of the -layer (see figure 10(a)) has dense connections to nodes in the dump plane (hubs) of the ′ ′ -layer as well as to many nodes in the wake region. Similarly, a node in the hub-region of the ′ ′ -layer (see figure  10(c)) connects to nodes in all regions, dump plane shear layer and the wake. On the other hand, a node selected in the wake region in either layer has significant inter-layer connections only to the hub region in the other layer or with some nodes in its neighboring spatial locations in the wake (see figure 10(b,d)).
In summary, using multilayer networks, we identify localized hub-regions in the recirculation zone behind the backward-facing step (i.e., the dump plane) where intense intersubsystem interactions emerge during the state of intermittency. Note that, such significant inter-subsystem interactions were believed to occur only during thermoacoustic instability when order is established in the spatio-temporal dynamics of the combustor (Poinsot et al. 1987;Schadow & Gutmark 1992). However, we show using multilayer network analysis that Figure 10: Visualization of inter-layer links derived from the inter-layer network during the state of intermittency for sample nodes. The green surface represents the -layer, while the red surface represents the ′ ′ -layer. The cut out region in each layer indicates the location of the bluff-body. Sample nodes are selected (a) in the recirculation zone in the -layer, (b) in the wake region in the -layer, (c) in the recirculation zone in the ′ ′ -layer, and (d) a node in the wake region in the ′ ′ -layer. The hub-region in each layer has dense connections to the hub-region as well as regions of low and moderate node strengths in the other layer, indicating the widespread influence of hubs on inter-subsystem activity during the state of intermittency.
such interactions emerge much prior to the onset of order, during the state of intermittency. The interactions between the two subsystems are intense in the dump plane but not restricted to this region.
The influence of the vorticity dynamics and thermoacoustic power generation in the dump plane is widespread throughout the turbulent reacting flow field of the combustor. Thus, attacking the spatial collection of hubs of the inter-layer network identified in the dump plane is a suggested passive control strategy for mitigating thermoacoustic instability. Due to the large spatially extended inter-layer connections of the hub nodes in the dump plane, any perturbation introduced in the hub-region will disrupt the inter-subsystem interactions throughout the combustion chamber and thus mitigate the emergence of order in the spatiotemporal dynamics of the combustor.

Thermoacoustic instability
During the state of thermoacoustic instability, order is established in the spatio-temporal dynamics of the combustor. We obtain high-amplitude periodic oscillations in the acoustic pressure and heat release rate signals (figure 11-(I)). Coherent heat release occurs and ′ ′ is generated in large regions of the combustion chamber periodically (see figure 11-(II)A,C). Such sources of high thermoacoustic power occur predominantly in the recirculation zones in the dump plane and the wake of the bluff-body (such as in figure 11-(II)C). The value of thermoacoustic power generated in the combustion chamber during thermoacoustic instability (refer figure 11-(II)) is an order of magnitude higher than that obtained during periodic epochs of intermittency (refer figure 8-(II)), and about three orders of magnitude higher than that during the state of combustion noise (refer figure 4-(II)). Also, strong thermoacoustic power sinks are formed in the combustor during such dynamics (see figure  11-(II)D).
The vorticity field during thermoacoustic instability (figure 11-(III)) exhibits periodic spatio-temporal patterns, and the magnitude of the vorticity fluctuations are higher than those during the states of combustion noise or during aperiodic epochs of intermittency. We can observe a large region of high positive vorticity in the recirculation zone in the dump plane (see figure 11-(III)A). From experiments, we understand that very large vortices are formed upstream of the bluff-body; these vortices span the dump plane and the outer shear layer regions. These large coherent structures are shed periodically at a frequency close to the acoustic frequency (George et al. 2018) and subsequently impinge on the walls of the combustor or the bluff-body. Such large vortices then break down into smaller vortices downstream in the inner shear layer and wake, causing further mixing of unburnt reactants and hot products. Also, the impingement of vortices causes burning of the fuel-air mixture contained therein, leading to large coherent heat release in the combustor. Figure 12 shows the results from the analysis of the inter-layer network for the state of thermoacoustic instability. In the ′ ′ -layer ( figure 12 (a)), the node strengths are moderately high throughout the combustion chamber (in the range 0.65 to 0.7). The node strengths are highest in the dump plane followed by high node strengths in a large patch in the wake region The inter-layer link-rank distribution of links between nodes ranked into bins in the ′ ′ -layer and the -layer. The colorbar indicates the percentage of links. Almost all nodes in the dump plane and the wake in the ′ ′ -layer are low-ranking nodes (high-node strengths). These nodes predominantly connect to low-ranking nodes in the -layer that occur in the shear layer regions and a small patch in the dump plane. of the ′ ′ -layer. In the -layer, the node strengths are very high (around 0.8) in the outer and inner shear layer as well as in a patch in the recirculation zone in the dump plane. Nodes in the wake region have moderately low node strengths (around 0.65).
Further, we find that the inter-layer degree correlation function ( ) does not exhibit a monotonic increase or decrease with the inter-layer degree ( ) and the degree correlation coefficient is found to be close to zero ( = −0.02) ( figure 12(c)). The inter-layer network topology is neutral, that is neither assortative nor disassortative, owing to the fact that almost all nodes in the ′ ′ -layer and several nodes in the -layer have similar values of node strengths.
Moreover, using the inter-layer link-rank distribution (figure 12(d)), we find that inter-layer links are predominantly low-to-low and low-to-high ( ′ ′ − to − ) type (refer schematic in figure 3). Note that, almost all nodes in the ′ ′ -layer have very high node strengths and are thus binned together as low ranking nodes. Thus, from the spatial distribution of node strengths and the link-rank distribution, we infer that the inter-layer links are densely Figure 13: Visualization of inter-layer links derived from inter-layer network during the state of thermoacoustic instability for sample nodes. The green surface represents the -layer, while the red surface represents the ′ ′ -layer. The cut out region in each layer indicates the location of the bluff-body. Sample nodes are selected (a) in the hub region in recirculation zone of the -layer, (b) in the farther downstream end of the inner shear layer of the -layer, (c) in the recirculation zone of the ′ ′ -layer, and (d) a node in the wake region of the ′ ′ -layer. Dense inter-layer links exist between the nodes in the shear layer regions of the -layer and the nodes in the dump plane and wake regions of the ′ ′ -layer. concentrated in the dump plane and the shear layer regions of the -layer but spread throughout the flow field in the ′ ′ -layer.
Inter-layer network analysis shows that, the vorticity dynamics in the dump plane and the shear layer regions is strongly correlated with the thermoacoustic power generated in the recirculation zone in the dump plane and the wake of the bluff-body. We understand that, the large vortices shed in the dump plane enhance the strain rates in the fluid and thus enhance fine-scale mixing in this region (Poinsot et al. 1987;Schadow & Gutmark 1992). As a result, most of the unburnt reactants and the hot products mix well and the unburnt fuel-air mixture undergoes combustion in the dump plane. Also, note that during thermoacoustic instability, vortices are being shed and convected at a high frequency (close to the acoustic frequency). Thus, these large vortices also convect some of the unburnt fuel-air mixture downstream before the mixture undergoes combustion in the dump plane. Then, these large vortices eventually break into smaller vortices that enhance mixing in the flow, leading to the combustion of the unburnt reactants in the inner shear layer and the wake. The thermoacoustic power generated in the dump plane and the wake of the bluff-body are therefore strongly correlated to the vorticity dynamics in regions where large coherent structures are shed (in the dump plane) and regions where such vortices impinge or break down into smaller vortices (in the shear layers), respectively. Also, moderately strong correlations exist between the thermoacoustic power field and the vorticity dynamics in the wake owing to the smaller vortices that dominate the hydrodynamics of the wake. Figure 13 shows the inter-layer links of sample nodes upstream and downstream of the bluff-body in both the layers. A node in the hub region in the dump plane of the -layer ( figure 13(a)) has inter-layer connections predominantly to the hub region in the dump plane of the ′ ′ -layer and also some inter-layer links to the wake region. Figure 13(b) shows that a hub-node from the -layer in the inner shear layer has links predominantly to the wake region and a few links to the dump plane in the ′ ′ -layer. Note, we conjectured that the breaking down of large coherent structures into smaller vortices in the shear layer and the impingement to the walls lead to coherent thermoacoustic power generation in the wake. Thus, we see that the vorticity dynamics at this node in the inner shear layer of the -layer ( figure 13(b)) connects predominantly to the wake region in the ′ ′ -layer. Further, a node in the dump plane in the ′ ′ -layer (figure 13(c)) has predominant inter-layer connections to nodes in the dump plane and the shear layer regions of the -layer (high node strength regions, see figure 12(b)). Also, figure 13(d) shows that a node in the wake region in the ′ ′ -layer has inter-layer connections with nodes in the shear layer regions as well as the recirculation zone in the -layer. Thus, visualization of the inter-layer links supports our inferences about the topology of the inter-layer network.
The results obtained from multilayer network analysis during the various dynamical states in a turbulent thermoacoustic system are summarized in table 1. We note that, the large coherent vortices shed during the occurrence of thermoacoustic instability play a pivotal role in the emergence of the large coherent source of thermoacoustic power generation in the combustion chamber. Our findings thus support the results of Poinsot et al. (1987) that the reaction zone and coherent heat release patches trail behind the large coherent structures during the state of thermoacoustic instability. Thus, using multilayer network analysis, we conclude that the role of coherent structures is critical in fostering inter-subsystem interactions between the hydrodynamic, acoustic and combustion subsystems during thermoacoustic instability.

Discussion
Results derived from multilayer network analysis in the current study complement and potentially explain the observations from some previous studies. For example, we are able to identify a critical region of interaction between thermoacoustic power and vorticity field in the dump plane during intermittency, which is very similar to the critical region identified by Premchand et al. (2019) corresponding to the dominant acoustic mode using flow decomposition techniques. Also, a similar critical region was identified by Krishnan et al. (2021) as the hubs of vorticity networks during the state of thermoacoustic instability. Further, Roy et al. (2021) found a similar region in the dump plane where intense turbulent velocity fluctuations occur during the occurrence of thermoacoustic instability.
Further, Krishnan et al. (2021) and Roy et al. (2021) found significant suppression of thermoacoustic instability even for low flow rates of targeted microjet injections in the recirculation zone behind the bluff-body. However, perturbing any other region such as patches identified at the top of the bluff-body or in the wake, does not readily aid the mitigation of thermoacoustic instability. While the authors (Krishnan et al. 2021;Roy et al. 2021) intended to target regions of intense hydrodynamic activity, they inadvertently perturbed the region of significant inter-subsystem interactions, which we identify, occurs in the recirculation zone. Thus, we conjecture that targeted attack in regions where dense hub-to-hub inter-layer connections exist can lead to the mitigation of thermoacoustic instability. Such perturbations would disrupt the strongest feedback interactions between the multiple subsystems such as the hydrodynamic field and the combustion-driven acoustics. To summarize, the use of multilayer networks potentiates physical insight into the occurrence as well as the mitigation of thermoacoustic instability.

Conclusions
In this work, we investigate the interactions between coherent structures in the hydrodynamic field and the generation of thermoacoustic power during the emergence of order amidst chaos in the spatio-temporal dynamics of a turbulent thermoacoustic system. We construct a multilayer network comprising two layers, the -layer and the ′ ′ -layer, where the inter-layer connections are derived from the short-window correlation between the vorticity fluctuations ( ) and the thermoacoustic power generated ( ′ ′ ) in a turbulent bluff-body stabilized combustor. The spatial distribution of node strengths derived from inter-layer networks helps us identify regions of intense inter-subsystem interactions (hubs of the inter-layer network) during various dynamical states. Further, we interpret the topology of the interlayer connections between the two layers by examining the inter-layer degree correlations. Also, by studying the distribution of links between nodes ranked according to their interlayer degree (link-rank distribution), we identify if significant hub-to-hub inter-layer links exist. The topology of the inter-layer network helps decipher the structure of inter-subsystem interactions in the system.
The interactions between the underlying hydrodynamic flow and the thermoacoustic power generated in the combustion chamber are remarkably different during distinct dynamical states. Such interactions are spread across the combustion chamber and are non-localised during combustion noise. During intermittency, a collection of nodes with very high node strengths (hubs) emerges in the dump plane in both the layers of the inter-layer network. Dense hub-to-hub connections exist between the two layers indicating intense interactions between large coherent structures and the acoustically coupled heat release rate fluctuations in the dump plane of the combustor. Further, hubs in one layer also have significant inter-layer connections to lower node strength regions of the combustor. We infer that, the influence of the coherent structures shed and thermoacoustic power generated in the dump plane strongly influences the dynamics downstream in the combustion chamber. Finally, during thermoacoustic instability, inter-layer network analysis shows that vorticity in the dump plane and the shear layer region is strongly correlated with the thermoacoustic power generated in the dump plane and the wake of the bluff-body.
Observing the simultaneous emergence of coherent structures and intense heat release in the combustor during the state of thermoacoustic instability, Poinsot et al. (1987) and Schadow & Gutmark (1992) conjectured that strong interactions occur between the hydrodynamic, acoustic and combustion subsystems during this state. Our analysis using multilayer network shows that such inter-subsystem interactions indeed exist and are intense when order is established in the dynamics of the combustor, i.e., during the state of thermoacoustic instability. Moreover, we show that such intense inter-subsystem interactions emerge much prior to the onset of order in the spatio-temporal dynamics of the combustor, i.e., during the state of intermittency. In summary, the approach of multilayer networks helped reveal the rich pattern of inter-subsystem interactions in turbulent thermoacoustic systems during various dynamical states.