# Impact of precessing vortex core dynamics on the thermoacoustic
instabilities in a swirl stabilized combustor

Ashwini Karmarkar  Department of Mechanical Engineering, Pennsylvania State
University, University Park, PA  Saarthak Gupta  Department of Aerospace
Engineering, Indian Institute of Science, Bengaluru, India  Isaac Boxx  DLR -
German Aerospace Center, Stuttgart, Germany  Santosh Hemchandra  Department of
Mechanical Engineering, Pennsylvania State University, University Park, PA
Jacqueline O’Connor  Department of Aerospace Engineering, Indian Institute of
Science, Bengaluru, India

###### Abstract

Global instabilities in swirling flows can significantly alter the flame and
flow dynamics of swirl-stabilized flames, such as those in modern power
generation gas turbine engines. In this study, we characterize the interaction
between the precessing vortex core (PVC), which is the consequence of a global
hydrodynamic instability, and thermoacoustic instabilities, which are the
result of a resonant coupling between combustor acoustics and the unsteady
heat release rate of combustion. This characterization is performed using
experimental data obtained from a model gas turbine combustor system employing
two concentric swirling nozzles of air, separated by a ring of fuel injectors,
operating at 5 bar pressure. The flow split between the two streams is
systematically varied to observe the impact of flow structure variation on the
flow and flame dynamics. High-speed stereoscopic particle image velocimetry,
OH planar laser-induced fluorescence, and acetone planar laser-induced
fluorescence are used to obtain information about the velocity fields, flame,
and fuel flow behavior, respectively. Spectral proper orthogonal decomposition
and spatial frequency analysis are used to identify and characterize the
dominant oscillation mechanisms driving the system. Three dominant modes are
seen: two thermoacoustic modes and the precessing vortex core. Our results
show that in the cases where the frequency of the PVC overlaps with either of
the thermoacoustic modes, the thermoacoustic modes are suppressed. A weakly
nonlinear asymptotic analysis shows that the suppression of the axisymmetric
shear layer shedding, and hence thermoacoustic mode, is the result of a
nonlinear coupling between the PVC and the axisymmetric mode of the swirling
jet. Evolution equations for both the symmetric and PVC modes are derived to
show the controlling parameters that drive this suppression. We conclude by
discussing ways in which thermoacoustic instability suppression can be
achieved through combustor flow field design.

##  1 Introduction

For modern power generation gas turbines, increasingly stringent requirements
on $NO_{x}$ emissions have driven a technology shift towards lean premixed
burners. While fuel lean, premixed systems allow combustors to achieve
significantly lower $NO_{x}$ emissions, these flames are inherently less
stable and hence more susceptible to acoustic disturbances, which can be a
precursor to the occurrence of combustion instabilities (Correa, 1998;
Lefebvre, 1998). Combustion instability is a thermoacoustic instability that
arises from a coupling between the resonant acoustic modes of the combustor
and the unsteady heat release rate of the flame. Occurrence of combustion
instability is a limiting factor to engine performance and operability, and,
in extreme cases, can lead to severe damage of engine hardware (Lieuwen &
Yang, 2005; Candel, 2002; Correa, 1998).

The acoustic modes of a combustor can couple with the unsteady rate of heat
release through fluctuations in pressure, velocity, and/or equivalence ratio
(Lieuwen, 2012). Velocity coupling mechanisms can be driven by large-scale
coherent structures in the flow that distort the flame and consequently alter
the heat release rate (Paschereit et al., 1999; Gonzalez et al., 2005; Poinsot
et al., 1987; Shanbhogue et al., 2009; Kang et al., 2007). For instance, it
has been shown that vortices in the flow can entrain fresh reactants and cause
the flame to roll-up, which sharply increases the flame surface area, thereby
creating a heat release pulse. This coherent oscillation in heat release feeds
back into the acoustic pressure fluctuations and creates a self-sustaining
thermoacoustic instability (Renard et al., 1999; Lee et al., 1993; Ghoniem &
Givi, 1988; Poinsot et al., 1987). A detailed review of the different
mechanisms by which flames interact with vortical structures is provided by
Renard et al. (2000).

Equivalence ratio fluctuations have also been shown to be drivers of
combustion instabilities in lean premixed combustors (Lieuwen & Zinn, 1998;
Kim et al., 2010; Shreekrishna et al., 2010; Auer et al., 2006; Bluemner et
al., 2019). Pressure oscillations in the combustor interact with the fuel
supply and cause coherent oscillations in the equivalence ratio of the
reactant mixture, which when convected to the flame, can cause pulsations in
the heat release rate. Lee et al. (2000) have shown, using an infrared
absorption technique to measure equivalence ratio fluctuations, that the
equivalence ratio fluctuations are strongly linked to the heat release
fluctuations in an unstable combustor and can thereby play a significant role
in driving the instability.

In partially-premixed, swirl-stabilized systems used in gas turbine
combustors, such as the combustor considered in this work (Stöhr et al.,
2019), complex hydrodynamic and thermo-chemical processes are involved, which
may lead to multiple coupling pathways being present in the system. The
dynamics of swirling flows are especially complex because swirling flows can
exhibit several modes of hydrodynamic instability. Swirling flows are
typically characterized by a non-dimensional parameter known as the swirl
number, which is the ratio of the axial flux of tangential momentum to the
axial flux of axial momentum (Gupta et al., 1984). Strongly swirling flows,
such as those seen in gas turbine combustors, are characterized by the
formation of a central recirculation zone, or a ‘vortex breakdown bubble’
(Harvey, 1962; Hall, 1972). This recirculation zone constantly supplies hot
products to the base of the flame, thereby enhancing its static stability.

Increasing the degree of swirl can cause this recirculation zone to precess
about the axis of symmetry, forming what is known as the precessing vortex
core (PVC) (Syred, 2006). PVCs are the the manifestation of a global
instability obtained through a super-critical Hopf bifurcation induced by
vortex breakdown (Manoharan et al., 2020; Oberleithner et al., 2011).
Manoharan et al. (2020) have shown, using experimental data and results from a
weakly nonlinear stability analysis, that beyond a critical swirl number, the
frequency and square of the amplitude of the PVC oscillations scale linearly
with swirl number. They have also confirmed from ab initio estimates of the
coefficients of the associated Stuart-Landau equation for the oscillation
amplitude of the PVC that it is a stable limit cycle flow oscillation that
results from a supercritical Hopf bifurcation in the flow state.

In combustor systems, the occurrence of a PVC is a function of swirl number,
flame shape, fuel/air mixing, combustor configuration, and equivalence ratio
(Syred, 2006). PVC dynamics can significantly influence the flame and flow
dynamics in a number of ways (Candel et al., 2014). The shear layer
oscillations due to the perturbations caused by the precession of the vortex
breakdown bubble have been shown to interact with the flame and cause flame
surface distortion (Stöhr et al., 2012, 2015). The occurrence of PVCs can
also, in some cases, impact the stability and structure of swirl stabilized
flames (Oberleithner et al., 2015). It has also been shown that PVCs can
enhance fuel/air mixing (Freitag et al., 2005; Stöhr et al., 2015).

The primary focus of this study is to characterize the impact of the PVC on
the thermoacoustic modes of the combustor in a partially-premixed flame at
elevated pressure. The interaction of hydrodynamic instabilities, such as the
PVC and symmetric acoustic modes, is of significant interest as both these
modes can be simultaneously observed in many combustor configurations.
Hemchandra et al. (2018) proposed two possible coupling pathways between self-
excited hydrodynamic instability modes and the combustor acoustic modes that
can lead to combustion instability. These results indicate that the strength
of the coupling between the hydrodynamic and acoustic modes depends on how
close the instability frequency is to the combustor acoustic eigenfrequency,
in addition to how effectively acoustic modes can excite the flow in the
regions of hydrodynamic mode receptivity.

Steinberg et al. (2010) showed that in a case where both the PVC and
thermoacoustic mode are simultaneously present, the PVC can undergo axial
extension and contraction at the thermoacoustic frequency. Moeck et al. (2012)
and Steinberg et al. (2010) reported the presence of a spectral peak at an
interaction frequency that corresponds to the difference between the PVC and
thermoacoustic frequencies and suggested that this interaction frequency is
likely a consequence of nonlinear phenomena. Lückoff & Oberleithner (2019)
studied the impact of actuation of helical flow modes using a centerbody on
their interaction with thermoacoustic modes. Their results showed that, in the
case of a partially-premixed flame, this actuation can significantly reduce
thermoacoustic instability. They suggested that this suppression of
thermoacoustic instability modes is likely a consequence of suppression of
equivalence ratio fluctuation due to enhanced mixing caused by the helical
flow mode excited. By contrast, in some cases, the self-excited PVC mode is
suppressed in reacting flows (Oberleithner et al., 2015, 2013). Taamallah et
al. (2016) showed that the presence of a PVC in a premixed swirl-stabilized
combustor can provide a low velocity path within the combustor flow field
resulting in a flame shape change. This new flame shape, in turn, causes the
combustor acoustic mode to couple with the flame, resulting in thermoacoustic
oscillations.

Previous work from some of the present authors showed, through a spectral
analysis of a swirling flow, that response of the shear layers to acoustic
forcing in an isothermal swirling jet can be suppressed in the presence of a
PVC (Mathews et al., 2016). This result was explained by means of a linear
stability analysis, showing that the formation of the PVC results in increased
shear layer thickness. This increased thickness leads to a progressive
weakening of the receptivity of the axisymmetric Kelvin-Helmholtz mechanism to
imposed forcing (Frederick et al., 2018).

Put together, the existing literature shows that the dynamics of the PVC can
have a critical impact on the performance of swirl-stabilized combustion
systems. Furthermore, the PVC can interact with the symmetric acoustic modes
present in the combustor in a number of ways, depending on the operating
conditions and combustor configuration. This study extends these concepts to a
more realistic configuration with a dual-annular swirl injector, commonly
found in gas turbine combustors (Joshi et al., 1998), partially-premixed fuel
injection, and elevated pressures. Using an optically accessible model swirl-
stabilized gas turbine combustor configuration, operating at an elevated
pressure of 5 bar, we study the interaction between the dominant oscillation
modes present in the system and observe the evolution of this interaction at a
range of flow conditions. Our results show that, in the configuration studied,
two dominant thermoacoustic modes are present. In cases where all three modes
are simultaneously present, the PVC oscillation mode is generally the
governing mode in the flow field. Furthermore, whenever the frequency of the
PVC approaches that of either thermoacoustic mode, the thermoacoustic mode is
subsequently suppressed.

In order to understand the interaction between the PVC and the thermoacoustic
modes, we extend the weakly nonlinear analysis of (Manoharan et al., 2020) to
variable density flows and include a stable hydrodynamic axisymmetric mode
weakly forced by thermoacoustic oscillations near its natural frequency.
Similar analysis has been performed recently by Rigas et al. (2017) in the
context of axisymmteric wakes. We obtain a closed-form expression for the
steady-state amplitude of the axisymmetric mode oscillating at the
thermoacoustic oscillation frequency. This steady-state amplitude has a strong
dependence on the difference between the thermoacoustic and the natural PVC
frequencies. The outcome of this dependence is that as the PVC frequency
approaches the thermoacoustic oscillation frequency, it causes the suppression
of the axisymmetric hydrodynamic, which can then result in the suppression of
the thermoacoustic oscillation. The analysis also shows that the nonlinear
coupling between the PVC and the axisymmetric hydrodynamic mode results in
helical oscillations at the interaction frequencies, i.e., the sum and
difference of the thermoacoustic and PVC oscillation frequencies. We restrict
our theoretical analysis in the present paper to a qualitative description of
the interaction between the PVC and the axisymmetric modes. The time-averaged
velocity and density fields are required to quantify the various contributions
to the steady-state solution. While it would be interesting to analyze various
contributions and the influence of nonlinear coupling, the spatial
distribution of the time-averaged density field, which has a significant
impact on the hydrodynamic modes (Manoharan et al., 2015; Emerson et al.,
2012; Manoharan & Hemchandra, 2014), is difficult to obtain from experiments
of the present kind with sufficient quantitative accuracy. Therefore,
quantitative confirmation of the theoretical analysis is beyond the scope of
this paper.

##  2 Experimental Methods

###  2.1 Burner and operating conditions

The burner used in the current study is shown in Fig. 1; its geometry has been
described in previous publications (Geigle et al., 2015, 2017). For
completeness, it is briefly described again here. The injector consists of a
pair of annular air swirl nozzles separated by a ring of 60 fuel channels,
each with an area of 0.5$\times$0.4 $mm^{2}$, for the introduction of gaseous
fuel, ethylene. The central air nozzle has a diameter of 12.3 mm. The annular
outer air nozzle measures 19.8 mm in diameter. The combustion chamber has a
square cross section measuring 68$\times$68 $mm^{2}$. It is 120 $mm$ long and
enclosed by 3-mm-thick quartz windows for optical access. The water-cooled
dome of the combustion chamber has a cylindrical exhaust hole - diameter 40
mm, length 24 $mm$ \- linked to the combustion chamber by a curved tube. The
combustor was mounted in an optically accessible pressure vessel.

Figure 1: Experimental Setup

The flows of air (central, annular and seeding) and fuel were supplied by
separate mass ﬂow controllers (Bronkhorst). The flow rate of air is 841.7 SLM
and flow rate of fuel is 37.9 SLM, resulting in a thermal power of 38.4 kW.
The flow controllers were calibrated in-house, resulting in an accuracy of
better than 1.0% of the controllers’ maximum ﬂow. The majority of the ethylene
flow was bubbled through a temperature-controlled reservoir (T=$25^{\circ}$C,
298 K) filled with acetone, which provided an acetone-saturated ethylene flow
of the desired composition. For the present study, the combustor was operated
with ethylene/acetone fuel at a pressure of 5 bars and $\phi=0.67$. The air
split, defined as the fraction of total air flow passing through the center
nozzle, was varied from 0.2 – 0.5. The Reynolds number based on the diameter
of the outer nozzle is $19,490$.

###  2.2 Measurement techniques

Time-resolved stereo particle image velocimetry (sPIV), planar laser-induced
fluorescence of OH (OH-PLIF), and acetone-PLIF were applied simultaneously
with a repetition rate of 10 kHz to capture flame/vortex/fuel interaction. Run
duration was 1 s in the technically premixed ethylene/air swirl flame
stabilized in the model of the aero-engine combustor. This diagnostic system
has been described previously in the scientific literature (Litvinov et al.,
2021). For completeness, a brief description is included below.

####  2.2.1 Particle Image Velocimetry

Velocity fields were measured using a two-camera, stereoscopic particle image
velocimetry system. The sPIV system used a dual-cavity, diode-pumped, solid
state laser (Edgewave, IS200-2-LD, up to 9 mJ/pulse at 532 nm) and a pair of
highspeed CMOS cameras (Phantom V1212). The cameras were mounted on opposite
sides of the laser sheet, looking down into the combustor. PIV image pairs
were acquired at 10 kHz and 640$\times$ 800 pixel resolution. The PIV
measurement domain spans across most of the width of the combustion chamber
($-27mm<x<27mm$). In the axial direction, the domain spans $5mm<y<27mm$, as
shown in Fig. 1. Time separation of the laser pulses was $\delta t=10\mu s$.
The beam was formed into a sheet using a pair of cylindrical lenses ($f=-38$
mm and 250 mm) and thinned to a waist using a third cylindrical lens ($f=700$
mm). Both combustion air flows were seeded with titanium dioxide ($TiO_{2}$)
particles of nominal diameter 0.5 $\mu m$. Image mapping, calibration, and
particle cross-correlations were completed using a commercial, multi-pass
adaptive window offset cross-correlation algorithm (LaVision DaVis 8.4). Final
interrogation window size and overlap were 24x24 pixels and 50%, respectively,
for a spatial resolution of 1.9 mm and vector spacing of 0.95 mm. The absolute
values of uncertainty of instantaneous velocities based on the correlation
statistics in DaVis were estimated to be 0.1 $m/s$—1.1 $m/s$ for the in-plane
components and 0.3 $m/s$—2.3 $m/s$ for the out-of-plane component. On average,
uncertainty values for the in-plane components are in the range of 0.04
$m/s$—0.45 $m/s$ and the uncertainty values for the out-of-plane component are
in the range of 0.10 $m/s$—0.9 $m/s$.

####  2.2.2 OH- and acetone-PLIF

The OH-/acetone-PLIF imaging system is based on a frequency-doubled dye laser,
pumped by high-speed, pulsed Nd:YAG laser (Edgewave IS400-2-L, 135 W at 532 nm
and 10kHz) and a pair of intensified high-speed CMOS camera systems. The dye
laser system (Sirah Credo) produced 5.3 – 5.5 W at 283 nm and 10 kHz
repetition rate (0.53-0.55 mJ/pulse). The dye laser was tuned to excite the
Q1(9) and Q2(8) lines of the
$A^{2}\Sigma^{+}-X^{2}\Pi(v^{\prime}=1,v^{\prime\prime}=0)$ band. These
transitions merge at high pressure due to increased collisional line
broadening, which mitigates to some degree fluorescence signal loss due to
collisional line broadening. The laser wavelength was monitored continuously
throughout the experiments using a photomultiplier tube (PMT) mounted to a 10
cm monochrometer and a premixed, laminar reference flame. The 283 nm PLIF
excitation beam is formed into a sheet approximately 40 mm (high) x 0.2 mm
(thick) using three fused-silica, cylindrical lenses (all anti-reflective
coated to maximize transmission). The laser sheets of the OH-/acetone-PLIF and
PIV systems were overlapped by passing the (green, 532 nm) PIV sheet through
the final OH-PLIF turning mirror.

The same laser was used to excite fluorescence signal in both the OH- and the
acetone-PLIF imaging systems. For acetone-PLIF, fluorescence signal was imaged
via a CMOS camera (LaVision HSS8), an external two-stage intensifier (LaVision
HS-IRO), equipped with 85 mm focal length, f/1.2 (Canon) objective and a band-
pass interference filter. The filter (LOT, 450FS40-50) was centered at 450 nm
for detection of the acetone fluorescence and had a bandpass of $\pm 20$ nm.
OH-PLIF fluorescence signal was imaged using a similar highspeed CMOS camera
(LaVision HSS8) and external two-stage intensifier (LaVision HS-IRO) from the
opposite side of the combustor. The OH-PLIF camera was equipped with 64 mm
focal length, f/2 (Halle) UV-objective and a high transmission bandpass
interference filter. The OH-/acetone-PLIF measurement domain was slightly
larger than that of the PIV system, spanning ($-27mm<x<27mm$, $0mm<y<30mm$),
as shown in Fig. 1.

Previous studies have demonstrated the feasibility of using acetone as fuel
tracer in a dual-swirl GT model combustor similar to that used in the present
study (Stöhr et al., 2015, 2017). These studies indicated that replacement of
10% vol of the original (methane) fuel by acetone vapor did not significantly
change the flame shape or dynamics. Although the present study uses ethylene
rather than methane, it is expected to have a similarly negligible effect on
the flame shape and dynamics. Data from the acetone PLIF signal is only
considered until downstream distances of $z=15$ mm, as the fuel is largely
consumed by that point. Signal downstream of that point is likely the result
of fluorescence of polycyclic aromatic hydrocarbons or even weak soot
luminescence, and so is not analyzed in the frequency analysis discussed in
Sec. 5.

###  2.3 Data analysis

In this study, we predominantly use spectral proper orthogonal decomposition
(SPOD) to characterize the frequency domain behavior of the flame and flow
field. Proper orthogonal decomposition (POD) (Berkooz et al., 1993) is an
eigenvalue decomposition that yields energy-ordered modes that can be used to
extract coherent structures from a flow field. Spectral POD is derived from a
space-time POD problem for statistically stationary flows and leads to modes
that each oscillate at a single frequency. The SPOD method used here is taken
from Towne et al. (2018). We use SPOD to extract high-energy coherent
structures rather than using a traditional snapshot POD method, which provides
spatially-orthogonal, energy-ordered modes of the data, because POD does not
discern motions based on their spatiotemporal evolution, just their spatial
correlations. As a result, tonal but spatially uncorrelated or low-energy
oscillations are not necessarily separated from other motions in the POD, but
can be identified in the SPOD. SPOD achieves this by calculating the
eigenvalue problem on the cross-spectral density tensor, which results in
spatial modes analogous to this from snapshot POD, but also provides
frequency-resolved energy information about each of these modes. Towne et al.
(2018) rigorously show that SPOD is an optimized dynamic mode decomposition
(DMD) (Schmid, 2010) method for stationary flows.

##  3 Experimental Results

###  3.1 Time-averaged flow and flame profiles

The time-averaged streamwise velocity contours with streamlines, obtained from
the sPIV measurement, for all air split conditions are shown in Fig. 2. In all
cases, a region of negative axial velocity, the central recirculation zone,
can be seen close to the axis of symmetry. As air split is increased, this
recirculation becomes stronger and more well defined.

Figure 2: Time-averaged stream-wise velocity component with streamlines.

Fig. 3 shows the variation of the time-averaged axial velocity along the
centerline and swirl number, derived from experimental data, $S_{d}$, as
defined in eq. 1, with downstream distance at every air split condition. To
calculate $S_{d}$, the axial flux of azimuthal velocity is radially integrated
from the $u_{z}=0$ contour at the outer boundary of the recirculation zone, at
$r=r_{c}$, to the edge of the frame of data at $r=R$, at each downstream
distance.

$S_{d}=\frac{\int_{r_{c}}^{R}u_{\theta}(z,r)u_{z}(z,r)rdr}{\int_{r_{c}}^{R}u_{z}(z,r)^{2}rdr}$
(1)

The plots show that as air split increases, the strength of recirculation
increases since the minimum velocity inside the vortex breakdown bubble
systematically decreases. Furthermore, the streamwise position of the
recirculation zone stays largely constant with varying air split. However, the
swirl number is relatively constant for all the cases, which indicates that
the change in the recirculation zone strength is a result of the change in the
critical swirl number for vortex breakdown with air split, rather than a
change in swirl number itself; this result is discussed further in Section 4,
as is the impact of the strength and location of the recirculation zone on the
hydrodynamic stability of the flow.

Figure 3: Centerline time-averaged axial velocity (left) and swirl number
(right) as a function of downstream distance for every air split condition.

The time-averaged flame structure (obtained from time-averaging the OH-PLIF
images) is shown in Fig. 4 for an air split of 0.4. Since this is a partially-
premixed flame, and the operating pressure is relatively high, the OH-PLIF has
a low signal-to-noise ratio and so the instantaneous flame shape is not always
apparent The M-flame shape is evident from the time-averaged image in Fig. 4,
but the figure should not be interpreted quantitatively. The time-averaged
fuel distribution (obtained from acetone-PLIF measurements) is shown also in
Fig. 4 at the same air split of 0.35. The time-averaged profiles for the flame
and fuel distribution do not change much with air split and so only one
condition is provided here. In this case, as in all cases, it’s evident that
the flame has an M-flame shape that is fueled by the acetone jet located
between the two swirl passages.

Figure 4: Time-averaged flame structure (left) and time-averaged fuel
distribution profile (right) for air split=0.4 case

###  3.2 Spectral characterization of flow fields

The time-averaged flow profiles suggest that increasing air split leads to
strengthening of recirculation in the vortex breakdown region. This increase
in recirculation strength can lead to changes in the hydrodynamic instability
characteristics of the flow. To gain more insight into these dynamics of the
flow, and to identify the governing oscillation modes in the system, the
frequency domain characteristics are analyzed.

The combustor acoustic pressure spectra, obtained by performing an ensemble-
averaged Fourier transform and calculating the power spectral density of the
pressure data at each air split condition, is shown in Fig. 5. Pressure is
measured near the outlet of the combustor. The power spectral density shows
multiple distinctive peaks at each condition, indicating the presence of
multiple sources of oscillation, which is characteristic of reacting systems
exhibiting hydrodynamic and thermoacoustic instability. In the lower air-split
cases, a strong peak can be seen at around 600 Hz, which is identified to be
the thermoacoustic instability, marked by a circle. As air split increases,
the peak at around 600 Hz persists, but is not as strong, and new peaks can be
seen at around 800 Hz, which is later identified as as second thermoacoustic
mode, marked with a triangle. Furthermore, harmonics of the dominant peaks can
also be seen. In order to identify coherent modes in the flow, the spatial
profiles of the flow fields at each of these modes are calculated using SPOD.

Figure 5: Acoustic pressure spectra at each air split condition.

In this study, we use spectral proper orthogonal decomposition (Towne et al.,
2018) to characterize the frequency domain behavior of the flow and to extract
the high energy coherent motion in the flow field. The eigenvalue spectra of
the velocity data at every air split condition is shown in Fig. 6. The red
curve in each plot depicts the highest energy mode and the gray curves depict
subsequent modes in order of decreasing energy. In all cases, most of the
coherent content is captured by the first mode, which is why the mode shapes
of the first mode are analyzed. Similar to the acoustic pressure spectra,
multiple peaks can be seen at each condition, which indicates that there are
multiple coherent oscillations in the system. In order to understand the
nature of the oscillation corresponding to every peak, the spatial mode shapes
of the first mode at the respective peak frequencies are reconstructed at
different operating conditions.

Figure 6: Eigenvalue spectra of velocity field obtained from SPOD.

Fig. 7 shows the spatial mode shapes of the radial velocity at the dominant
peak frequencies at the lowest air split (air split=0.2) case. The spatial
structure of the first mode at $f=664.1$ Hz shows concentrated regions of
fluctuation on either side of the central axis. Since this is a mode shape of
radial velocity, this spatial structure corresponds to a symmetric
oscillation, which is characteristic of the flow response to a longitudinal
thermoacoustic oscillation; we will refer to this mode as an $m=0$ mode, where
$m$ is the azimuthal wavenumber of the flow disturbance. In this case, the
longitudinal acoustic oscillation provides a symmetric disturbance to the flow
and so the flow responds symmetrically, as shown in previous swirling flow
studies by O’Connor & Lieuwen (n.d.). The mode at the second frequency,
$f=1309$ Hz, is the first harmonic of that mode, which can be seen both by the
doubling in frequency as well as the halving of the spatial wavelength of the
oscillations in the mode shape. These mode shapes show that, in the lowest air
split case, the main oscillation mode is the thermoacoustic oscillation. No
PVC is seen to be present in this case, and this is supported by the weak
recirculation zone seen in time-averaged velocity profiles in Fig. 2.

Figure 7: Radial mode shapes of the first mode for air split=0.2.

As the air split increases, the strongest peak changes frequency and the mode
shape of this new mode has a distinctively anti-symmetric spatial structure,
indicative of a PVC mode shape (Manoharan et al., 2020). While the dominant
peak now corresponds to the PVC mode, it is important to note that the
thermoacoustic peaks are still present, but significantly weaker than the PVC.
Figure 8 shows, for the case with the highest air split of 0.5, the radial
velocity mode shape of the strongest mode at a frequency of 762 Hz. This
spatial mode shape shows a strong precessing center at the base and then
helical vortical structures in the shear layer, both of which are
characteristic of a $m=1$ PVC mode. The peak seen at 586 Hz still corresponds
to a symmetric mode and so does the small peak seen at 918 Hz. Both these
peaks always exhibit a symmetric spatial structure, characteristic of $m=0$
thermoacoustic modes. The spatial modes at the peak frequencies show that in
the lower air split cases, the symmetric mode corresponding to thermoacoustic
oscillation is dominant, but as the air split increases, the thermoacoustic
modes are still present but the PVC gets stronger and dominates the flow
field.

Figure 8: Radial mode shapes of the first mode for air split=0.5.

The SPOD mode shapes are reconstructed at each of the frequency peaks in all
the cases and the peaks are characterized by the spatial structure of the
radial mode shape. From the mode shapes, it is seen that there are two
symmetric $m=0$ modes in the flow corresponding to frequencies near 600 Hz and
800 Hz, both of these modes are thermoacoustic oscillation modes. In the
lowest air split case, only the 600 Hz mode is seen, but subsequent cases show
the second thermoacoustic mode as well, where the 600 Hz mode is always
several orders of magnitude stronger than the 800 Hz thermoacoustic mode. The
PVC $m=1$ mode strengthens with increasing air split and its frequency
increases with linearly with increasing air split. The frequencies and
amplitudes of three oscillation modes from the SPOD for each condition are
depicted in the plot shown in Fig. 9.

Figure 9: Frequency (left) and amplitude (right) characterization of the three
oscillation modes seen at each air-split condition from SPOD.

From the plots in Fig. 9 and the eigenvalue spectra in Fig.6, two important
characteristics of the oscillation modes in the system can be observed. First,
while the frequency of the PVC increases linearly with increasing air-split,
the two symmetric frequencies remain largely the same at all air split
conditions; this is further confirmation that the 600 and 800 Hz modes are
thermoacoustic as they depend on the geometry of the system, not the flow
conditions. Second, increasing amplitude of the PVC corresponds to decreasing
amplitude of the first thermoacoustic mode; the second thermoacoustic mode is
weak at all conditions. A key observation of the experimental portion of this
study is that it appears that the global $m=1$ hydrodynamic instability, the
precessing vortex core, suppresses the $m=0$ flow field response to
thermoacoustic oscillation if the PVC amplitude is large enough. The
theoretical formulation described next explains the reason for this observed
behavior.

##  4 Theoretical Formulation

We formulate an asymptotic solution to the variable-density, low Mach number
Navier-Stokes equations to gain insight into the underlying mechanisms that
determine the coherent flow oscillation behavior observed in the experimental
results. The SPOD analysis presented in the previous section suggests that
acoustic forcing from the first thermoacoustic mode excites an axisymmetric
hydrodynamic response at a fixed frequency. The response is hydrodynamic
because it is driven by the dynamics of the vorticity field when perturbed by
acoustic velocity oscillations. Therefore, it is reasonable to expect that
these hydrodynamic oscillations are strong in regions of high time-averaged
vorticity. In regions where the time-averaged vorticity is small, flow
oscillations are simply acoustic velocity oscillations generated by the
acoustic pressure field. The discussion in prior sections has shown that the
amplitude of axisymmetric flow oscillations are strong in the flame region and
smaller in regions away from the flame. This shows that the thermoacoustic
pressure oscillation imposes a weak forcing on the time-averaged flow, where
‘weak’ implies that the amplitude of the hydrodynamic response is much larger
than that of the imposed forcing. This is common in acoustically-forced shear
flows (Mathews et al., 2016).

Manoharan et al. (2020) showed, using a high-$Re$, variable-$S$ jet, that PVC
oscillations are a result of linear instability at a critical swirl number
$S_{c}$ that induces a stable limit cycle oscillation in the flow,
characterized by helical precession of the vortex breakdown bubble around the
flow central axis. This precession, in turn, results in the helical
perturbation and rollup of the shear layers. A large number of prior
experiments support this observation (Oberleithner et al., 2011; Anacleto et
al., 2003; Escudier & Keller, 1985). The value of $S_{c}$ depends on the
evolution of the time-averaged flow field as $S$ varies. Therefore, in the
case of the present experiment, varying the air split between the inner and
outer flow circuits has the impact of varying $S_{c}$ while keeping $S$
nominally constant, as can be seen in Fig. 3 (right). Accordingly, for an air
split of 0.3 and above, Fig. 6 shows that $S>S_{c}$ due to the emergence of
the PVC.

Following Manoharan et al. (2020), we derive an asymptotic solution for the
flow state at a given $S$ assuming a small departure from $S_{c}$ as follows,

$S=S_{c}+\epsilon^{2}\Delta_{s}$ (2)

where $\epsilon$ is a small number and $\Delta_{s}\sim O(1)$. The governing
equations are the nonlinear, variable density Navier-Stokes (NS) equations in
the low Mach number limit, formulated in a cylindrical coordinate
$(r,\theta,z)$ system with the $z$ and $r$ axes aligned like the $x$ and $y$
axes, respectively, as shown in Fig. 1. Axial $(u_{z})$ and radial $(u_{r})$
velocity components are non-dimensionalized using a suitably chosen reference
velocity, $U_{z,ref}$. The azimuthal $(u_{\theta})$ velocity component is non-
dimensionalized using $U_{\theta,ref}=SU_{z,ref}$, where $S$ is the flow swirl
number. All lengths are normalized by the burner outer diameter. The governing
equations for the coherent flow component,
$\boldsymbol{\tilde{q}}(r,\theta,z,t)=\big{[}\tilde{\rho},\tilde{u}_{r},\tilde{u}_{\theta},\tilde{u}_{z},\tilde{p}\big{]}^{T}$,
can now be written as follows,

$\begin{split}\Bigg{(}\mathscr{B}+\mathscr{B}^{S}\\{\boldsymbol{\tilde{q}}\\}\Bigg{)}\frac{\partial\boldsymbol{\tilde{q}}}{\partial
t}=-\mathscr{N}\\{\boldsymbol{\tilde{q}}\\}\boldsymbol{\tilde{q}}-S\mathscr{N}^{S}\\{\boldsymbol{\tilde{q}}\\}\boldsymbol{\tilde{q}}-S^{2}\mathscr{N}^{SS}\\{\boldsymbol{\tilde{q}}\\}\boldsymbol{\tilde{q}}\\\
+\mathcal{L}_{T}\boldsymbol{\tilde{q}}+S\mathcal{L}_{T}^{S}\boldsymbol{\tilde{q}}+\epsilon^{3}\boldsymbol{\hat{q}_{a}}(r,z)\cos{(\omega_{a}t)}\end{split}$
(3)

where, the $O(\epsilon^{3})$ term on the right is the forcing term capturing
the weak forcing imposed by the axisymmetric thermoacoustic pressure
oscillation with a frequency $\omega_{a}$ and spatial amplitude distribution
$\boldsymbol{\hat{q}_{a}}$.

In eq. 3, the operators $\mathscr{B}$ and
$\mathscr{B}^{S}\\{\boldsymbol{\tilde{q}}\\}$ are diagonal matrices with
elements $\\{1,0,0,0,0\\}$ and $\\{0,\rho,S\rho,\rho,0\\}$ respectively and
$\mathscr{N}\\{\boldsymbol{\tilde{q}}\\}$,
$\mathscr{N}^{S}\\{\boldsymbol{\tilde{q}}\\}$ and
$\mathscr{N}^{SS}\\{\boldsymbol{\tilde{q}}\\}$ are nonlinear operators arising
from convective terms. $\mathcal{L}_{T}$ and $\mathcal{L}^{S}_{T}$ in eq. 3
are the linear operators arising from the pressure gradient term and molecular
and turbulent transport terms. An eddy viscosity model along with the
Boussinesq assumption is used to close the turbulent transport term as in
prior studies (Manoharan et al., 2020; Tammisola & Juniper, 2016; Oberleithner
et al., 2015). The various operators in eq. 3 are functions of the quantities
presented within the associated braces, $\\{~{}\\}$. Also, solutions to eq. 3
must be periodic in $\theta$ and obey no-slip conditions at the walls and
kinematic compatibility conditions on the flow centerline (Batchelor & Gill,
1962). The details of the all the operators used in eq. 3 are algebraically
complex expressions and are therefore presented in the appendix (eqs. 23, 26
and 27).

We note from the discussion in prior sections that the flow state is
characterized by flow oscillations with constant amplitudes. Accordingly, we
now derive a solution to eq. 3 accurate up to $O(\epsilon^{3})$, composed of
oscillatory components whose amplitudes grow slowly to their steady-state
oscillation values, using the method of multiple scales (Nayfeh, 2008). Thus,
we introduce a ‘fast’ time scale $t_{1}=t$, a ‘slow’ time scale
$t_{2}=\epsilon^{2}t$, and an expansion for $\boldsymbol{\tilde{q}}$ in terms
of $\epsilon$ as follows,

$\boldsymbol{\tilde{q}}(r,\theta,z,t_{1},t_{2})=\boldsymbol{q_{0}}(r,z)+\epsilon\boldsymbol{q_{1}}(r,\theta,z,t_{1},t_{2})+\epsilon^{2}\boldsymbol{q_{2}}(r,\theta,z,t_{1},t_{2})+\epsilon^{3}\boldsymbol{q_{3}}(r,\theta,z,t_{1},t_{2})+...$
(4)

where, $\boldsymbol{q_{0}}$ is the time-averaged flow state at $S_{c}$ and
each term in the expansion is assumed to depend independently on fast and slow
time scales. All the operators in eq. 3 are also expanded in powers of
$\epsilon$ (see eqs. 31, 36 and 41 for details). We then substitute these
expansions into eq. 3 and compare terms that are coefficients of the same
power of $\epsilon$ on both sides. This yields a sequence of equations for
$\boldsymbol{q_{1}}$, $\boldsymbol{q_{2}}$, $\dots$, etc. Although not shown
here, boundary conditions are treated similarly and yield homogeneous boundary
conditions that must be obeyed by the coefficients in eq. 4 at each order of
$\epsilon$.

At $O(\epsilon)$, the unsteady, linearized Navier-Stokes (LNS) equations for
$q_{1}$ are obtained as follows,

$\left(\mathscr{B}_{0}\frac{\partial}{\partial
t_{1}}+\mathscr{L}\right)\boldsymbol{q_{1}}=0$ (5)

where, $\mathscr{L}$ is an operator representing the spatial derivative terms
in the LNS equations and depends on $\boldsymbol{q_{0}}$ and $S_{c}$. The
details of the operators $\mathscr{B}_{0}$ and $\mathscr{L}$ are presented in
the appendix (eqs. 22 and 46 respectively). We decompose $\boldsymbol{q_{1}}$
into Fourier modes in the azimuthal direction and write the solution to eq. 5
as follows,

$\boldsymbol{q_{1}}(r,\theta,z,t_{1},t_{2})=A_{1}(t_{2})\boldsymbol{\hat{q}_{1}}(r,z)e^{i(\theta-\omega_{1}t_{1})}+A_{0}(t_{2})\boldsymbol{\hat{q}_{0}}(r,z)e^{-i\omega_{0}t_{1}}+c.c.$
(6)

where, for the present flow, contributions from an axisymmetric mode
($\boldsymbol{\hat{q}_{0}}$) and a helical mode ($\boldsymbol{\hat{q}_{1}}$)
have been retained. Thus, $\omega_{0}$ and $\omega_{1}$ are their associated
natural frequencies, which are given by the solution to the following
eigenvalue problem,

$\left(-i\omega_{m}\mathscr{B}_{0}+\mathscr{L}_{m}\right)\boldsymbol{\hat{q}_{m}}(r,z)=0$
(7)

where, the operator $\mathscr{L}_{m}$ is obtained from $\mathscr{L}$ by the
transformation $\partial/\partial\theta\rightarrow im$ and $m=0,1$ for the
axisymmetric and helical contributions in eq. 6, respectively. The fact that
the axisymmetric hydrodynamic mode oscillation frequency does not depend on
air split, as the data analysis in prior sections shows, reveals that this
mode is not self-excited, i.e., $\omega_{0i}<0$ \- the subscript ‘$i$’
represents the imaginary part. Accordingly, the thermoacoustic frequency,
$\omega_{a}$ can now be written in terms of $\omega_{0}$ and the difference
between $S$ and $S_{c}$ as follows,

$\omega_{a}=\omega_{0}+\epsilon^{2}\Omega_{a}$ (8)

where, $\Omega_{a}\sim O(1)$ and quantifies the extent of de-tuning between
the natural oscillation frequency $\omega_{0}$ and the thermoacoustic
frequency $\omega_{a}$. Next, since $S_{c}$ is the swirl number at which PVC
oscillations originate, $\omega_{1}$ is neutrally stable, i.e. $\omega_{1i}=0$
(Manoharan et al., 2020).

The solution for $\boldsymbol{q_{2}}$ in eq. 4 is determined from the equation
obtained by comparing $O(\epsilon^{2})$ terms on both sides of eq. 3 as
follows,

$\begin{split}\boldsymbol{q_{2}}(r,\theta,z,t_{1},t_{2})=\Delta_{s}\boldsymbol{\hat{q}_{\Delta}}+|A_{1}|^{2}\boldsymbol{\hat{q}_{A_{1}A_{1}^{*}}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\Big{(}A_{1}A_{0}\boldsymbol{\hat{q}_{A_{1}A_{0}}}e^{i\theta}e^{-i(\omega_{1}+\omega_{0})t_{1}}+A_{1}A^{*}_{0}\boldsymbol{\hat{q}_{A_{1}A_{0}^{*}}}e^{i\theta}e^{-i(\omega_{1}-\omega^{*}_{0})t_{1}}+c.c.\Big{)}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+\Big{(}(A_{1})^{2}\boldsymbol{\hat{q}_{A_{1}A_{1}}}e^{2i(\theta-\omega_{1}t_{1})}+(A_{0})^{2}\boldsymbol{\hat{q}_{A_{0}A_{0}}}e^{-2i\omega_{0}t_{1}}+c.c.\Big{)}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+\Big{(}B_{1}\boldsymbol{\hat{q}_{1}}e^{i(\theta-\omega_{1}t_{1})}+B_{0}\boldsymbol{\hat{q}_{0}}e^{-i\omega_{0}t_{1}}+c.c.\Big{)}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(9)

where, the first term quantifies the change in the time-averaged flow from its
state at $S_{c}$, the second term is the time-averaged distortion (to leading
order) due to finite amplitude helical oscillations. The third set of terms
are due to nonlinear coupling between the helical and axisymmetric modes. Note
that this results in helical oscillations at frequencies corresponding to the
sum and difference of the thermoacoustic and PVC frequencies. The fourth set
of terms are the first harmonic of the helical and axisymmetric oscillations.
The last set of terms are the complementary functions of the homogeneous form
of eq. 3 at $O(\epsilon^{2})$. These terms are included here for mathematical
completeness and will not contribute to the final steady state oscillatory
solution up to $O(\epsilon^{3})$.

The amplitudes $A_{0}(t_{2})$ and $A_{1}(t_{2})$ are determined from the
equation for $\boldsymbol{q_{3}}$ (coefficient of $\epsilon^{3}$ in eq. 4).
This equation has oscillatory source terms of the form
$e^{\\-i\omega_{0}t_{1}}$ and $e^{-i\omega_{1}t_{1}}$ which in general yield
terms in the solution for $q_{3}$ whose value grows exponentially with time.
Thus, a bounded solution is possible only when the coefficients of these
secular terms vanish. This condition yields evolution equations for
$A_{0}(t_{2})$ and $A_{1}(t_{2})$ as follows,

$\frac{dA_{0}}{dt_{2}}=\Delta_{s}\alpha_{A_{0}}A_{0}-\beta_{A_{0}A_{1}}|A_{1}|^{2}A_{0}+\beta_{A_{0}f}e^{-i\Omega_{a}t_{2}}$
(10)
$\frac{dA_{1}}{dt_{2}}=\Delta_{s}\alpha_{A_{1}}A_{1}-\beta_{A_{1}A_{1}}|A_{1}|^{2}A_{1}$
(11)

where, the coefficients in eqs. 10 and 11 are determined from inner products
between the adjoint modes $\boldsymbol{\hat{q}_{0}^{\dagger}}$ and
$\boldsymbol{\hat{q}_{1}^{\dagger}}$, and expressions involving various the
functions appearing in eqs. 6 and 9. These expressions are given in the
appendix (eqs. 57 \- 61). Equation 11 is equivalent to that derived by
Manoharan et al. (2020) for the evolution of helical mode oscillations in an
unforced flow with constant density. This shows that in the limit of weak
axisymmetric forcing, the evolution of $A_{1}$, i.e., helical oscillations, is
independent of $A_{0}$, i.e., the hydrodynamic response to the thermoacoustic
forcing. On the other hand, the second term in eq. 10 shows that helical
oscillations can influence the evolution of the axisymmetric hydrodynamic
response. The efficiency with which this occurs is determined by the value of
$\beta_{A_{0}A_{1}}$.

The steady-state oscillation amplitudes are determined as follows. First, we
introduce $A_{0}(t_{2})=D_{0}(t_{2})\exp[i\phi_{0}(t_{2})]$ and
$A_{1}(t_{2})=D_{1}(t_{2})\exp[i\phi_{1}(t_{2})]$ in eqs. 10 and 11, and
equate real and imaginary parts on both sides. The real parts define evolution
equations for $D_{1}$ and $D_{0}$. Requiring time derivatives of $D_{0}$ and
$D_{1}$ to vanish in these equations yields solutions for the steady-state
amplitudes. The imaginary parts then yields expressions for the natural
frequencies of the axisymmetric and helical modes at $S$ to leading order in
$S-S_{c}$.

The solutions for the steady-state amplitude of the helical mode and its
characteristic frequency are as follows,

$A_{PVC}=\sqrt{\frac{(S-S_{c})\alpha_{A_{1}r}}{\beta_{A_{1}A_{1}r}}}$ (12)

and,

$\omega_{PVC}=\omega_{1}+(S-S_{c})\bigg{(}\alpha_{A_{1}r}\frac{\beta_{A_{1}A_{1}i}}{\beta_{A_{1}A_{1}r}}-\alpha_{A_{1}i}\bigg{)}$
(13)

where, the letters ‘$r$’ and ‘$i$’ in the subscripts denote real and imaginary
parts. Note that as may be expected, the results in eqs. 12 and 13 are
analogous to those derived by Manoharan et al. (2020) for a constant density
flow.

Manoharan et al. (2020) show in the case of their variable-swirl, non-reacting
jet that helical oscillations occur due to the inception of vortex breakdown
at $S=S_{c}$, with the breakdown bubble being increasingly well-defined
thereafter. For the present case, $S_{c}$ quantifies the minimum swirl
intensity that the nozzle needs to generate for a given air split in order to
induce vortex breakdown in the flow. Thus, the appearance of a well-defined
vortex breakdown bubble on the flow centreline with increasing air split as
Fig. 3 (left) shows, with $S$ remaining nearly the same across air splits (see
Fig. 3 (right)), is a consequence of a corresponding reduction in $S_{c}$.
Therefore, eq. 13 suggests the linear rise in characteristic PVC frequency
with air split, as Fig. 9 (left) shows, is due to the increase in $S-S_{c}$.
Likewise, eq. 12 then shows that the increase in PVC amplitude seen in Fig. 9
(right) is because of the same reason. The slight decrease in the amplitude
for an air split of 0.5 in Fig. 9 (right) is possibly due to the fact that the
bubble moves upstream into the nozzle, resulting in a part of the oscillation
field moving out of the field of view of the sPIV cameras. This blocked view
causes only a part of the total energy of the oscillation to be captured by
the SPOD for this case alone. Note that the characteristic PVC frequency
continues to increase linearly even for this case; see Fig. 9 (left).

For the axisymmetric mode, the steady-state oscillation frequency is simply
the thermoacoustic frequency $\omega_{a}$. The steady state amplitude,
$A_{TH}$ of the axisymmetric mode is given by the following,

$A_{TH}=\frac{|\beta_{A_{0}f}|}{\Big{|}i\big{[}\omega_{a}-\omega_{0}-i(S-S_{c})\alpha_{A_{0}}\big{]}-\beta_{A_{0}A_{1}}A_{PVC}^{2}\Big{|}}$
(14)

where $\beta_{A_{0}f}$ is given as follows:

$\beta_{A_{0}f}=\frac{\langle\boldsymbol{\hat{q}_{0}^{\dagger}},\boldsymbol{\hat{q}_{a}}\rangle}{2\langle\boldsymbol{\hat{q}_{0}^{\dagger}},\mathscr{B}_{0}\boldsymbol{\hat{q}_{0}}\rangle}$
(15)

and $\beta_{A_{0}A_{1}}$ is given by the following expression:

$\begin{split}\beta_{A_{0}A_{1}}=\frac{1}{\langle\boldsymbol{\hat{q}_{0}^{\dagger}},\mathscr{B}_{0}\boldsymbol{\hat{q}_{0}}\rangle}\bigg{\langle}\boldsymbol{\hat{q}^{0\dagger}},\Big{[}\mathbb{I}_{(-1,1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}^{*}}\\}\boldsymbol{\hat{q}_{A_{1}A_{0}}}+\mathbb{I}_{(1,-1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{A_{1}A_{0}}}\\}\boldsymbol{\hat{q}_{1}^{*}}+\mathbb{I}_{(-1,1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}\\}\boldsymbol{\hat{q}_{1}}\\\
+~{}\mathbb{I}_{(1,-1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}+\mathbb{I}_{(0,0)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{A_{1}A_{1}^{*}}}\\}\boldsymbol{\hat{q}_{0}}+\mathbb{S}_{(0,0)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}}\\}\boldsymbol{\hat{q}_{A_{1}A_{1}^{*}}}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}\mathbb{Q}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{\hat{q}_{0}}+\mathbb{Q}_{1}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{0}}\\}\boldsymbol{\hat{q}_{1}}+\mathbb{Q}_{-1}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{0}}\\}\boldsymbol{\hat{q}_{1}^{*}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+\mathbb{T}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{A_{1}A_{1}^{*}}},\boldsymbol{\hat{q}_{A_{1}A_{0}}},\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}\\}\boldsymbol{q_{0}}\Big{]}\bigg{\rangle}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(16)

where, the inner product between two functions $\boldsymbol{\hat{p}_{1}}(r,z)$
and $\boldsymbol{\hat{p}_{2}}(r,z)$ is defined as follows,

$\langle\boldsymbol{\hat{p}_{1}},\boldsymbol{\hat{p}_{2}}\rangle=\int_{0}^{Z}\int_{0}^{R}(\boldsymbol{\hat{p}_{1}})^{H}\boldsymbol{\hat{p}_{2}}\,r\,dr\,dz$
(17)

where, the superscript ‘$H$’ denotes the transpose conjugate.

In eq. 14, $\beta_{A_{0}f}$ quantifies the receptivity of the axisymmetric
mode to the forcing imposed by the thermoacosutic oscillations. The
denominator in eq. 14 quantifies two separate sources of de-tuning. The first
term in box brackets in the denominator of eq. 14 is the difference between
the thermoacoustic forcing and the natural frequency of the axisymmetric mode
at $S$. Thus, as may be expected, a larger difference between these two
frequencies results in a smaller response amplitude. The last term in
denominator of eq. 14 shows that the presence of helical oscillations can
additionally reduce $A_{TH}$ if $|\beta_{A_{0}A_{1}}|$ is large, i.e., the
nonlinear coupling between helical and axisymmetric oscillations is large.
From the expression for $\beta_{A_{0}A_{1}}$ in eq. 16, the contribution from
nonlinear coupling between the helical and axisymmetric mode depends on
$\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}$ and
$\boldsymbol{\hat{q}_{A_{1}A_{0}}}$. These quantities are components of the
solution for $\boldsymbol{q_{2}}$ (see eq. 9) that oscillate with frequencies
corresponding to the difference and sum of the helical and axisymmetric
oscillations. The equation that determines
$\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}$ is as follows,

$\begin{split}\Big{(}i\Delta\omega\mathscr{B}_{0}+\mathscr{L}_{-1}\Big{)}\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}=-\Big{(}\mathbb{I}_{(0,-1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}}\\}\boldsymbol{\hat{q}_{1}^{*}}+\mathbb{I}_{(-1,0)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}^{*}}\\}\boldsymbol{\hat{q}_{0}}+\mathbb{Q}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{1}^{*}}\\}\boldsymbol{q_{0}}\Big{)}\end{split}$
(18)

where, $\Delta\omega=\omega_{1}-\omega_{0}$ and the various operators on the
right are given in the appendix (eqs. 67 and 66). The solution to eq 18 can be
written in terms of the eigenfunctions of $\mathscr{L}_{-1}$,
$\boldsymbol{\hat{q}_{-1,k}}$, and their adjoints as follows,

$\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}={\sum}_{k}\Bigg{[}\frac{\Big{\langle}\boldsymbol{\hat{q}^{\dagger}_{-1,k}},i\Big{(}\mathbb{I}_{(0,-1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}}\\}\boldsymbol{\hat{q}_{1}^{*}}+\mathbb{I}_{(-1,0)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}^{*}}\\}\boldsymbol{\hat{q}_{0}}+\mathbb{Q}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{1}^{*}}\\}\boldsymbol{q_{0}}\Big{)}\Big{\rangle}}{(\Delta\omega+\omega_{-1,k})\langle\boldsymbol{\hat{q}^{\dagger}_{-1,k}},\mathscr{B}_{0}\boldsymbol{\hat{q}_{-1,k}}\rangle}\boldsymbol{\hat{q}_{-1,k}}\Bigg{]}$
(19)

where, the summation is over the entire eigenspace of $\mathscr{L}_{-1}$.
Thus, eq. 19 shows that when $\Delta\omega$ is small, i.e., when the
axisymmetric and helical mode natural frequencies are close, the magnitude of
$\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}$ is large. Therefore, when the forcing
is near-resonant, i.e., $\omega_{a}\sim\omega_{0}$, $\beta_{A_{0}A_{1}}$
becomes large when $\omega_{1}\sim\omega_{a}$, resulting in an efficient
reduction in the amplitude of the axisymmetric hydrodynamic response as eq. 14
suggests and Fig. 9 (right) shows for the present experiment.

Thus, the asymptotic solution for the stationary flow state, up to
$\epsilon^{3}$, can be written by combining eqs. 6, 9 12, 13 and 14 as
follows,

$\begin{split}\boldsymbol{\tilde{q}}(r,\theta,z,t)=\boldsymbol{q_{0}}+(S-S_{c})\boldsymbol{q_{\Delta}}+A_{PVC}^{2}\boldsymbol{q_{A_{1}A^{*}_{1}}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+\Big{(}A_{PVC}\boldsymbol{\hat{q}_{1}}e^{i\theta}e^{-i\omega_{PVC}t}+A_{TH}\boldsymbol{\hat{q}_{0}}e^{-i\omega_{a}t}+c.c.\Big{)}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+\Big{(}A_{PVC}^{2}\boldsymbol{\hat{q}_{A_{1}A_{1}}}e^{2i\theta}e^{-2i\omega_{PVC}t}+A_{TH}^{2}\boldsymbol{\hat{q}_{A_{0}A_{0}}}e^{-2i\omega_{a}t}+c.c.\Big{)}~{}~{}~{}~{}\\\
+\Big{(}A_{PVC}A_{TH}\boldsymbol{\hat{q}_{A_{1}A_{0}}}e^{i\theta}e^{-i(\omega_{PVC}+\omega_{a})t}+c.c.\Big{)}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+\Big{(}A_{PVC}A_{TH}\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}e^{-i\theta}e^{i(\omega_{PVC}-\omega_{a})t}+c.c.\Big{)}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+O\big{(}(S-S_{c})^{3/2}\big{)})~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(20)

where, $t_{1}$ and $t_{2}$ have been replaced by $t$. It is now evident from
eq. 20 that in general, the flow solution must have oscillating components at
the thermoacoustic frequency, the PVC frequency, and their sum and difference.
This result is consistent with the SPOD spectra for the velocity field shown
in Fig. 6. The suppression of the axisymmetric hydrodynamic response to
acoustic forcing, due to the presence of a PVC, can potentially result in the
suppression of the component of global heat release oscillation due to burning
area oscillations in the case of an axisymmetric flame (Acharya et al., 2012;
Moeck et al., 2012). Therefore, ensuring that the characteristic frequency of
PVC oscillations in a combustor nozzle matches those of potentially unstable
thermoacoustic modes during initial stages of design might prove beneficial
because the presence of a suitably designed PVC could suppress global heat
release rate oscillations through the nonlinear mechanism described in this
section as the qualitative agreement between flow dynamics in the present
experimental study and the theory shows.

##  5 Impact of hydrodynamics on flame behavior

The time-averaged fuel and flame contours, such as those shown in Fig. 4, do
not show significant variations with air split. The instantaneous images,
however, distinctly show a shift from largely symmetric to asymmetric with
increasing air split. Figure 10 shows instantaneous images of the flame,
obtained from the OH-PLIF images (depicted by the red contour), and the fuel
distribution, obtained from the acetone-PLIF images (depicted by the green
contour), at a low air split of 0.25 and a high air split of 0.5.

Figure 10: Instantaneous images of superimposed OH-PLIF (red) and acetone-PLIF
(green) profiles at air splits 0.25 (left) and 0.5 (right). Videos of these
conditions are provided in the supplemental material.

The instantaneous images illustrate the transition from symmetric motion to
anti-symmetric motion as the air split is increased. In the low air split
case, symmetric vortex roll-up of the flame surface can be seen. The fuel jets
interact with the symmetric vortical structures, leading up to the onset of
thermoacoustic instability as a result of periodic fluctuations in both the
air and fuel flow rates. The high air split case, on the other hand, shows an
anti-symmetric structure, which resembles the helical vortex structures
generated by the PVC. These instantaneous images illustrate the influence that
the PVC has on the flame structure and fuel distribution in the system. In the
absence of the PVC, the symmetric thermoacoustic mode governs the flame
response, but as the PVC amplitude increases, the PVC begins to govern the
fuel flow and flame response. This transition would suggest that the flame
response is strongly coupled to the velocity field in the flow, indicating
that in the presence of a sufficiently strong global hydrodynamic instability,
the flame is less receptive to thermoacoustic oscillations.

The results from the frequency-domain characterization of the flow field show
that as the air split is increased, the PVC strengthens and eventually
dominates the flow field dynamics. The instantaneous flame structure and fuel
distribution profiles indicate that the flame and fuel response are coupled
with the velocity response. In order to draw conclusions about the impact of
the PVC on the fuel/air mixing and flame dynamics, a more quantitative
analysis is needed.

We perform spectral proper orthogonal decomposition on the acetone- and OH-
PLIF image time series. Figure 11 shows the eigenvalue spectra of the flame
oscillation, obtained from the OH-PLIF images, at each air split. The peak
frequencies in the OH-PLIF spectra are similar to the peaks seen in the
velocity spectra, where the thermoacoustic modes and PVC are indicated with
the same symbols as in Figs. 5 and 6. In the low air split cases, the
thermoacoustic mode dominates the flame behavior as the first thermoacoustic
mode amplitude is the highest peak. As the air split is increased, the PVC
peak grows in strength, indicating that the PVC fluctuations drive the flame
behavior. Furthermore, similar to the flow field response, whenever the PVC
amplitude is high, the thermoacoustic mode is suppressed. The eigenvalue
spectra of the OH-PLIF fields indicate that the flame is strongly coupled with
the flow field.

Figure 11: Eigenvalue spectra of OH-PLIF field obtained from spectral POD.

Figure 12 shows the SPOD spectra of the acetone-PLIF images. Compared to the
eigenvalue spectra of the velocity fields, the eigenvalue spectra of the PLIF
data shows more peaks. This is likely because the acetone-PLIF data is
inherently more noisy as compared to the velocity data. Unlike the OH-PLIF
eigenvalue spectra, the acetone-PLIF eigenvalue spectra are distinctly
different from the velocity field. While the peak frequencies are largely the
same, the relative peak amplitudes are significantly different. Just like in
the flame and velocity spectra, as the air split increases, the strength of
the PVC mode increases. However, in the acetone spectra, the thermoacoustic
mode is just as strong, if not stronger, than the PVC mode, even in the high
air split cases, indicating that, unlike flame surface oscillation, fuel
fluctuation is not dominated purely by the PVC. This is likely due to the
distance between the base of the recirculation zone and the location of fuel
injection at the dump plane, so the PVC does not impact fuel injection. As
seen in Fig. 3, the center of the recirculation zone is located at a
downstream location of 10 mm, whereas the fuel is injected at the dump plane.
The fuel fluctuation, therefore, is impacted more by the pressure
oscillations, which likely have an anti-node at the dump plane, than the
oscillation of the PVC.

Figure 12: Eigenvalue spectra of the acetone field obtained from spectral POD.

The eigenvalue spectra of the acetone-PLIF images suggest that the fuel
fluctuation is not as sensitive to the PVC as the flame fluctuation. However,
since the spectral POD method considers the energy of the entire flow field,
it is important to be careful when interpreting the results of the eigenvalue
decomposition. While the acetone burns out fairly close to the nozzle, the LIF
signal further downstream could be signal from either flame or soot
luminosity. In order to correctly interpret the results of the spectral
characterization, the oscillation modes at different spatial locations in the
flow need to be observed. To that end, we perform a Fourier transform on the
time series of the PLIF images at every pixel location in the field of view.
At each location, the peak frequency is extracted and plotted as a contour in
space.

Figure 13 shows the peak frequency contours of the OH field. The white contour
represents the $v_{z,mean}=0$ contour, which is indicative of the edge
recirculation zone; this inner shear layer is where the flame stabilizes in
this configuration. The spatial peak frequency contours of the OH field, much
like its SPOD spectra, show that the OH field is largely driven by the PVC
once the PVC is of sufficient strength. The peak frequency contours of the
velocity fields (not shown) are very similar to that of the OH field.

Figure 13: Peak frequency contours of OH field. The white contour represents
the recirculation zone border, ($v_{z,mean}=0$). Figure 14: Peak frequency
contours of acetone field. The white contour represents the recirculation zone
border, ($v_{z,mean}=0$).

Figure 14 shows the peak frequency contours of the acetone field. The acetone
contours show a significantly different spatial distribution in the peak
frequency contours compared to the frequency contours of the OH and velocity
fields. While the PVC oscillation still dominates regions close to the shear
layer, these regions are substantially smaller when compared to the peak
frequency contours of the OH field. The results from this spatial frequency
characterization support the results obtained from the eigenvalue spectra
obtained using SPOD. As air split increases, the PVC strengthens and
progressively dominates the flame and flow response. However, it is seen that,
due to the relative offset between the downstream locations of the
recirculation zone and the fuel injection location, the fuel fluctuation is
less sensitive to the PVC oscillation when compared to the flame and flow
fields.

These differences in dominant frequency have implications for thermoacoustic
coupling processes. In this complex, realistic combustor flow field, operating
at an elevated pressure, the thermoacoustic feedback loop could involve both
velocity coupling and equivalence ratio coupling mechanisms. As such, it is
important to characterize not only how these coupling mechanisms evolve with
changing flow conditions but also how the two coupling pathways interact with
one another. The results presented here indicate that the hydrodynamic
instability and equivalence ratio coupling pathways do not strongly interact
with each other since the recirculation zone is located sufficiently further
downstream of the fuel injection location.

##  6 Conclusions

This study showed using both experimental and theoretical analysis that a
precessing vortex core can suppress longitudinal thermoacoustic instability in
a partially-premixed, swirl-stabilized flame at elevated pressure. When the
PVC frequency, which increases with air split, approaches that of the
thermoacoustic instability, a nonlinear coupling, shown by eq. 10, between the
$m=1$ PVC motion and the $m=0$ shear layer response leads to suppression of
the $m=0$ mode and hence the velocity-coupled thermoacoustic instability.

In addition to this main conclusion, this work has identified a number of new
results in the PVC literature. First, varying the air split in this injector
does not change the swirl number of the flow as the geometric swirl number of
the inner and outer swirl passages are relatively similar. However, increases
in air split result in stronger central recirculation in the vortex breakdown
bubble, although the location of the bubble does not change significantly with
air split. This result is explained by the variation in the characteristic
critical swirl number of the burner that results from changes in air split.
This critical swirl number is the minimum swirl number needed at a given air
split to initiate vortex breakdown in the flow. The results suggest that
increasing the air split decreases the critical swirl number as a result of
structural changes to the flow field. As the offset between the operating
swirl number of the experiment and critical swirl number increases, the vortex
breakdown bubble recirculation strengthens. This result has been seen in
previous studies by increasing swirl number (Liang & Maxworthy, 2005; Escudier
& Keller, 1985), but, to the best of the authors’ knowledge, not seen by
decreasing the critical swirl number at constant swirl.

Second, the impact of the PVC on shear layer behavior is consistent with many
previous studies (Oberleithner et al., 2011; Frederick et al., 2018; Moeck et
al., 2012), where the shear layers surrounding the vortex breakdown region
respond to the wavemaker region at the base of the recirculation zone driven
by the $m=1$ precession mode of the breakdown bubble. The spatial plots of
peak frequency in Fig. 14 show, however, that the influence of the PVC is
localized and does not drive fluctuations in the fuel at the dump plane, which
is located at a distance from the upstream stagnation point of the
recirculation zone. This result indicates that the influence of the PVC is
local, unlike the dynamics of the thermoacoustic mode, which are felt
throughout the combustor as they are driven by a pressure mode shape inside
the combustor.

Finally, the insights obtained from this analysis have significant
implications for the suppression of thermoacoustic oscillations in gas turbine
combustors. Work by first Mathews et al. (2016) and then Frederick et al.
(2018) showed, in a non-reacting, variable-swirl-number jet configuration,
that the presence of the PVC suppressed the response of the shear layer to
external forcing. Calculation of the spatial growth rate of the $m=0$ shear
layer instability in Frederick et al. (2018) indicated that the PVC caused the
growth rate to become negative, resulting in no response at the forcing
frequency. We had suggested, at the time, that these results may be extended
to reacting circumstances as long as the PVC could be sustained. The present
results confirm this hypothesis by showing the nonlinear coupling mechanism
that leads to $m=0$ mode suppression, and hence suppression of the velocity-
coupling pathway for thermoacoustic instability. This result suggests that
“designing” a PVC in a combustor flow field at frequencies close to those of a
longitudinal thermoacoustic mode could suppress oscillations.

##  7 Acknowledgements

This work was supported by National Science Foundation grant CBET-1749679 and
Air Force Office of Scientific Research grant FA9550-16-1-0044. Any opinions,
findings, and conclusions or recommendations expressed in this material are
those of the authors and do not necessarily reflect the views of the National
Science Foundation or the Air Force Office of Scientific Research. The
assistance of Dr. Jisu Yoon (ADD, Republic of Korea), Dr. Carrie Noren (AFRL,
USA) and Dr. Klaus-Peter Geigle (DLR, Germany) with the setup and operation of
the experimental apparatus described in this paper are gratefully
acknowledged.

##  8 Declaration of Interests

The authors report no conflict of interest.

##  Appendix A Governing equations in the operator form

The governing equations for
$\boldsymbol{\tilde{q}}=\big{[}\tilde{\rho},\tilde{u}_{r},\tilde{u}_{\theta},\tilde{u}_{z},\tilde{p}\big{]}^{T}$
i.e. the coherent flow components are represented in the operator form as
given in Eq. 3. The operators $\mathscr{B}$ and
$\mathscr{B}^{S}\\{\boldsymbol{\tilde{q}}\\}$ are diagonal matrices with
elements
$\mathscr{B}=\begin{bmatrix}1~{}~{}0~{}~{}0~{}~{}0~{}~{}0\end{bmatrix}$ and
$\mathscr{B}^{s}\\{\boldsymbol{\tilde{q}}\\}=\begin{bmatrix}0~{}~{}\tilde{\rho}~{}~{}S\tilde{\rho}~{}~{}\tilde{\rho}~{}~{}0\end{bmatrix}$.
The operator $\mathscr{B}^{S}\\{\boldsymbol{\tilde{q}}\\}$ can then be
expanded as follows:

$\mathscr{B}^{S}\\{\boldsymbol{\tilde{q}}\\}=\mathscr{B}^{S}_{1}\\{\boldsymbol{\tilde{q}}\\}+\epsilon^{2}\Delta_{s}\mathscr{B}^{S}_{2}\\{\boldsymbol{\tilde{q}}\\}$
(21)

where, $\mathscr{B}^{S}_{1}\\{\boldsymbol{\tilde{q}}\\}$ and
$\mathscr{B}^{S}_{2}\\{\boldsymbol{\tilde{q}}\\}$ are diagonal matrices with
elements
$\mathscr{B}^{S}_{1}\\{\boldsymbol{\tilde{q}}\\}=\begin{bmatrix}0~{}~{}\tilde{\rho}~{}~{}S_{c}\tilde{\rho}~{}~{}\tilde{\rho}~{}~{}0\end{bmatrix}$
and
$\mathscr{B}^{S}_{2}\\{\boldsymbol{\tilde{q}}\\}=\begin{bmatrix}0~{}~{}0~{}~{}\tilde{\rho}~{}~{}0~{}~{}0\end{bmatrix}$.
The matrix $\mathscr{B}_{0}$ is defined as:

$\mathscr{B}_{0}=\mathscr{B}+\mathscr{B}^{S}_{1}\\{\boldsymbol{q_{0}}\\}$ (22)

The nonlinear operators $\mathscr{N}\\{\boldsymbol{\tilde{q}}\\}$,
$\mathscr{N}^{S}\\{\boldsymbol{\tilde{q}}\\}$ and
$\mathscr{N}^{SS}\\{\boldsymbol{\tilde{q}}\\}$ in eq. 3 acting on the vector
field $\boldsymbol{\tilde{q}}$, are defined as follows:

$\mathscr{N}\\{\boldsymbol{\tilde{q}}\\}=\mathscr{N}_{1}\\{\boldsymbol{\tilde{q}}\\}+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{\tilde{q}}\\}\mathscr{N}_{2}\\{\boldsymbol{\tilde{q}}\\}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\tilde{q}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\tilde{q}}\\}\mathscr{N}_{3}\\{\boldsymbol{\tilde{q}}\\}$
(23)

where,

$\mathscr{N}_{1}\\{\boldsymbol{\tilde{q}}\\}=\begin{bmatrix}\tilde{\Delta}_{rz}&&\tilde{\rho}\left(\frac{1}{r}+\frac{\partial}{\partial
r}\right)&&0&&\tilde{\rho}\frac{\partial}{\partial z}&&0\\\ \\\
0&&0&&0&&0&&0\\\ \\\ 0&&0&&0&&0&&0\\\ \\\ 0&&0&&0&&0&&0\\\ \\\
\frac{1}{RePr}\left[\tilde{\rho}\Lambda_{r\theta
z}-2\left(\frac{\partial\tilde{\rho}}{\partial r}\frac{\partial}{\partial
r}+\frac{1}{r^{2}}\frac{\partial\tilde{\rho}}{\partial\theta}\frac{\partial}{\partial\theta}+\frac{\partial\tilde{\rho}}{\partial
z}\frac{\partial}{\partial z}\right)\right]&&0&&0&&0&&0\end{bmatrix}$ (24)
$\mathscr{N}_{2}\\{\boldsymbol{\tilde{q}}\\}=\begin{bmatrix}0&&0&&0&&0&&0\\\
\\\ 0&&\tilde{\Delta}_{rz}&&0&&0&&0\\\ \\\ 0&&0&&0&&0&&0\\\ \\\
0&&0&&0&&\tilde{\Delta}_{rz}&&0\\\ \\\
0&&0&&0&&0&&0\end{bmatrix}~{}~{}~{}~{}~{}\mathscr{N}_{3}\\{\boldsymbol{\tilde{q}}\\}=\begin{bmatrix}0&&0&&0&&0&&0\\\
\\\ 0&&0&&0&&0&&0\\\ \\\ 0&&0&&0&&0&&0\\\ \\\ 0&&0&&0&&0&&0\\\ \\\
0&&\tilde{\rho}\left(\frac{\partial}{\partial
r}+\frac{1}{r}\right)&&0&&\tilde{\rho}\frac{\partial}{\partial
z}&&0\end{bmatrix}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$ (25)

Similarly, we can write $\mathscr{N}^{S}\\{\boldsymbol{\tilde{q}}\\}$ and
$\mathscr{N}^{SS}\\{\boldsymbol{\tilde{q}}\\}$ as follows:

$\mathscr{N}^{S}\\{\boldsymbol{\tilde{q}}\\}=\mathscr{N}^{S}_{1}\\{\boldsymbol{\tilde{q}}\\}+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{\tilde{q}}\\}\mathscr{N}^{S}_{2}\\{\boldsymbol{\tilde{q}}\\}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\tilde{q}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\tilde{q}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{\tilde{q}}\\}$
(26)
$\mathscr{N}^{SS}\\{\boldsymbol{\tilde{q}}\\}=\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{\tilde{q}}\\}\mathscr{N}^{SS}_{2}\\{\boldsymbol{\tilde{q}}\\}$
(27)

where,

$\mathscr{N}^{S}_{1}\\{\boldsymbol{\tilde{q}}\\}=\begin{bmatrix}\frac{\tilde{u}_{\theta}}{r}\frac{\partial}{\partial\theta}&&0&&\frac{\tilde{\rho}}{r}\frac{\partial}{\partial\theta}&&0&&0\\\
\\\ 0&&0&&0&&0&&0\\\ \\\ 0&&0&&0&&0&&0\\\ \\\ 0&&0&&0&&0&&0\\\ \\\
0&&0&&0&&0&&0\end{bmatrix}~{}~{}~{}~{}~{}~{}~{}~{}\mathscr{N}^{S}_{3}\\{\boldsymbol{\tilde{q}}\\}=\begin{bmatrix}0&&0&&0&&0&&0\\\
\\\ 0&&0&&0&&0&&0\\\ \\\ 0&&0&&0&&0&&0\\\ \\\ 0&&0&&0&&0&&0\\\ \\\
0&&0&&\frac{\tilde{\rho}}{r}\frac{\partial}{\partial\theta}&&0&&0\end{bmatrix}$
(28)
$\mathscr{N}^{S}_{2}\\{\boldsymbol{\tilde{q}}\\}=\begin{bmatrix}0&&0&&0&&0&&0\\\
\\\ 0&&\frac{\tilde{u}_{\theta}}{r}\frac{\partial}{\partial\theta}&&0&&0&&0\\\
\\\ 0&&0&&\big{(}\tilde{\Delta}_{rz}+\frac{\tilde{u}_{r}}{r}\big{)}&&0&&0\\\
\\\ 0&&0&&0&&\frac{\tilde{u}_{\theta}}{r}\frac{\partial}{\partial\theta}&&0\\\
\\\ 0&&0&&0&&0&&0\end{bmatrix}$ (29)
$\mathscr{N}^{SS}_{2}\\{\boldsymbol{\tilde{q}}\\}=\begin{bmatrix}0&&0&&0&&0&&0\\\
\\\ 0&&0&&-\frac{\tilde{u}_{\theta}}{r}&&0&&0\\\ \\\
0&&0&&\frac{\tilde{u}_{\theta}}{r}\frac{\partial}{\partial\theta}&&0&&0\\\ \\\
0&&0&&0&&0&&0\\\ \\\ 0&&0&&0&&0&&0\end{bmatrix}$ (30)

where we have defined $\Lambda_{r\theta z}=\left(\frac{\partial^{2}}{\partial
r^{2}}+\frac{1}{r}\frac{\partial}{\partial
r}+\frac{1}{r^{2}}\frac{\partial^{2}}{\partial\theta^{2}}+\frac{\partial^{2}}{\partial
z^{2}}\right)$ and
$\tilde{\Delta}_{rz}=\left(\tilde{u}_{r}\frac{\partial}{\partial
r}+\tilde{u}_{z}\frac{\partial}{\partial z}\right)$. Also, we define matrix
$\mathscr{R}^{\rho}_{i,j,...}\\{\boldsymbol{\tilde{q}}\\}$ as a diagonal
matrix with non-zero entries as $\tilde{\rho}$ on the $i^{th},j^{th},...$
indices. For example, $\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{\tilde{q}}\\}$
is a diagonal matrix with elements:
$\begin{bmatrix}0&&\tilde{\rho}&&0&&\tilde{\rho}&&0\end{bmatrix}$.

The nonlinear operators can be expanded in powers of $\epsilon$ as follows:

$\mathscr{N}\\{\boldsymbol{\tilde{q}}\\}=\mathbb{N}^{1}\\{\boldsymbol{q_{0}}\\}+\epsilon\mathbb{N}^{\epsilon}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}}\\}+\epsilon^{2}\mathbb{N}^{\epsilon^{2}}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}},\boldsymbol{q_{2}}\\}+\epsilon^{3}\mathbb{N}^{\epsilon^{3}}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}},\boldsymbol{q_{2}},\boldsymbol{q_{3}}\\}+...$
(31)

where,

$\mathbb{N}^{1}\\{\boldsymbol{q_{0}}\\}=\mathscr{N}_{1}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{q_{0}}\\}\mathscr{N}_{2}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{0}}\\}$
(32)

$\begin{split}\mathbb{N}^{\epsilon}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}}\\}=\mathscr{N}_{1}\\{\boldsymbol{q_{1}}\\}+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{q_{1}}\\}\mathscr{N}_{2}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{q_{0}}\\}\mathscr{N}_{2}\\{\boldsymbol{q_{1}}\\}~{}~{}~{}~{}~{}~{}~{}\\\
+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{1}}\\}\end{split}$
(33)

$\begin{split}\mathbb{N}^{\epsilon^{2}}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}},\boldsymbol{q_{2}}\\}=\mathscr{N}_{1}\\{\boldsymbol{q_{2}}\\}+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{q_{2}}\\}\mathscr{N}_{2}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{q_{0}}\\}\mathscr{N}_{2}\\{\boldsymbol{q_{2}}\\}~{}~{}~{}~{}~{}~{}~{}\\\
+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{2}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{2}}\\}\\\
+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{q_{1}}\\}\mathscr{N}_{2}\\{\boldsymbol{q_{1}}\\}+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{1}}\\}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{0}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(34)

$\begin{split}\mathbb{N}^{\epsilon^{3}}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}},\boldsymbol{q_{2}},\boldsymbol{q_{3}}\\}=\mathscr{N}_{1}\\{\boldsymbol{q_{3}}\\}+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{q_{3}}\\}\mathscr{N}_{2}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{q_{0}}\\}\mathscr{N}_{2}\\{\boldsymbol{q_{3}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{3}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{3}}\\}~{}~{}\\\
+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{q_{1}}\\}\mathscr{N}_{2}\\{\boldsymbol{q_{2}}\\}+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{2}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{2}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{q_{2}}\\}\mathscr{N}_{2}\\{\boldsymbol{q_{1}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{2}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{1}}\\}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{2}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{0}}\\}~{}~{}\\\
+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{1}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(35)

Similarly,

$\mathscr{N}^{S}\\{\boldsymbol{\tilde{q}}\\}=\mathbb{N}^{S,1}\\{\boldsymbol{q_{0}}\\}+\epsilon\mathbb{N}^{S,\epsilon}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}}\\}+\epsilon^{2}\mathbb{N}^{S,\epsilon^{2}}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}},\boldsymbol{q_{2}}\\}+\epsilon^{3}\mathbb{N}^{S,\epsilon^{3}}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}},\boldsymbol{q_{2}},\boldsymbol{q_{3}}\\}+...$
(36)

where,

$\mathbb{N}^{S,1}\\{\boldsymbol{q_{0}}\\}=\mathscr{N}^{S}_{1}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{S}_{2}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{0}}\\}$
(37)

$\begin{split}\mathbb{N}^{S,\epsilon}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}}\\}=\mathscr{N}^{S}_{1}\\{\boldsymbol{q_{1}}\\}+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{q_{1}}\\}\mathscr{N}^{S}_{2}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{S}_{2}\\{\boldsymbol{q_{1}}\\}~{}~{}\\\
+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{1}}\\}\end{split}$
(38)

$\begin{split}\mathbb{N}^{S,\epsilon^{2}}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}},\boldsymbol{q_{2}}\\}=\mathscr{N}^{S}_{1}\\{\boldsymbol{q_{2}}\\}+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{q_{2}}\\}\mathscr{N}^{S}_{2}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{S}_{2}\\{\boldsymbol{q_{2}}\\}~{}~{}\\\
+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{2}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{2}}\\}\\\
+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{q_{1}}\\}\mathscr{N}^{S}_{2}\\{\boldsymbol{q_{1}}\\}+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{1}}\\}~{}~{}~{}~{}~{}~{}\\\
+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{0}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(39)

$\begin{split}\mathbb{N}^{S,\epsilon^{3}}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}},\boldsymbol{q_{2}},\boldsymbol{q_{3}}\\}=\mathscr{N}^{S}_{1}\\{\boldsymbol{q_{3}}\\}+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{q_{3}}\\}\mathscr{N}^{S}_{2}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{S}_{2}\\{\boldsymbol{q_{3}}\\}~{}~{}~{}~{}\\\
+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{3}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{3}}\\}~{}~{}\\\
+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{q_{1}}\\}\mathscr{N}^{S}_{2}\\{\boldsymbol{q_{2}}\\}+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{2}}\\}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{2}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{q_{2}}\\}\mathscr{N}^{S}_{2}\\{\boldsymbol{q_{1}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{2}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{1}}\\}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{2}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{0}}\\}~{}\\\
+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{1}}\\}\mathscr{N}^{S}_{3}\\{\boldsymbol{q_{1}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(40)

$\mathscr{N}^{SS}\\{\boldsymbol{\tilde{q}}\\}$ is also expanded in similar
way:

$\mathscr{N}^{SS}\\{\boldsymbol{\tilde{q}}\\}=\mathbb{N}^{SS,1}\\{\boldsymbol{q_{0}}\\}+\epsilon\mathbb{N}^{SS,\epsilon}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}}\\}+\epsilon^{2}\mathbb{N}^{SS,\epsilon^{2}}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}},\boldsymbol{q_{2}}\\}+\epsilon^{3}\mathbb{N}^{SS,\epsilon^{3}}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}},\boldsymbol{q_{2}},\boldsymbol{q_{3}}\\}+...$
(41)

where,

$\begin{split}\mathbb{N}^{SS,1}\\{\boldsymbol{q_{0}}\\}=\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{SS}_{2}\\{\boldsymbol{q_{0}}\\}\end{split}$
(42)

$\begin{split}\mathbb{N}^{SS,\epsilon}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}}\\}=\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{SS}_{2}\\{\boldsymbol{q_{1}}\\}+\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{q_{1}}\\}\mathscr{N}^{SS}_{2}\\{\boldsymbol{q_{0}}\\}\end{split}$
(43)

$\begin{split}\mathbb{N}^{SS,\epsilon^{2}}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}},\boldsymbol{q_{2}}\\}=\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{SS}_{2}\\{\boldsymbol{q_{2}}\\}+\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{q_{2}}\\}\mathscr{N}^{SS}_{2}\\{\boldsymbol{q_{0}}\\}+\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{q_{1}}\\}\mathscr{N}^{SS}_{2}\\{\boldsymbol{q_{1}}\\}\end{split}$
(44)

$\begin{split}\mathbb{N}^{SS,\epsilon^{3}}\\{\boldsymbol{q_{0}},\boldsymbol{q_{1}},\boldsymbol{q_{2}},\boldsymbol{q_{3}}\\}=\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{SS}_{2}\\{\boldsymbol{q_{3}}\\}+\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{q_{3}}\\}\mathscr{N}^{SS}_{2}\\{\boldsymbol{q_{0}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{q_{1}}\\}\mathscr{N}^{SS}_{2}\\{\boldsymbol{q_{2}}\\}+\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{q_{2}}\\}\mathscr{N}^{SS}_{2}\\{\boldsymbol{q_{1}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(45)

We define the linearized Navier-Stokes operator i.e. $\mathscr{L}$ acting on
any vector field $\boldsymbol{p}$ as follows:

$\begin{split}\mathscr{L}\boldsymbol{p}=\big{(}\mathbb{N}^{1}\\{\boldsymbol{q_{0}}\\}\boldsymbol{p}+\mathbb{N}^{\epsilon}\\{\boldsymbol{q_{0}},\boldsymbol{p}\\}\boldsymbol{q_{0}}\big{)}+S_{c}\big{(}\mathbb{N}^{S,1}\\{\boldsymbol{q_{0}}\\}\boldsymbol{p}+\mathbb{N}^{S,\epsilon}\\{\boldsymbol{q_{0}},\boldsymbol{p}\\}\boldsymbol{q_{0}}\big{)}\\\
+~{}S_{c}^{2}\big{(}\mathbb{N}^{SS,1}\\{\boldsymbol{q_{0}}\\}\boldsymbol{p}+\mathbb{N}^{SS,\epsilon}\\{\boldsymbol{q_{0}},\boldsymbol{p}\\}\boldsymbol{q_{0}}\big{)}-\big{(}\mathcal{L}_{T}\boldsymbol{p}+S_{c}\mathcal{L}_{T}^{S}\boldsymbol{p}\big{)}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(46)

We define the vector field
$\boldsymbol{p}(r,\theta,z,t_{1},t_{2})=A_{p}(t_{2})\boldsymbol{\hat{p}_{m}}(r,z)e^{i(m\theta-\omega
t_{1})}$. The operator $\mathscr{L}_{m}$ is obtained from $\mathscr{L}$ by the
transformation $\partial/\partial\theta\to im$ and the operation
$\mathscr{L}_{m}\boldsymbol{\hat{p}_{m}}$ can be written as follows:

$\begin{split}\mathscr{L}_{m}\boldsymbol{\hat{p}_{m}}=\big{(}\mathbb{N}^{1}_{(0,m)}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{p}_{m}}+\mathbb{N}^{\epsilon}_{(m,0)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}_{m}}\\}\boldsymbol{q_{0}}\big{)}+S_{c}\big{(}\mathbb{N}^{S,1}_{m}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{p}_{m}}+\mathbb{N}^{S,\epsilon}_{m}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}_{m}}\\}\boldsymbol{q_{0}}\big{)}\\\
+~{}S_{c}^{2}\big{(}\mathbb{N}^{SS,1}_{m}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{p}_{m}}+\mathbb{N}^{SS,\epsilon}_{m}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}_{m}}\\}\boldsymbol{q_{0}}\big{)}-\big{(}\mathcal{L}_{T,m}\boldsymbol{\hat{p}_{m}}+S_{c}\mathcal{L}_{T,m}^{S}\boldsymbol{\hat{p}_{m}}\big{)}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(47)

We now define following vector fields to explain the various operator-vector
operations used in eq. 47:

$\begin{split}\boldsymbol{r}(r,\theta,z,t_{1},t_{2})=A_{r}(t_{2})\boldsymbol{\hat{r}}(r,z)e^{i(r\theta-\omega_{r}t_{1})}\\\
\boldsymbol{s}(r,\theta,z,t_{1},t_{2})=A_{s}(t_{2})\boldsymbol{\hat{s}}(r,z)e^{i(s\theta-\omega_{s}t_{1})}\end{split}$
(48)

$\mathbb{N}^{1}\\{\boldsymbol{q_{0}}\\}\boldsymbol{s}$ operation can be
expanded as follows:

$\mathbb{N}^{1}\\{\boldsymbol{q_{0}}\\}\boldsymbol{s}=A_{s}e^{is\theta}e^{-i\omega_{s}t_{1}}\mathbb{N}^{1}_{(0,s)}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}$
(49)

Using eq. 32,

$\mathbb{N}^{1}_{(0,s)}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}=\mathscr{N}_{1}^{(0,s)}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{q_{0}}\\}\mathscr{N}_{2}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}$
(50)

Similarly,

$\mathbb{N}^{\epsilon}\\{\boldsymbol{q_{0}},\boldsymbol{r}\\}\boldsymbol{s}=A_{r}A_{s}e^{i(r+s)\theta}e^{-i(\omega_{r}+\omega_{s})t_{1}}\mathbb{N}^{\epsilon}_{(r,s)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}$
(51)

Using eq. 33,

$\begin{split}\mathbb{N}^{\epsilon}_{(r,s)}\\{\boldsymbol{q_{0}},\boldsymbol{r}\\}\boldsymbol{s}=\mathscr{N}_{1}^{(r,s)}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{q_{0}}\\}\mathscr{N}_{2}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}+\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{\hat{r}}\\}\mathscr{N}_{2}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}\\\
+~{}2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{r}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{N}_{3}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}\end{split}$
(52)

The matrix-vector operation
$\mathscr{N}_{1}^{(r,s)}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}$ is
obtained by setting $\partial\tilde{\rho}/\partial\theta=ir\tilde{\rho}$ in
the first element of last row and $\partial/\partial\theta=is$ in the
expression for $\mathscr{N}_{1}\\{\boldsymbol{\tilde{q}}\\}$ in eq. 23.

Similarly, we expand the various other operations used in eq. 47 using eqs.
37, 38, 42 and 43 as follows:

$\mathbb{N}^{S,1}_{s}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}=\mathscr{N}^{S,s}_{1}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{p}}\\}\mathscr{N}^{S,s}_{2}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{S,s}_{3}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}$
(53)
$\begin{split}\mathbb{N}^{S,\epsilon}_{s}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}=\mathscr{N}^{S,s}_{1}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{\hat{r}}\\}\mathscr{N}^{S,s}_{2}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}+\mathscr{R}^{\rho}_{2,3,4}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{S,s}_{2}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}\\\
+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{r}}\\}\mathscr{N}^{S,s}_{3}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}+\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{S,s}_{3}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}\end{split}$
(54)

$\mathbb{N}^{SS,1}_{s}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}=\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{SS,s}_{2}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}$
(55)
$\begin{split}\mathbb{N}^{SS,\epsilon}_{s}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}=\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{q_{0}}\\}\mathscr{N}^{SS,s}_{2}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}+\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{\hat{r}}\\}\mathscr{N}^{SS,s}_{2}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}\end{split}$
(56)

##  Appendix B Coefficients in the evolution equations of $A_{1}$ and $A_{0}$

The coefficients in the evolution equation for $A_{1}(t_{2})$ i.e. eq. 11 are
defined as follows:

$\begin{split}\alpha_{A_{1}}=\frac{1}{\langle\boldsymbol{\hat{q}_{1}^{\dagger}},\mathscr{B}_{0}\boldsymbol{\hat{q}_{1}}\rangle}\bigg{\langle}\boldsymbol{\hat{q}_{1}^{\dagger}},-\Big{[}(-i\omega_{1})\mathscr{B}_{1}^{S}\\{\boldsymbol{\hat{q}_{\Delta}}\\}\boldsymbol{\hat{q}_{1}}+\mathbb{S}_{(0,1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{\Delta}}\\}\boldsymbol{\hat{q}_{1}}+\mathbb{S}_{(1,0)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{\hat{q}_{\Delta}}\\\
+~{}\mathbb{R}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{\Delta}}\\}\boldsymbol{q_{0}}+\mathbb{R}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{\Delta}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{q_{0}}+\mathbb{N}_{1}^{S,1}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{q}_{1}}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}2S_{c}\mathbb{N}_{1}^{SS,1}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{q}_{1}}-\mathcal{L}_{T,0}^{S}\boldsymbol{\hat{q}_{1}}+(-i\omega_{1})\mathscr{B}^{S}_{2}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{q}_{1}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}\mathbb{N}^{S,\epsilon}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{q_{0}}+2S_{c}\mathbb{N}^{SS,\epsilon}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{q_{0}}\Big{]}\bigg{\rangle}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(57)
$\begin{split}\beta_{A_{1}A_{1}}=\frac{1}{\langle\boldsymbol{\hat{q}_{1}^{\dagger}},\mathscr{B}_{0}\boldsymbol{\hat{q}_{1}}\rangle}\bigg{\langle}\boldsymbol{\hat{q}_{1}^{\dagger}},\Big{[}\mathbb{I}_{(-1,2)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}^{*}}\\}\boldsymbol{\hat{q}_{A_{1}A_{1}}}+\mathbb{I}_{(0,1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{A_{1}A_{1}^{*}}}\\}\boldsymbol{\hat{q}_{1}}+\mathbb{I}_{(2,-1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{A_{1}A_{1}}}\\}\boldsymbol{\hat{q}_{1}}\\\
+~{}\mathbb{S}_{(1,0)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{\hat{q}_{A_{1}A_{1}^{*}}}+\mathbb{R}_{-1}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{\hat{q}_{1}^{*}}+\mathbb{Q}_{1}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{\hat{q}_{1}}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}\mathbb{Q}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{A_{1}A_{1}^{*}}}\\}\boldsymbol{q_{0}}+\mathbb{Q}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{A_{1}A_{1}}}\\}\boldsymbol{q_{0}}+\mathbb{P}^{\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}^{*}}\\}\boldsymbol{q_{0}}~{}~{}~{}~{}\\\
+~{}\mathbb{P}^{\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{q_{0}}+\mathbb{P}^{\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{q_{0}}+S_{c}\mathbb{P}_{0}^{S,\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}^{*}}\\}\boldsymbol{q_{0}}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}S_{c}\mathbb{P}_{0}^{S,\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{q_{0}}+S_{c}\mathbb{P}_{0}^{S,\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{q_{0}}\Big{]}\bigg{\rangle}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(58)

The coefficients in the evolution equation for $A_{0}(t_{2})$ i.e. eq. 10 are
defined as follows:

$\begin{split}\alpha_{A_{0}}=\frac{1}{\langle\boldsymbol{\hat{q}_{0}^{\dagger}},\mathscr{B}_{0}\boldsymbol{\hat{q}_{0}}\rangle}\bigg{\langle}\boldsymbol{\hat{q}_{0}^{\dagger}},-\Big{[}(-i\omega_{0})\mathscr{B}_{1}^{S}\\{\boldsymbol{\hat{q}_{\Delta}}\\}\boldsymbol{\hat{q}_{0}}+\mathbb{S}_{(0,0)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{\Delta}}\\}\boldsymbol{\hat{q}_{0}}+\mathbb{S}_{(0,0)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}}\\}\boldsymbol{\hat{q}_{\Delta}}\\\
+~{}\mathbb{R}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{\Delta}}\\}\boldsymbol{q_{0}}+\mathbb{R}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{\Delta}},\boldsymbol{\hat{q}_{0}},\\}\boldsymbol{q_{0}}+\mathbb{N}_{1}^{S,1}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{q}_{0}}~{}~{}~{}~{}~{}~{}\\\
+~{}2S_{c}\mathbb{N}_{1}^{SS,1}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{q}_{0}}-\mathcal{L}_{T,0}^{S}\boldsymbol{\hat{q}^{0}}+(-i\omega_{0})\mathscr{B}^{S}_{2}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{q}_{0}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}\mathbb{N}^{S,\epsilon}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}}\\}\boldsymbol{q_{0}}+2S_{c}\mathbb{N}^{SS,\epsilon}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}}\\}\boldsymbol{q_{0}}\Big{]}\bigg{\rangle}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(59)
$\beta_{A_{0}f}=\frac{\langle\boldsymbol{\hat{q}_{0}^{\dagger}},\boldsymbol{\hat{q}_{a}}\rangle}{2\langle\boldsymbol{\hat{q}_{0}^{\dagger}},\mathscr{B}_{0}\boldsymbol{\hat{q}_{0}}\rangle}$
(60)
$\begin{split}\beta_{A_{0}A_{1}}=\frac{1}{\langle\boldsymbol{\hat{q}_{0}^{\dagger}},\mathscr{B}_{0}\boldsymbol{\hat{q}_{0}}\rangle}\bigg{\langle}\boldsymbol{\hat{q}^{0\dagger}},\Big{[}\mathbb{I}_{(-1,1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}^{*}}\\}\boldsymbol{\hat{q}_{A_{1}A_{0}}}+\mathbb{I}_{(1,-1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{A_{1}A_{0}}}\\}\boldsymbol{\hat{q}_{1}^{*}}+\mathbb{I}_{(-1,1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}\\}\boldsymbol{\hat{q}_{1}}\\\
+~{}\mathbb{I}_{(1,-1)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}+\mathbb{I}_{(0,0)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{A_{1}A_{1}^{*}}}\\}\boldsymbol{\hat{q}_{0}}+\mathbb{S}_{(0,0)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}}\\}\boldsymbol{\hat{q}_{A_{1}A_{1}^{*}}}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}\mathbb{Q}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{1}}\\}\boldsymbol{\hat{q}_{0}}+\mathbb{Q}_{1}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{0}}\\}\boldsymbol{\hat{q}_{1}}+\mathbb{Q}_{-1}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{0}}\\}\boldsymbol{\hat{q}_{1}^{*}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+\mathbb{T}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{A_{1}A_{1}^{*}}},\boldsymbol{\hat{q}_{A_{1}A_{0}}},\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}\\}\boldsymbol{q_{0}}\Big{]}\bigg{\rangle}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(61)

##  Appendix C Operator definitions

We define following vector fields for describing the various matrix operators
used in this paper:

$\begin{split}\boldsymbol{\tilde{p}}(r,\theta,z,t_{1},t_{2})=A_{p}(t_{2})\boldsymbol{\hat{p}}(r,z)e^{i(p\theta-\omega_{p}t_{1})}\\\
\boldsymbol{\tilde{r}}(r,\theta,z,t_{1},t_{2})=A_{r}(t_{2})\boldsymbol{\hat{r}}(r,z)e^{i(r\theta-\omega_{r}t_{1})}~{}\\\
\boldsymbol{\tilde{s}}(r,\theta,z,t_{1},t_{2})=A_{s}(t_{2})\boldsymbol{\hat{s}}(r,z)e^{i(s\theta-\omega_{s}t_{1})}~{}\end{split}$
(62)
$\begin{split}\mathbb{T}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{A_{1}A_{1}^{*}}},\boldsymbol{\hat{q}_{A_{1}A_{0}}},\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}\\}\boldsymbol{q_{0}}=\Big{[}\mathbb{Q}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{A_{1}A_{1}^{*}}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}\mathbb{Q}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{A_{1}^{*}A_{0}}}\\}+\mathbb{Q}_{0}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{A_{1}A_{0}}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}\mathbb{P}^{\epsilon^{3}}\\{\boldsymbol{\hat{q}^{0}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}^{*}}\\}+\mathbb{P}^{\epsilon^{3}}\\{\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{1}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}\mathbb{P}^{\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{1}^{*}}\\}+\mathbb{P}^{\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{1}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}\mathbb{P}^{\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{0}}\\}+\mathbb{P}^{\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{0}}\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}S_{c}\big{(}\mathbb{P}_{0}^{S,\epsilon^{3}}\\{\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}^{*}}\\}+\mathbb{P}_{0}^{S,\epsilon^{3}}\\{\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{1}}\\}\big{)}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}S_{c}\big{(}\mathbb{P}_{0}^{S,\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{1}^{*}}\\}+\mathbb{P}_{0}^{S,\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{0}},\boldsymbol{\hat{q}_{1}}\\}\big{)}~{}~{}~{}~{}~{}~{}~{}\\\
+~{}S_{c}\big{(}\mathbb{P}_{0}^{S,\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}^{0}}\\}+\mathbb{P}_{0}^{S,\epsilon^{3}}\\{\boldsymbol{\hat{q}_{1}^{*}},\boldsymbol{\hat{q}_{1}},\boldsymbol{\hat{q}_{0}}\\}\big{)}\Big{]}\boldsymbol{q_{0}}~{}~{}~{}\end{split}$
(63)

where,

$\begin{split}\mathbb{P}^{\epsilon^{3}}\\{\boldsymbol{\hat{p}},\boldsymbol{\hat{r}},\boldsymbol{\hat{s}}\\}\boldsymbol{q_{0}}=\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{p}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{r}}\\}\mathscr{N}_{3}\\{\boldsymbol{\hat{s}}\\}\boldsymbol{q_{0}}\end{split}$
(64)
$\begin{split}\mathbb{P}^{S,\epsilon^{3}}_{0}\\{\boldsymbol{\hat{p}},\boldsymbol{\hat{r}},\boldsymbol{\hat{s}}\\}\boldsymbol{q_{0}}=\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{p}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{r}}\\}\mathscr{N}^{S,0}_{3}\\{\boldsymbol{\hat{s}}\\}\boldsymbol{q_{0}}\end{split}$
(65)
$\begin{split}\mathbb{Q}_{s}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}=\mathbb{R}_{s}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}+\mathbb{R}_{s}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{r}},\boldsymbol{\hat{p}}\\}\boldsymbol{\hat{s}}\end{split}$
(66)
$\begin{split}\mathbb{I}_{(p,r)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}}\\}\boldsymbol{\hat{r}}=(-i\omega_{r})\mathscr{B}^{S}_{1}\\{\boldsymbol{\hat{p}}\\}\boldsymbol{\hat{r}}+\mathbb{S}_{(p,r)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}}\\}\boldsymbol{\hat{r}}\end{split}$
(67)
$\begin{split}\mathbb{R}_{s}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}=\mathbb{P}^{\epsilon^{2}}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}+S_{c}\mathbb{P}^{S,\epsilon^{2}}_{s}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}+S^{2}_{c}\mathbb{P}^{SS,\epsilon^{2}}_{s}\\{\boldsymbol{\hat{p}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}\end{split}$
(68)
$\begin{split}\mathbb{S}_{(p,r)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}}\\}\boldsymbol{\hat{r}}=\mathbb{N}^{\epsilon}_{(p,r)}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}}\\}\boldsymbol{\hat{r}}+S_{c}\mathbb{N}^{S,\epsilon}_{r}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{r}}+S^{2}_{c}\mathbb{N}^{SS,\epsilon}_{r}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}}\\}\boldsymbol{\hat{r}}\end{split}$
(69)
$\begin{split}\mathbb{P}^{\epsilon^{2}}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}=\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{\hat{p}}\\}\mathscr{N}_{2}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{p}}\\}\mathscr{N}_{3}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}\\\
+~{}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{p}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{r}}\\}\mathscr{N}_{3}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(70)
$\begin{split}\mathbb{P}^{S,\epsilon^{2}}_{s}\\{\boldsymbol{q_{0}},\boldsymbol{\hat{p}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}=\mathscr{R}^{\rho}_{2,4}\\{\boldsymbol{\hat{p}}\\}\mathscr{N}^{S,s}_{2}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}+2\mathscr{R}^{\rho}_{5}\\{\boldsymbol{q_{0}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{p}}\\}\mathscr{N}^{S,s}_{3}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}\\\
+~{}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{p}}\\}\mathscr{R}^{\rho}_{5}\\{\boldsymbol{\hat{r}}\\}\mathscr{N}^{S,s}_{3}\\{\boldsymbol{q_{0}}\\}\boldsymbol{\hat{s}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\end{split}$
(71)
$\begin{split}\mathbb{P}^{SS,\epsilon^{2}}_{s}\\{\boldsymbol{\hat{p}},\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}=\mathscr{R}^{\rho}_{2,3}\\{\boldsymbol{\hat{p}}\\}\mathscr{N}^{SS,s}_{2}\\{\boldsymbol{\hat{r}}\\}\boldsymbol{\hat{s}}\end{split}$
(72)

## References

  * (1)
  * Acharya et al. (2012) Acharya, V., Shreekrishna, Shin, D.-H. & Lieuwen, T. (2012), ‘Swirl effects on harmonically excited, premixed flame kinematics’, Combustion and Flame 159(3), 1139 – 1150.
http://www.sciencedirect.com/science/article/pii/S0010218011003002

  * Anacleto et al. (2003) Anacleto, P. M., Fernandes, E. C., Heitor, M. V. & Shtork, S. I. (2003), ‘Swirl flow structure and flame characteristics in a model lean premixed combustor’, Combustion Science and Technology 175(8), 1369–1388.
  * Auer et al. (2006) Auer, M., Hirsch, C. & Sattelmayer, T. (2006), Influence of air and fuel mass flow fluctuations in a premix swirl burner on flame dynamics, in ‘ASME Turbo Expo 2006: Power for Land, Sea, and Air’, American Society of Mechanical Engineers Digital Collection.
  * Batchelor & Gill (1962) Batchelor, G. K. & Gill, A. E. (1962), ‘Analysis of the stability of axisymmetric jets’, Journal of Fluid Mechanics 14(4), 529–551.
  * Berkooz et al. (1993) Berkooz, G., Holmes, P. & Lumley, J. L. (1993), ‘The proper orthogonal decomposition in the analysis of turbulent flows’, Annual Review of Fluid Mechanics 25(1), 539–575.
  * Bluemner et al. (2019) Bluemner, R., Paschereit, C. O. & Oberleithner, K. (2019), ‘Generation and transport of equivalence ratio fluctuations in an acoustically forced swirl burner’, Combustion and Flame 209, 99–116.
  * Candel (2002) Candel, S. (2002), ‘Combustion dynamics and control: Progress and challenges’, Proceedings of the Combustion Institute 29(1), 1–28.
  * Candel et al. (2014) Candel, S., Durox, D., Schuller, T., Bourgouin, J.-F. & Moeck, J. P. (2014), ‘Dynamics of swirling flames’, Annual Review of Fluid Mechanics 46, 147–173.
  * Correa (1998) Correa, S. M. (1998), Power generation and aeropropulsion gas turbines: From combustion science to combustion technology, in ‘Symposium (International) on Combustion’, Vol. 27, pp. 1793–1807.
  * Emerson et al. (2012) Emerson, B., O’Connor, J., Juniper, M. & Lieuwen, T. (2012), ‘Density ratio effects on reacting bluff-body flow field characteristics’, Journal of Fluid Mechanics 706, 219–250.
  * Escudier & Keller (1985) Escudier, M. P. & Keller, J. J. (1985), ‘Recirculation in swirling flow - a manifestation of vortex breakdown’, AIAA Journal 23(1), 111–116.
  * Frederick et al. (2018) Frederick, M., Manoharan, K., Dudash, J., Brubaker, B., Hemchandra, S. & O’Connor, J. (2018), ‘Impact of precessing vortex core dynamics on shear layer response in a swirling jet’, Journal of Engineering for Gas Turbines and Power 140(6), 061503.
  * Freitag et al. (2005) Freitag, M., Klein, M., Gregor, M., Nauert, A., Geyer, D., Schneider, C., Dreizler, A. & Janicka, J. (2005), Mixing analysis of a swirling recirculating flow using dns and experimental data, in ‘Fourth International Symposium on Turbulence and Shear Flow Phenomena’.
  * Geigle et al. (2017) Geigle, K. P., Hadef, R., Stöhr, M. & Meier, W. (2017), ‘Flow field characterization of pressurized sooting swirl flames and relation to soot distributions’, Proceedings of the Combustion Institute 36(3), 3917–3924.
  * Geigle et al. (2015) Geigle, K. P., Köhler, M., O’Loughlin, W. & Meier, W. (2015), ‘Investigation of soot formation in pressurized swirl flames by laser measurements of temperature, flame structures and soot concentrations’, Proceedings of the Combustion Institute 35(3), 3373–3380.
  * Ghoniem & Givi (1988) Ghoniem, A. F. & Givi, P. (1988), ‘Lagrangian simulation of a reacting mixing layer at low heat release’, AIAA Journal 26(6), 690–697.
  * Gonzalez et al. (2005) Gonzalez, E., Lee, J. & Santavicca, D. (2005), A study of combustion instabilities driven by flame-vortex interactions, in ‘41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit’.
  * Gupta et al. (1984) Gupta, A. K., Lilley, D. G. & Syred, N. (1984), Swirl flows, Abacus press.
  * Hall (1972) Hall, M. (1972), ‘Vortex breakdown’, Annual Review of Fluid Mechanics 4(1), 195–218.
  * Harvey (1962) Harvey, J. (1962), ‘Some observations of the vortex breakdown phenomenon’, Journal of Fluid Mechanics 14(4), 585–592.
  * Hemchandra et al. (2018) Hemchandra, S., Shanbhogue, S., Hong, S. & Ghoniem, A. F. (2018), ‘Role of hydrodynamic shear layer stability in driving combustion instability in a premixed propane-air backward-facing step combustor’, Physical Review Fluids 3(6), 063201.
  * Joshi et al. (1998) Joshi, N. D., Mongia, H. C., Leonard, G., Stegmaier, J. W. & Vickers, E. C. (1998), Dry low emissions combustor development, in ‘Turbo Expo: Power for Land, Sea, and Air’, Vol. 78644, American Society of Mechanical Engineers, p. V003T06A027.
  * Kang et al. (2007) Kang, D., Culick, F. & Ratner, A. (2007), ‘Combustion dynamics of a low-swirl combustor’, Combustion and Flame 151(3), 412–425.
  * Kim et al. (2010) Kim, K. T., Lee, J. G., Quay, B. D. & Santavicca, D. (2010), ‘Response of partially premixed flames to acoustic velocity and equivalence ratio perturbations’, Combustion and Flame 157(9), 1731–1744.
  * Lee et al. (2000) Lee, J. G., Kim, K. & Santavicca, D. (2000), ‘Measurement of equivalence ratio fluctuation and its effect on heat release during unstable combustion’, Proceedings of the Combustion Institute 28(1), 415–421.
  * Lee et al. (1993) Lee, T.-W., Lee, J., Nye, D. & Santavicca, D. (1993), ‘Local response and surface properties of premixed flames during interactions with kármán vortex streets’, Combustion and Flame 94(1-2), 146–160.
  * Lefebvre (1998) Lefebvre, A. H. (1998), Gas Turbine Combustion, CRC press.
  * Liang & Maxworthy (2005) Liang, H. & Maxworthy, T. (2005), ‘An experimental investigation of swirling jets’, Journal of Fluid Mechanics 525, 115–159.
  * Lieuwen (2012) Lieuwen, T. C. (2012), Unsteady combustor physics, Cambridge University Press.
  * Lieuwen & Yang (2005) Lieuwen, T. C. & Yang, V. (2005), Combustion instabilities in gas turbine engines: operational experience, fundamental mechanisms, and modeling, American Institute of Aeronautics and Astronautics.
  * Lieuwen & Zinn (1998) Lieuwen, T. & Zinn, B. T. (1998), The role of equivalence ratio oscillations in driving combustion instabilities in low nox gas turbines, in ‘Symposium (International) on Combustion’, Vol. 27, pp. 1809–1816.
  * Litvinov et al. (2021) Litvinov, I., Yoon, J., Stöhr, M., Boxx, I. & Geigle, K. (2021), ‘Time-resolved study of mixing and reaction in an aero-engine model combustor at increased pressure’, Combustion and Flame under review.
  * Lückoff & Oberleithner (2019) Lückoff, F. & Oberleithner, K. (2019), ‘Excitation of the precessing vortex core by active flow control to suppress thermoacoustic instabilities in swirl flames’, International Journal of Spray and Combustion Dynamics 11, 1756827719856237.
  * Manoharan et al. (2020) Manoharan, K., Frederick, M., Clees, S., O’Connor, J. & Hemchandra, S. (2020), ‘A weakly nonlinear analysis of the precessing vortex core oscillation in a variable swirl turbulent round jet’, Journal of Fluid Mechanics 884.
  * Manoharan et al. (2015) Manoharan, K., Hansford, S., O’Connor, J. & Hemchandra, S. (2015), Instability mechanism in a swirl flow combustor: precession of vortex core and influence of density gradient, in ‘ASME Turbo Expo 2015: Turbine Technical Conference and Exposition’, American Society of Mechanical Engineers.
  * Manoharan & Hemchandra (2014) Manoharan, K. & Hemchandra, S. (2014), ‘Absolute/Convective Instability Transition in a Backward Facing Step Combustor: Fundamental Mechanism and Influence of Density Gradient’, Journal of Engineering for Gas Turbines and Power 137(2).  021501\.
https://doi.org/10.1115/1.4028206

  * Mathews et al. (2016) Mathews, B., Hansford, S. & O’Connor, J. (2016), Impact of swirling flow structure on shear layer vorticity fluctuation mechanisms, in ‘ASME Turbo Expo 2016: Turbomachinery Technical Conference and Exposition’, American Society of Mechanical Engineers.
  * Moeck et al. (2012) Moeck, J. P., Bourgouin, J.-F., Durox, D., Schuller, T. & Candel, S. (2012), ‘Nonlinear interaction between a precessing vortex core and acoustic oscillations in a turbulent swirling flame’, Combustion and Flame 159(8), 2650–2668.
  * Nayfeh (2008) Nayfeh, A. H. (2008), Perturbation methods, John Wiley & Sons.
  * Oberleithner et al. (2011) Oberleithner, K., Sieber, M., Nayeri, C., Paschereit, C., Petz, C., Hege, H.-C., Noack, B. & Wygnanski, I. (2011), ‘Three-dimensional coherent structures in a swirling jet undergoing vortex breakdown: stability analysis and empirical mode construction’, Journal of Fluid Mechanics 679, 383–414.
  * Oberleithner et al. (2015) Oberleithner, K., Stöhr, M., Im, S. H., Arndt, C. M. & Steinberg, A. M. (2015), ‘Formation and flame-induced suppression of the precessing vortex core in a swirl combustor: Experiments and linear stability analysis’, Combustion and Flame 162(8), 3100–3114.
  * Oberleithner et al. (2013) Oberleithner, K., Terhaar, S., Rukes, L. & Paschereit, C. O. (2013), ‘Why nonuniform density suppresses the precessing vortex core’, Journal of Engineering for Gas Turbines and Power 135(12), 121506.
  * O’Connor & Lieuwen (n.d.) O’Connor, J. & Lieuwen, T. (n.d.), ‘Further characterization of the disturbance field in a transversely excited swirl-stabilized flame’, Journal of Engineering for Gas Turbines and Power 134(1), 011501.
  * Paschereit et al. (1999) Paschereit, C. O., Gutmark, E. & Weisenstein, W. (1999), ‘Coherent structures in swirling flows and their role in acoustic combustion control’, Physics of Fluids 11(9), 2667–2678.
  * Poinsot et al. (1987) Poinsot, T. J., Trouve, A. C., Veynante, D. P., Candel, S. M. & Esposito, E. J. (1987), ‘Vortex-driven acoustically coupled combustion instabilities’, Journal of Fluid Mechanics 177, 265–292.
  * Renard et al. (1999) Renard, P.-H., Rolon, J. C., Thévenin, D. & Candel, S. (1999), ‘Investigations of heat release, extinction, and time evolution of the flame surface, for a nonpremixed flame interacting with a vortex’, Combustion and Flame 117(1-2), 189–205.
  * Renard et al. (2000) Renard, P.-H., Thevenin, D., Rolon, J.-C. & Candel, S. (2000), ‘Dynamics of flame/vortex interactions’, Progress in Energy and Combustion Science 26(3), 225–282.
  * Rigas et al. (2017) Rigas, G., Morgans, A. S. & Morrison, J. F. (2017), ‘Weakly nonlinear modelling of a forced turbulent axisymmetric wake’, Journal of Fluid Mechanics 814, 570–591.
  * Schmid (2010) Schmid, P. J. (2010), ‘Dynamic mode decomposition of numerical and experimental data’, Journal of Fluid Mechanics 656, 5–28.
  * Shanbhogue et al. (2009) Shanbhogue, S., Shin, D.-H., Hemchandra, S., Plaks, D. & Lieuwen, T. (2009), ‘Flame-sheet dynamics of bluff-body stabilized flames during longitudinal acoustic forcing’, Proceedings of the Combustion Institute 32(2), 1787–1794.
  * Shreekrishna et al. (2010) Shreekrishna, Hemchandra, S. & Lieuwen, T. (2010), ‘Premixed flame response to equivalence ratio perturbations’, Combustion Theory and Modelling 14(5), 681–714.
  * Steinberg et al. (2010) Steinberg, A. M., Boxx, I., Stöhr, M., Carter, C. D. & Meier, W. (2010), ‘Flow–flame interactions causing acoustically coupled heat release fluctuations in a thermo-acoustically unstable gas turbine model combustor’, Combustion and Flame 157(12), 2250–2266.
  * Stöhr et al. (2015) Stöhr, M., Arndt, C. M. & Meier, W. (2015), ‘Transient effects of fuel–air mixing in a partially-premixed turbulent swirl flame’, Proceedings of the Combustion Institute 35(3), 3327–3335.
  * Stöhr et al. (2012) Stöhr, M., Boxx, I., Carter, C. D. & Meier, W. (2012), ‘Experimental study of vortex-flame interaction in a gas turbine model combustor’, Combustion and Flame 159(8), 2636–2649.
  * Stöhr et al. (2019) Stöhr, M., Geigle, K., Hadef, R., Boxx, I., Carter, C., Grader, M. & Gerlinger, P. (2019), ‘Time-resolved study of transient soot formation in an aero-engine model combustor at elevated pressure’, Proceedings of the Combustion Institute 37(4), 5421–5428.
  * Stöhr et al. (2017) Stöhr, M., Yin, Z. & Meier, W. (2017), ‘Interaction between velocity fluctuations and equivalence ratio fluctuations during thermoacoustic oscillations in a partially premixed swirl combustor’, Proceedings of the Combustion Institute 36(3), 3907–3915.
  * Syred (2006) Syred, N. (2006), ‘A review of oscillation mechanisms and the role of the precessing vortex core (pvc) in swirl combustion systems’, Progress in Energy and Combustion Science 32(2), 93–161.
  * Taamallah et al. (2016) Taamallah, S., Shanbhogue, S. J. & Ghoniem, A. F. (2016), ‘Turbulent flame stabilization modes in premixed swirl combustion: Physical mechanism and karlovitz number-based criterion’, Combustion and Flame 166, 19–33.
  * Tammisola & Juniper (2016) Tammisola, O. & Juniper, M. P. (2016), ‘Coherent structures in a swirl injector at re = 4800 by nonlinear simulations and linear global modes’, Journal of Fluid Mechanics 792, 620–657.
  * Towne et al. (2018) Towne, A., Schmidt, O. T. & Colonius, T. (2018), ‘Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis’, Journal of Fluid Mechanics 847, 821–867.