## 1. Introduction

Thermoacoustic instability is the result of the mutual coupling between flow dynamics, the unsteady heat release produced by a flame and the surrounding acoustic environment (Dowling & Stow Reference Dowling and Stow2003; Dowling & Morgans Reference Dowling and Morgans2005; Lieuwen & Yang Reference Lieuwen and Yang2005; Culick Reference Culick2006; Poinsot Reference Poinsot2017). Thermoacoustic instability is a problem of major concern for the development of gas turbines that reliably work under a wide range of operating conditions, while producing reduced levels of carbon dioxide and $\textrm {NO}_{x}$ emissions that comply with environmental regulations. During thermoacoustic instability, large-amplitude pressure fluctuations develop inside the combustion chamber and affect the entire engine as undesired vibrations. These vibrations affect the normal operation of the system and reduce the lifespan of the engine. In extreme cases, thermoacoustic instability may induce flashback of the flame, causing severe damage to the system elements (Lieuwen & Yang Reference Lieuwen and Yang2005). Quantitative stability prediction and analysis of thermoacoustic systems require the calculation of complex-valued eigenvalues and their associated eigenvectors. Thermoacoustic eigenvalues can be found by solving a nonlinear eigenvalue problem, which is often derived from the non-homogeneous Helmholtz equation including a feedback term that represents the flame response to acoustic perturbations (e.g. Nicoud *et al.* Reference Nicoud, Benoit, Sensiau and Poinsot2007). This calculation may be computationally demanding if systems with millions of degrees of freedom are considered. In order to calculate the drift of eigenvalues and eigenvectors due to changes in parameters at an affordable computational cost, high-order adjoint-based perturbation theory can instead be used (Mensah, Orchini & Moeck Reference Mensah, Orchini and Moeck2020).

### 1.1. Thermoacoustic eigenvalues: classification and origin

In this study, we consider a finite-dimensional nonlinear eigenvalue problem of the form

where $s$ is the eigenvalue and $\hat {\boldsymbol {p}}$ is the associated eigenvector. Nonlinear eigenproblems appear in various applications in science and engineering beyond thermoacoustics, for example, in vibrations of structures, fluid–structure interactions, nanotechnology (quantum dots), time-delay systems and control theory, to name a few (Friedman & Shinbrot Reference Friedman and Shinbrot1968; Mehrmann & Voss Reference Mehrmann and Voss2004; Betcke *et al.* Reference Betcke, Higham, Mehrmann, Schröder and Tisseur2013; Güttel & Tisseur Reference Güttel and Tisseur2017). The classification of the eigenvalues according to their algebraic and geometric multiplicity, and their thermoacoustic physical origin, is essential, because it reflects certain physical properties of the system, such as symmetry and sensitivity. In the following, we briefly recall the relevant definitions.

An eigenvalue has algebraic multiplicity $a$ if $\partial ^j/\partial s^j\det \left (\boldsymbol{\mathsf{L}}\right )=0$ and $\partial ^a/\partial s^a\det \left (\boldsymbol{\mathsf{L}}\right )\neq 0$, where $j=0,1,\ldots ,a-1$. The geometric multiplicity $g$ of an eigenvalue $s$ is the dimension of the null space of $\boldsymbol{\mathsf{L}}(s)$, i.e. $g\equiv \dim \operatorname {null}\boldsymbol{\mathsf{L}}(s)$. The geometric multiplicity is always less than or equal to the algebraic multiplicity. An eigenvalue is semi-simple if $a=g$, it is defective if $a>g$ and it is called simple if $a=g=1$. Eigenvalues that are not simple are degenerate. Degenerate semi-simple eigenvalues are of relevance in several applications with spatial symmetries, including thermoacoustics. For example, rotationally symmetric annular and can-annular combustors, which are are common in thermoacoustics, feature degenerate semi-simple eigenvalues. An important class of defective eigenvalues are branch-point solutions of the characteristic function, which are known as exceptional points (EPs) (Heiss Reference Heiss2004). As recently shown, these spectral singularities are general features of thermoacoustic systems (Mensah *et al.* Reference Mensah, Magri, Silva, Buschmann and Moeck2018; Orchini *et al.* Reference Orchini, Silva, Mensah and Moeck2020). At EPs, the eigenvalues have infinite sensitivity to infinitesimal perturbations to the system.

From a physical point of view, thermoacoustic eigenvalues can be classified according to the feedback loop between unsteady heat release and acoustics. Unsteady heat release rate generates acoustic waves, which propagate away from the flame until they reach the system's boundaries. After reflection at the boundaries, the acoustic waves impinge on the flame and modulate it, thereby generating new fluctuations in the heat release rate. In this work, we refer to the eigenvalues associated with this feedback loop as of acoustic origin. In recent years, another feedback loop in thermoacoustic systems was discovered: the intrinsic thermoacoustic (ITA) feedback loop (Hoeijmakers *et al.* Reference Hoeijmakers, Kornilov, Arteaga, de Goey and Nijmeijer2014; Bomberg, Emmert & Polifke Reference Bomberg, Emmert and Polifke2015). In the ITA loop, the upstream-travelling acoustic waves produced by the flame directly modulate the upstream velocity (without reflection from the boundaries), which, in turn, causes new fluctuations in the unsteady heat release rate, which closes the loop. The ITA loop is independent of the surrounding acoustic boundaries. It exists in both anechoic environments (Silva *et al.* Reference Silva, Emmert, Jaensch and Polifke2015; Hoeijmakers *et al.* Reference Hoeijmakers, Kornilov, Arteaga, de Goey and Nijmeijer2016) and reflective environments (Emmert *et al.* Reference Emmert, Bomberg, Jaensch and Polifke2017; Mukherjee & Shrira Reference Mukherjee and Shrira2017; Silva *et al.* Reference Silva, Merk, Komarek and Polifke2017; Buschmann, Mensah & Moeck Reference Buschmann, Mensah and Moeck2020*a*; Orchini *et al.* Reference Orchini, Silva, Mensah and Moeck2020). We refer to the eigenvalues associated with this feedback mechanism as of ITA origin.

### 1.2. Adjoint-based methods in thermoacoustics

Thermoacoustic systems may be exceedingly sensitive to small variations in system parameters (Juniper & Sujith Reference Juniper and Sujith2018; Magri Reference Magri2019). For the accurate calculation of these sensitivities, adjoint methods proved to be efficient mathematical and computational tools, as reviewed by Magri (Reference Magri2019). Adjoint methods for thermoacoustic eigenvalue sensitivity analysis were developed for design parameter and passive control by Magri & Juniper (Reference Magri and Juniper2013), and subsequently applied to more complex flames in Magri & Juniper (Reference Magri and Juniper2014); Orchini & Juniper (Reference Orchini and Juniper2016). Rigas *et al.* (Reference Rigas, Jamieson, Li and Juniper2015) tested experimentally adjoint-based predictions, showing that the eigenvalue shift was predicted accurately by adjoint sensitivity analysis. The sensitivity information provided by adjoint methods can be embedded into a gradient-update optimization routine to optimally place and tune acoustic dampers in annular combustors (Mensah & Moeck Reference Mensah and Moeck2017).

The thermoacoustic eigenvalue problem is typically nonlinear in the eigenvalue $s$. Existing methodologies for the solution of nonlinear thermoacoustic eigenproblems utilize iterative schemes (Nicoud *et al.* Reference Nicoud, Benoit, Sensiau and Poinsot2007). This may be expensive, for example, for Helmholtz solvers with tens or hundreds of thousands of degrees of freedom, which makes parametric studies computationally demanding. Adjoint methods can also be exploited to simplify the solution of nonlinear eigenvalue problems. Using them, it is possible to map nonlinear eigenproblems, which are difficult to solve, into a series of linear non-homogeneous equations, which are easier to solve, to approximate eigensolutions to any desired order. For simple eigenvalues, general formulae based on the high-order expansion of the eigenvalue problem have been derived by the thermoacoustic community (Mensah *et al.* Reference Mensah, Orchini and Moeck2020). The same level of generality, however, has not been reached for degenerate thermoacoustic modes, which are often found in practise due to the rotational symmetries of annular and can-annular combustors in gas turbines. In this study, higher-order eigenvalue perturbation expansions of thermoacoustic eigenvalues are extended to the degenerate case, which has challenging mathematical complications, as explained in §§ 2 and 3.

### 1.3. Exceptional points

In the past few years, the theory of EPs has been widely employed to explain physical phenomena, e.g. in quantum mechanics and optics (Heiss Reference Heiss2012; Miri & Alù Reference Miri and Alù2019). In thermoacoustics, recent studies have shown that, in certain areas of the complex-frequency plane, small variations in a parameter of a thermoacoustic system lead to a significant change in the eigenvalues in the complex plane (Silva & Polifke Reference Silva and Polifke2019; Sogaro, Schmid & Morgans Reference Sogaro, Schmid and Morgans2019). These studies, however, did not explain what caused the observed high sensitivities. As highlighted in the present work, high sensitivity is a manifestation of the existence of EPs in the thermoacoustic spectrum. As a matter of fact, the sensitivity to infinitesimal changes in parameters is infinite at EPs (Kato Reference Kato1980). Practically, the existence of EPs is observed via strong veering (due to high sensitivity) of the eigenvalue trajectories in their vicinity.

In recent work, the existence of EPs in the spectrum of a one-dimensional Rijke tube has been shown (Mensah *et al.* Reference Mensah, Magri, Silva, Buschmann and Moeck2018). Orchini *et al.* (Reference Orchini, Silva, Mensah and Moeck2020) extended these findings to realistic configurations, by investigating the thermoacoustic modes associated with the acoustic and ITA loop in three-dimensional longitudinal and annular combustors with an $n$–$\tau$ flame response model. The corresponding thermoacoustic eigenvalues of acoustic and ITA origin were studied in the complex plane for systematic variations of $n$ and $\tau$. It was shown that eigenvalues of acoustic origin can coalesce with eigenvalues of ITA origin, manifesting in EPs. Furthermore, in an annular combustor, an EP may also originate from the coalescence of two eigenvalues of acoustic origin. These eigenvalues were found to be associated with two azimuthal modes, one dominant in the plenum and the other in the combustion chamber. This has analogies with the EPs arising due to the acoustic coupling between a cavity and an acoustic damper. In this respect, Bourquard & Noiray (Reference Bourquard and Noiray2019) experimentally demonstrated the existence of such EPs, and showed that both Helmholtz resonators and quarter-wave tube dampers achieve optimal damping performance when tuned to operate at the EPs of the closed-loop coupled acoustic system. In this study, we relate EPs in the spectrum of thermoacoustic systems with (i) the high sensitivity experienced by some thermoacoustic eigenvalues and (ii) the limits of validity of high-order perturbation methods.

### 1.4. Scope

All the concepts introduced in the introduction – high-order perturbation theory, ITA modes and EPs – have been independently shown to be relevant to thermoacoustic models in recent years, and thoroughly studied. They are, however, strongly interconnected. The objectives of this article are to reveal these interconnections, develop an efficient and accurate method for the calculation of thermoacoustic eigenvalue variations and gain physical insight into the properties and behaviour of the thermoacoustic spectrum.

The article is structured as follows. In § 2 a general theory for high-order adjoint-based perturbation expansion of degenerate semi-simple thermoacoustic eigenvalues and eigenvectors is presented. The role of EPs in relation to high-order perturbation theory is discussed in § 3. It is shown how EPs can be identified numerically, exploiting the high-order perturbation theory presented in § 2. Furthermore, we discuss how knowledge of the location of EPs determines well-defined ranges of convergence for the eigenvalues estimated by perturbation theory. In § 4 we apply the presented high-order perturbation theory to two thermoacoustic cases: a simple eigenvalue of an axial combustor and a degenerate semi-simple eigenvalue of an annular combustor. Lastly, in § 5 we discuss how perturbation theory of degenerate thermoacoustic eigenvalues can also be used at EPs by means of Puiseux series expansions. This highlights the difference between symmetry-induced degenerate modes, which are semi-simple and have finite sensitivity with respect to parameter perturbations, and degenerate modes at EPs, which are defective and have an infinite sensitivity.

## 2. High-order adjoint-based perturbation theory for degenerate thermoacoustic modes

In this section, we present a general formulation of high-order adjoint-based perturbation theory for twofold-degenerate semi-simple eigenvalues. This is the category of eigenvalues under which symmetry-induced degeneracies fall, as, for example, the thermoacoustic eigenvalues of rotationally symmetric annular combustors. With adjoint perturbation theory it is possible to (i) understand if a given perturbation unfolds the degeneracy or not; (ii) track the evolution of the split eigenvalues in the complex plane when a parameter is changed; and (iii) calculate the variation of the split eigenvectors when the parameter is changed. We indicate with

a nonlinear eigenvalue problem that depends on (a set of) parameters $\varepsilon$. Here $\boldsymbol{\mathsf{L}}$ is a linear operator acting on an eigenvector $\hat {\boldsymbol {p}}$. The pairs $(s,\hat {\boldsymbol {p}})$ for which (2.1) is satisfied represent the eigenvalues and eigenvectors of the operator. The operator $\boldsymbol{\mathsf{L}}$ is assumed to have an analytical dependence on the eigenvalue and the parameter(s). No further assumptions are made on the properties of the operator $\boldsymbol{\mathsf{L}}$, which, in general, can be non-self-adjoint or even non-normal. Its corresponding adjoint operator, $\boldsymbol{\mathsf{L}}^{\dagger }$, is defined via

where $\left \langle \left .\cdot \right |\cdot \right \rangle$ is an inner product and $\boldsymbol {f}$ and $\boldsymbol {g}$ are arbitrary complex-valued vectors in their relevant Hilbert spaces. In the following, we adopt the Hermitian inner product $\left \langle \left .\boldsymbol {g}\right |\boldsymbol {f}\right \rangle = \boldsymbol {g}^{H} \boldsymbol {f}$, where the superscript $H$ indicates conjugate transpose. Note that, according to this definition, the direct and adjoint operators have the same eigenvalues (López-Gómez & Mora-Corral Reference López-Gómez and Mora-Corral2007; Güttel & Tisseur Reference Güttel and Tisseur2017) and, in a discretized finite element framework as that used in § 4, the discrete adjoint operator is equivalent to the Hermitian transpose of the direct operator, i.e. $\boldsymbol{\mathsf{L}}^{\dagger }=\boldsymbol{\mathsf{L}}^{H}$. For the thermoacoustic problem, the eigenproblem that we are solving is (Dowling & Stow Reference Dowling and Stow2003; Culick Reference Culick2006; Nicoud *et al.* Reference Nicoud, Benoit, Sensiau and Poinsot2007)

coupled with a set of boundary conditions. The finite-dimensional operator $\boldsymbol{\mathsf{L}}$ and the eigenvector $\hat {\boldsymbol {p}}$ in (2.1) are, respectively, the discretization of the thermoacoustic operator and the acoustic pressure $\hat {p}$ in (2.3). In (2.3), which is valid in the zero mean Mach number limit, $c$ is the speed of sound (which may vary spatially), $\gamma$ is the ratio of specific heats, assumed to be homogeneous, and $\bar {\rho }$, $\bar {Q}$ and $\bar {U}$ are the mean density, heat release rate and flow velocity, respectively. The last term on the left-hand side represents the effect of unsteady heat release on the acoustics, modelled with a so-called $n$–$\tau$ model (Crocco Reference Crocco1965). According to this model, the flame response is proportional to the delayed axial acoustic velocity fluctuations at a reference location upstream of the flame. Together with non-trivial boundary conditions (Nicoud *et al.* Reference Nicoud, Benoit, Sensiau and Poinsot2007), the flame response causes the thermoacoustic operator $\boldsymbol{\mathsf{L}}$ to have non-orthogonal eigenvectors, thus exhibiting a non-normal response. Additionally, the delayed response of the flame also causes the thermoacoustic operator to be nonlinear in the eigenvalue $s$. In the present study, we consider both flame response coefficients, $n$ and $\tau$, as perturbation parameters. We highlight that the use of an $n$–$\tau$ model is not a limitation of the perturbative method that we discuss. The method can be applied to any flame model that is analytic in the eigenvalue. This was demonstrated in Mensah *et al.* (Reference Mensah, Magri, Orchini and Moeck2019) on the basis of an experimentally measured flame response expressed in state-space form.

For semi-simple eigenvalues, the eigenvalue and eigenvector dependence on a parameter $\varepsilon$ is expressed in terms of power series expansions of the form

*a*,

*b*)\begin{equation} s(\varepsilon) \approx s_0 + \sum_{j=1}^N \varepsilon^j s_j, \quad \boldsymbol{\hat p}(\varepsilon) \approx \boldsymbol{\hat p}_0 + \sum_{j=1}^N \varepsilon^j \boldsymbol{\hat p}_j, \end{equation}

where, without loss of generality, the perturbation parameter $\varepsilon$ is centred at a reference value $\varepsilon _0=0$. The coefficients $s_j$ and $\hat {\boldsymbol {p}}_j$ are the $j$th-order corrections to the eigenvalues and eigenvectors, respectively. The approximation symbols in (2.4*a*,*b*) indicate that the power series are truncated at order $N$. In thermoacoustics, arbitrarily high-order adjoint-based perturbation theory for non-degenerate eigenvalues has been presented in Mensah *et al.* (Reference Mensah, Orchini and Moeck2020). We report here the key ideas and results of the method since they serve as a starting point for the discussion of the degenerate case, which is the main focus of this study. It is convenient to define the shorthand

The power series approximations (2.4*a*,*b*) are substituted into the eigenvalue problem (2.1), which is then expanded into a Taylor series. By collecting the terms at every order of $\varepsilon$, one obtains a series of linear, non-homogeneous equations that need to be solved in ascending order:

We refer to Mensah *et al.* (Reference Mensah, Orchini and Moeck2020) for a detailed derivation of (2.6). The complexity of the equations is hidden in the $\boldsymbol {r}_j$ terms, which (i) contain all the possible ways of distributing $j$ derivatives between $s$, $\varepsilon$ and $\hat {\boldsymbol {p}}$ and (ii) are functions of the eigenvalue and eigenvector corrections $s_k$ and $\hat {\boldsymbol {p}}_k$ at orders $k<j$. Explicit expressions for the list of all the terms that compose $\boldsymbol {r}_j$ at any order can be analytically obtained. This helps the recursive implementation of perturbation theory (Mensah *et al.* Reference Mensah, Orchini and Moeck2020). In appendix A, we provide a general formula for $\boldsymbol {r}_j$ at any order, and its explicit expressions for $j=1,2$.

The solution strategy becomes straightforward: at any order, a solvability condition based on the Fredholm alternative is imposed, by projecting the right-hand side of (2.6) onto the adjoint eigenvector $\hat {\boldsymbol {p}}_0^{\dagger }$, defined by $\boldsymbol{\mathsf{L}}_{0,0}^{H}\,\hat {\boldsymbol {p}}_0^{\dagger }=\boldsymbol {0}$. This yields a general equation for the eigenvalue correction at order $j$:

Note that, at first order, for which $\boldsymbol {r}_1 = \boldsymbol{\mathsf{L}}_{0,1}\hat {\boldsymbol {p}}_0$, one finds $s_1 = \left . -\left \langle \left .\hat {\boldsymbol {p}}_0^{\dagger }\right |\boldsymbol{\mathsf{L}}_{0,1}\hat {\boldsymbol {p}}_0\right \rangle \right /\left \langle \left .\hat {\boldsymbol {p}}_0^{\dagger }\right |\boldsymbol{\mathsf{L}}_{1,0}\hat {\boldsymbol {p}}_0\right \rangle$, retrieving the known first-order sensitivity expression for nonlinear eigenvalue problems (Magri, Bauerheim & Juniper Reference Magri, Bauerheim and Juniper2016). Once the eigenvalue correction at order $j$ is known, it can be substituted back into the linear systems (2.6), which can then be solved with standard methods for $\hat {\boldsymbol {p}}_j$. Although its solution is not unique – (2.6) is an underdetermined system of equations – the ambiguity in the solution can always be addressed by choosing a normalization condition for the eigenvectors. With both eigenvalue and eigenvector corrections at order $j$, we can move to order $j+1$ and repeat the procedure, up to any desired order.

High-order expressions of thermoacoustic eigenvalue sensitivities have not been developed for the degenerate case. The state of the art is a second-order analysis for the eigenvalues only (Magri *et al.* Reference Magri, Bauerheim and Juniper2016; Mensah *et al.* Reference Mensah, Magri, Orchini and Moeck2019). Starting from the procedure outlined above, in the following we show how perturbation theory can be generalized to handle degenerate semi-simple eigenvalues. We show that, in order to develop a theory generalizable to arbitrarily high order, perturbation theory of degenerate semi-simple eigenvalues needs to be carried out in parallel on both members of the degenerate eigenvalue pair.

### 2.1. Baseline and adjoint degenerate solution

As in any perturbative method, we first require a baseline solution. We assume that the baseline solution, obtained for $\varepsilon = 0$, is degenerate with algebraic multiplicity 2 and semi-simple, so that the geometric multiplicity is also 2. We thus have an unperturbed eigenvalue $s_0$ with an associated two-dimensional subspace spanned by two eigenvectors, denoted $\hat {\tilde {\boldsymbol {p}}}_{0,1}$ and $\hat {\tilde {\boldsymbol {p}}}_{0,2}$, which are chosen to be orthonormal without loss of generality. The first subscript refers to the expansion order and the second distinguishes between the modes in the degenerate eigenvalue. The tilde symbol highlights that the choice of these vectors is not unique. We also need to calculate the associated adjoint eigenvectors, $\hat {\tilde {\boldsymbol {p}}}^{\dagger }_{0,1}$ and $\hat {\tilde {\boldsymbol {p}}}^{\dagger }_{0,2}$, which satisfy $\boldsymbol{\mathsf{L}}_{0,0}^{H}\,\hat {\tilde {\boldsymbol {p}}}^{\dagger }_{0,\zeta }=\boldsymbol {0}$, for $\zeta = 1,2$. As a convention, the subscripts of the following equations contain Latin letters to indicate the perturbation order and Greek letters to distinguish between the (two) degenerate modes.

For semi-simple eigenvalues, it can be shown that the direct and adjoint eigenvectors can always be chosen to satisfy the bi-orthonormalization condition (Güttel & Tisseur Reference Güttel and Tisseur2017)

with $\boldsymbol{\mathsf{L}}_{1,0}$ defined via (2.5). This condition is valid also for non-normal operators, and we will adopt it to simplify the perturbative equations.

### 2.2. Solvability conditions

Because the operator $\boldsymbol{\mathsf{L}}_{0,0}$ has a nullspace of dimension two (spanned by $\hat {\tilde {\boldsymbol {p}}}_{0,\zeta }$), each of the perturbative (2.6) requires two solvability conditions. More specifically, for the equations to admit solutions, their right-hand side must be orthogonal to the (two-dimensional) adjoint subspace spanned by $\hat {\tilde {\boldsymbol {p}}}^{\dagger }_{0,\zeta }$. Depending on whether eigenvalue splitting has occurred or not, different solution strategies need to be employed. This is discussed in the following sections.

#### 2.2.1. Case 1: degeneracy is not resolved

As long as the perturbation considered does not resolve the degeneracy (e.g. for perturbations that preserve the symmetry of the problem), the perturbed eigenvalues will remain degenerate, and the ambiguity in the choice of a basis in the nullspace of the perturbed operator will persist. Thus, at an arbitrary order $j$, we have that the two eigenvalue corrections at orders $0\leq k < j$ are identical, $s_{k,1} = s_{k,2} = s_k$, and the degenerate subspace correction is given by the linear combination

where the $\alpha _{\zeta }$ coefficients are, without loss of generality, chosen to be identical at every order $k$. We are therefore still tracking one degenerate eigenvalue, governed by (2.6). By imposing the two solvability conditions at order $j$, we obtain

*a*)\begin{gather} \left\langle\left.\hat{\tilde{\boldsymbol{p}}}_{0,1}^{\dagger}\right|\boldsymbol{r}_j\right\rangle + \left\langle\left.\hat{\tilde{\boldsymbol{p}}}_{0,1}^{\dagger}\right|s_j\boldsymbol{\mathsf{L}}_{1,0}\hat{\boldsymbol{p}}_0\right\rangle = 0, \end{gather}

*b*)\begin{gather} \left\langle\left.\hat{\tilde{\boldsymbol{p}}}_{0,2}^{\dagger}\right|\boldsymbol{r}_j\right\rangle + \left\langle\left.\hat{\tilde{\boldsymbol{p}}}_{0,2}^{\dagger}\right|s_j\boldsymbol{\mathsf{L}}_{1,0}\hat{\boldsymbol{p}}_0\right\rangle = 0. \end{gather}

Each of the terms contained in $\boldsymbol {r}_j$ is proportional to $\hat {\boldsymbol {p}}_k$ for some $k<j$, which can be expressed as (2.9). By indicating with $\tilde {\boldsymbol {r}}_{j,\zeta }$ the terms of $\boldsymbol {r}_j$ proportional to $\hat {\tilde {\boldsymbol {p}}}_{k,\zeta }$, the solvability conditions (2.10) can be written in matrix form as

where $\boldsymbol {X}_{j_{\eta ,\zeta }} \equiv - \left \langle \left .\hat {\tilde {\boldsymbol {p}}}_{0,\eta }^{\dagger }\right |\tilde {\boldsymbol {r}}_{j,\zeta }\right \rangle$ and $\boldsymbol {\alpha } \equiv [\alpha _1,\alpha _2]^{\textrm {T}}$. Equation (2.11) is a $2\times 2$ linear eigenvalue problem, which we refer to as the auxiliary eigenvalue problem. We need to distinguish two solution cases:

(i) If the two eigenvalues of (2.11) are identical, the problem remains degenerate at this order. We therefore cannot uniquely determine a basis for the eigenvector corrections, but it is convenient to choose them as the solutions of the linear equations

(2.12)\begin{equation} \boldsymbol{\mathsf{L}}_{0,0}\hat{\tilde{\boldsymbol{p}}}_{j,\zeta} = -\tilde{\boldsymbol{r}}_{j,\zeta} - s_j\boldsymbol{\mathsf{L}}_{1,0} \hat{\tilde{\boldsymbol{p}}}_{0,\zeta}, \quad \text{for } \zeta=1,2, \end{equation}so that (2.9) holds also at order $j$, and, at the next order, the same procedure outlined in this subsection can be applied.(ii) If the two eigenvalues of (2.11) are different, the degeneracy unfolds at this order. Together with the eigenvalue corrections $s_{j,\zeta }$, which have different values, we obtain the eigenvectors $\boldsymbol {\alpha }_{\zeta }$ that uniquely determine the directions along which the degenerate subspace of the problem unfolds at lower orders as

(2.13)\begin{equation} \hat{\boldsymbol{p}}_{k,\zeta} = \left[\hat{\tilde{\boldsymbol{p}}}_{k,1}, \hat{\tilde{\boldsymbol{p}}}_{k,2}\right] \boldsymbol{\cdot} \boldsymbol{\alpha}_{\zeta}, \quad \text{for } k = 0,\ldots,j-1. \end{equation}This is the appropriate basis with which to investigate the problem at higher orders because, from (2.4*a*,*b*), it ensures that $\hat {\boldsymbol {p}}_{\zeta }(\varepsilon )$ smoothly approaches $\hat {\boldsymbol {p}}_{0,\zeta }$ when $\varepsilon \rightarrow 0$. To each eigenvalue at order $j$ corresponds an eigenvector correction $\hat {\boldsymbol {p}}_{j,\zeta }$ defined by(2.14)\begin{equation} \boldsymbol{\mathsf{L}}_{0,0}\hat{\boldsymbol{p}}_{j,\zeta} = -\boldsymbol{r}_{j,\zeta} - s_{j,\zeta}\boldsymbol{\mathsf{L}}_{1,0} \hat{\boldsymbol{p}}_{0,\zeta}, \quad \text{for } \zeta=1,2. \end{equation}Note the differences between (2.12) and (2.14): in the latter the tilde symbols have been dropped because the basis is uniquely determined, and an additional index has been appended to the eigenvalue correction at order $j$, as the two eigenvalues now follow different trajectories.

Importantly, the system of equations for the eigenvector corrections, (2.12) or (2.14), admits solutions but is underdetermined since the matrix $\boldsymbol{\mathsf{L}}_{0,0}$ has a non-zero nullspace. Therefore, it admits an infinite number of solutions, which can be expressed as

where $\hat {\boldsymbol {p}}^{\bot }_{j,\zeta }$ is orthogonal to the unperturbed degenerate subspace and $c_{j,\zeta ,\eta }$ are undetermined coefficients (there are two coefficients ($\eta$) for each order ($\,j$) for each eigenvalue ($\zeta$)). As for the non-degenerate case, one coefficient associated with each eigenvector can be determined by imposing a normalization condition on the eigenvectors. The other coefficients, however, are uniquely determined by solvability conditions at higher orders if the eigenvalues split; this is discussed in § 2.2.2.

#### 2.2.2. Case 2: degeneracy is resolved

If at a certain order the perturbation resolves the degeneracy, the eigenvalues split, and we can identify the unique eigendirections along which this splitting occurs. Let us assume that the degeneracy is resolved at order $d$. At orders $n>d$, we are therefore tracking two branches (solutions), whose equations are governed by (2.14), and have as unknowns two eigenvalue and two eigenvector corrections. The four solvability conditions (two for each branch) in this case read

By exploiting the bi-orthonormality condition (2.8), these reduce to

*a*)\begin{gather} s_{j,\zeta} = -\left\langle\left.\hat{\boldsymbol{p}}_{0,\zeta}^{\dagger}\right|\boldsymbol{r}_{j,\zeta}\right\rangle \quad \text{if } \eta=\zeta, \end{gather}

*b*)\begin{gather} \left\langle\left.\hat{\boldsymbol{p}}_{0,\eta}^{\dagger}\right|\boldsymbol{r}_{j,\zeta}\right\rangle = 0 \quad \text{if } \eta \neq \zeta. \end{gather}

The solvability condition (2.17*a*) defines the eigenvalue corrections on each branch $\zeta$ at order $j$, and is identical to the non-degenerate (2.7) when the bi-orthonormalization condition (2.8) is considered. The second condition, (2.17*b*) instead, is new, and belongs to the degenerate case only. It has not been considered by the thermoacoustic community so far, which is why the current state of the art on perturbation theory (Magri *et al.* Reference Magri, Bauerheim and Juniper2016) is limited to second order. If not considered, the solvability conditions are not satisfied, which would then lead to incorrect results in the evaluation of the eigenvector corrections and higher-order coefficients. This fact was first mentioned by Mensah *et al.* (Reference Mensah, Magri, Orchini and Moeck2019) and is formally demonstrated in a complete form in the current study.

The degrees of freedom that can be leveraged to satisfy the conditions (2.17*b*) are the coefficients $c_{j-d,\zeta ,\eta }$. In fact, $\boldsymbol {r}_{j,\zeta }$ is a function of all the eigenvector corrections $\hat {\boldsymbol {p}}_{k,\zeta }$ that have been determined at orders $k<j$, and due to (2.15), it is a function of all the coefficients $c_{k,\zeta ,\eta }$. Analogous to the derivation outlined by Mensah *et al.* (Reference Mensah, Orchini and Moeck2020), it can be shown that all the coefficients with $k>j-d$ have no influence on the order $j$ conditions (2.17), and that the order $j-d$ coefficients that guarantee solvability are given by

where the terms in $\boldsymbol {r}_{j,\zeta }^{\bot }$ include all the information available at order $j$ on the eigenvectors $\hat {\boldsymbol {p}}_{k,\zeta }$ – specifically, the orthogonal components $\hat {\boldsymbol {p}}^{\bot }_{k,\zeta }$ and all the coefficients $c_{k,\zeta ,\eta }$ for $k<j-d$. A derivation of this equation in the case $d=1$, which is the most common scenario, is outlined in appendix B. The general case is treated in § 3 of the supplementary material available at https://doi.org/10.1017/jfm.2020.586. Once both the eigenvalue corrections $s_{j,\zeta }$ and the coefficients $c_{j,\zeta ,\eta }$ have been evaluated, (2.14) is guaranteed to be solvable, the eigenvector corrections $\hat {\boldsymbol {p}}_{j,\zeta }$ can be calculated and one can finally proceed to the next order.

Equation (2.18) is a theoretical contribution of this study, and is important for several reasons. It is inversely proportional to the eigenvalue split gap that occurred at order $d$; the numerator is formally equivalent to the eigenvalue correction equation, but with the adjoint eigenvector chosen to be that of the ‘other’ branch ($\eta \neq \zeta$); although it is obtained at order $j$, it contains no unknown terms at this order, and instead it defines coefficients at order $j-d$. This is consistent with the fact that the numerator is of order $ O(\varepsilon ^j)$, whereas the denominator is of order $O(\varepsilon ^d)$. As a consequence, in order to obtain perturbations accurate to order $N$, an expansion at order $N+d$ is needed. By repeatedly applying the equations contained in §§ 2.2.1 and 2.2.2, one can calculate the eigenvalue and eigenvector coefficients of power series expansions of degenerate, semi-simple eigenvalues to arbitrarily high orders.

To conclude this theoretical section, we highlight that, in most cases of practical relevance, perturbations unfold degenerate states at first order ($d=1$). This is known as complete regular splitting (Lancaster, Markus & Zhou Reference Lancaster, Markus and Zhou2004). The solution of the first-order equations ($j=1$) then follows what is described in § 2.2.1 and determines the eigenvalue splitting and the correct basis. At second order ($j=2$), the solution follows what is described in § 2.2.2, from which one can see that the expression for the eigenvalues is still exact (because no coefficients $c$ are evaluated at first order). However, also the coefficients $c_{1,\zeta ,\eta }$ need to be determined for solvability at second order; if these are ignored, all the higher-order coefficients for both eigenvalues and eigenvectors will be incorrect. Perturbations that unfold the degeneracy at first order were discussed by Magri *et al.* (Reference Magri, Bauerheim and Juniper2016), where only variations in the eigenvalues and not in the eigenvectors were considered; this explains why the perturbation theory that was outlined in that study was applicable up to second order only.

## 3. Radius of convergence and EPs

The theory introduced in § 2 yields approximations for the parametric dependence of simple and semi-simple eigenvalues and their associated eigenvectors. It provides explicit expressions for the coefficients of power series expansion up to arbitrary order. For a simple eigenvalue $s$, the function $s=s(\varepsilon )$ can always be locally expanded into a power series up to arbitrary order (Kato Reference Kato1980). For degenerate semi-simple eigenvalues, power series expansions at high orders can also generally be obtained, provided that the eigenvalue splitting is regular (Lancaster *et al.* Reference Lancaster, Markus and Zhou2004). However, regardless of the degeneracy of the eigenvalue of interest, power series expansions generally have a finite radius of convergence (Fisher Reference Fisher1999). There is, therefore, a fundamental question that needs to be addressed: In what region of the parameter space do these power series approximations of the eigenvalues and eigenvectors converge?

The limit in the convergence of a power series expansion is ruled by the closest singularity in parameter space, i.e. a point $\varepsilon _{{sng}}$ such that $s(\varepsilon _{sng})$ is singular. There may be two reasons for a singularity to exist. (i) The algebraic dependence $s(\varepsilon )$ explicitly contains a pole of the form $1/(\varepsilon -\varepsilon _{sng})$. A notable example for compressible fluid dynamics and (thermo)acoustic problems is the dependence of the governing equations on the boundary conditions expressed in terms of an impedance $Z$, which can appear at the denominator of the governing equations, and cause a pole singularity for sound soft boundary conditions, $Z=0$. (ii) For $\varepsilon =\varepsilon _{sng}$ the eigenvalue problem features a defective eigenvalue with infinite sensitivity, i.e. $\varepsilon _{sng}$ is an EP in parameter space (Kato Reference Kato1980). Fortunately, the closest singularity can be estimated directly from the power series coefficients, at no additional numerical cost. Using an approach that has been successfully applied in quantum mechanics (Fernández Reference Fernández2000), in the following we demonstrate how this can be achieved and how it enables us to identify the EPs closest to an eigenvalue of interest.

### 3.1. Estimating the radius of convergence from high-order perturbation coefficients

Close to a singularity located at $\varepsilon = \varepsilon _{sng}$ in the parameter space, the eigenvalue parameter dependence has to be of the form

where $k\in \mathbb {Q}\backslash \mathbb {N}$. If $k \in \mathbb {Z}^-$, the singularity corresponds to a pole; if $k\in \mathbb {Q}^+\backslash \mathbb {N}$, the singularity corresponds to a branch point. The values of $\varepsilon _{sng}$ and the exponent $k$ are unknown *a priori*. However, it can be shown that both quantities can be estimated from the coefficients $s_j$ of a power series that is expanded in the vicinity of (but not at) the singularity, using the relations

*a*)\begin{gather} \varepsilon_{sng}= \varepsilon_0 + \frac{s_j s_{j-1}}{(\,j+1)s_{j+1}s_{j-1}-js_j^2}, \end{gather}

*b*)\begin{gather} k=\frac{(\,j^2-1)s_{j+1}s_{j-1}-(\,js_j)^2}{(\,j+1)s_{j+1}s_{j-1}-js_j^2}. \end{gather}

We refer to the perturbation techniques explained in Fernández (Reference Fernández2000, chap. 6) for a detailed derivation of (3.2). These estimates are asymptotic, become increasingly more accurate with the perturbation order and converge to the closest singularity. This is numerically shown in § 4. When calculating the high-order sensitivity of eigenvalues around an unperturbed parameter $\varepsilon _0$, the series of eigenvalue correction coefficients will therefore converge to the actual value within a disc with radius

known as the radius of convergence. The value of $k$ aids in understanding the nature of the singularity: poles are identified for negative values of $k$, whereas EPs are identified by fractional values of $k$ of the form $1/a$, where $a$ is the algebraic multiplicity of the defective eigenvalue at the EP ($a=2$ for the cases considered in this article).

### 3.2. Locating EPs using perturbation theory

The closer the expansion point is to the singularity, the higher is the accuracy of the singularity estimated by (3.2*b*). This suggests a procedure that can be used to accurately locate EPs. Rather than performing a high-order expansion around $\varepsilon _0$, which becomes relatively time consuming at high orders, we adopt the following iterative scheme: (i) calculate the expansion coefficients $s_j$ of an eigenvalue up to about order $N=10$, using the theory of § 2; (ii) use these coefficients to estimate the closest singularity $\varepsilon _{sng}$ by means of (3.2*b*) at the highest available order; (iii) if the radius of convergence (3.3) is larger than a predefined threshold $\delta$, shift the expansion point to $\varepsilon _0 \leftarrow \varepsilon _0 + \varepsilon _{sng}$ and repeat from point (i). When $R_c < \delta$, the (shifted) expansion point coincides with the singular parameter, up to an error of order $O(\delta )$.

In general, the closest singular parameter $\varepsilon _{sng}$ will be a complex number, even if the associated physical parameter (e.g. a time delay or a length) is a real quantity. We refer to these as non-physically realizable EPs, because one cannot perform real-world experiments with complex-valued parameters. Real-valued EPs exist, but are unlikely to be found when varying only one parameter (Seyranian, Kirillov & Mailybaev Reference Seyranian, Kirillov and Mailybaev2005). A strategy to locate EPs while varying two parameters was suggested by Orchini *et al.* (Reference Orchini, Silva, Mensah and Moeck2020). Even if the singularity $\varepsilon _{sng}$ is found in the complex plane, it nonetheless limits the convergence of the power series. This also applies when considering only real values of the parameter $\varepsilon$. Furthermore, although not realizable, the presence of complex-valued singularities has an effect on the eigenvalue trajectories. In fact, in the vicinity of EPs, eigenvalues have extremely large first-order sensitivities (which become infinite at the EP). These large sensitivities cause steep variations of the eigenvalue in the complex plane, as recently analysed by Orchini *et al.* (Reference Orchini, Silva, Mensah and Moeck2020) and observed in, for example, Bauerheim *et al.* (Reference Bauerheim, Salas, Nicoud and Poinsot2014) and Sogaro *et al.* (Reference Sogaro, Schmid and Morgans2019). This phenomenon is known as mode veering (Seyranian *et al.* Reference Seyranian, Kirillov and Mailybaev2005), and is the fundamental cause of the strong sensitivity of some thermoacoustic eigenvalues, and the deviation of the eigenvalue trajectories away from the first-order sensitivity predictions.

## 4. Applications

We apply the methods developed in §§ 2 and 3 to two fundamental configurations for the investigation of thermoacoustic instabilities: a longitudinal combustor and an annular combustor. Both geometries correspond to existing experiments: respectively, the BRS combustor (Komarek & Polifke Reference Komarek and Polifke2010) and the MICCA annular combustor (Bourgouin *et al.* Reference Bourgouin, Durox, Moeck, Schuller and Candel2013). The nonlinear thermoacoustic eigenvalue problem (2.3) is solved for these configurations, using an $n$–$\tau$ model to reproduce the flame response at a frequency of interest.

### 4.1. Axial combustors: non-degenerate thermoacoustic modes

We consider a non-degenerate thermoacoustic eigenvalue in a longitudinal combustor. In addition to demonstrating the validity of non-degenerate high-order perturbation expansions, detailed in Mensah *et al.* (Reference Mensah, Orchini and Moeck2020) and summarized in (2.7), we use this simpler configuration to show how the theory of § 3 can be used to (i) quantify the convergence limit of high-order perturbation theory and (ii) identify the closest EP to an operating condition of interest, which in turn gives essential information on the thermoacoustic spectrum and the trajectory followed by the eigenvalues when a parameter is varied. The set-up we model is known as the BRS combustor; a detailed description of the geometry and the experiment is given by Komarek & Polifke (Reference Komarek and Polifke2010). The three-dimensional geometry of the model we solve is shown in figure 1. It consists of a cylindrical plenum, a premixing/swirling duct and a combustion chamber with rectangular cross-section.

The BRS combustor is one of the first configurations in which thermoacoustic instabilities at a frequency that is not directly related to an acoustic mode have been experimentally observed, at approximately 100 Hz (Tay-Wo-Chong *et al.* Reference Tay-Wo-Chong, Bomberg, Ulhaq, Komarek and Polifke2012). This instability was first generically attributed to ‘flame dynamics’. Later, it has been better understood and reproduced in a low-order network model by Emmert *et al.* (Reference Emmert, Bomberg, Jaensch and Polifke2017), and relabelled as ITA instability. Such ITA instabilities can be observed even in purely anechoic conditions, as they originate from the intrinsic feedback between the generation of acoustic waves by the flame and the sensitivity of the latter to upstream velocity fluctuations.

We discretize the thermoacoustic (2.3) on this geometry, imposing sound hard boundary conditions ($Z=\infty$) on all walls, except at the outlet, which is assumed to be sound soft ($Z=0$). A compact flame is located at the inlet of the combustion chamber. Across the flame, a temperature jump $T_2/T_1\approx 5$ is imposed. The flame response is modelled with an $n$–$\tau$ model, whose values are extracted from the flame transfer function (FTF) reported by Tay-Wo-Chong *et al.* (Reference Tay-Wo-Chong, Bomberg, Ulhaq, Komarek and Polifke2012) around a 100 Hz frequency. In particular, the time lag is assumed to be constant, as the experimentally determined FTF phase linearly decreases with frequency, with a slope $\tau \approx 4\ \textrm {ms}$. The flame gain $n$ was instead shown to be frequency dependent, with values up to 2. We choose to specify a constant value of the FTF gain, $n_0=1.68$, and consider it as a perturbation parameter.

The thermoacoustic eigenvalue problem is solved using the open-source thermoacoustic eigenvalue solver PyHoltz (PyHoltz2018). We first employ standard iterative Newton techniques to solve the eigenvalue problem. We identify a thermoacoustic mode with growth rate $\sigma =-150\ \textrm {s}^{-1}$ and frequency $f=65.2\ \textrm {Hz}$. This eigenvalue is close to that of a Helmholtz resonant mode of the combustor, in which the plenum acts as a cavity and the premixing tube as a neck (Emmert *et al.* Reference Emmert, Bomberg, Jaensch and Polifke2017). However, its eigenvector – Modeshape C in figure 1 – is not fully consistent with that of a Helmholtz mode: this mode is in fact active not only in the plenum, but also at the end of the premixing tube and at the inlet of the combustion chamber. The nature of this modeshape is clarified in the following.

By slowly decreasing the interaction index $n$ towards zero, with steps $\Delta n=-0.02$, we track the eigenvalue trajectory in the complex plane with a continuation-like method. This eigenvalue trajectory is shown in figure 2(*a*). We find that, in the limit $n\rightarrow 0$, the growth rate of this eigenvalue tends to negative infinity, and the frequency is consistent with that of the ITA mode identified by Emmert *et al.* (Reference Emmert, Bomberg, Jaensch and Polifke2017) and Orchini *et al.* (Reference Orchini, Silva, Mensah and Moeck2020). Modeshape B in figure 1 shows the magnitude of the pressure mode found when ${n=0.01}$. Its shape is indeed consistent with that of an ITA mode, as the magnitude of the eigenvector is high only around the flame (Courtine, Selle & Poinsot Reference Courtine, Selle and Poinsot2015).

We then repeat the eigenvalue tracking restarting from $n=n_0$ and increasing the interaction index $n$ up to 2.5, with steps $\Delta n=0.02$. The trajectory of the eigenvalue for large values of $n$ follows a highly nonlinear path; strong mode veering is observed. Mode veering may generally occur when a small variation in a parameter can cause two closely spaced eigenvalues to coalesce into a single degenerate eigenvalue. This degenerate eigenvalue is more likely to be an EP than a semi-simple one, unless the problem considered contains specific symmetries (Seyranian *et al.* Reference Seyranian, Kirillov and Mailybaev2005), which is not the case for the thermoacoustic system considered here. We can therefore exploit perturbation theory to identify the EP responsible for the eigenvalue veering.

We consider the interaction index $n$ as a perturbation parameter, choosing as a baseline solution the value $n_0=1.68$. We apply high-order perturbation theory for non-degenerate eigenvalues, (2.6), and calculate the coefficients of the Taylor expansion of the eigenvalue up to $30$th order. We then employ (3.2) to estimate the value of $n$ at which the closest singularity is found, $n_{sng}$, and its exponent, $k$. At the highest order considered, the exponent of this singularity is $k=0.49$, close to that of a square-root branch point. This supports the fact that the singularity is due to an EP with algebraic multiplicity $a=2$, at which $k$ should have the value $1/a$. We can also estimate the radius of convergence of the power series from (3.3). This is shown in figure 2(*b*) as a function of the perturbation order. The estimated value for $R_c$ strongly oscillates with estimates at low expansion orders, but converges to a constant value, $R_c=0.52$, at high expansion orders ($N>10$). Thus, we can determine that the Taylor expansion of the eigenvalue converges to the correct result in the entire range $n \in [1.16,2.2]$, which is a broad range for a flame gain parameter. More specifically, if we were to allow for complex values of the interaction index $n$, the eigenvalue power series expansion would converge for all $|n-n_0|<R_c$. This is shown in figure 3. The first-order sensitivity estimate (dashed line in figure 2*a*) correctly predicts the slope of the eigenvalue trajectory at the expansion point, but fails in identifying the mode veering. The trajectory reconstructed from the expansion up to $30$th order (solid black line in figure 2*a*), instead, captures the veering and is almost indistinguishable from the exact solution (thick, shaded line) within the convergence limits (black markers). Outside of the radius of convergence, the expansion is not valid, and eigenvalue estimates quickly diverge.

The value of $n$ for which the eigenproblem contains an EP is $n_{sng} = 2.19 - 0.067\mathrm {i}$, highlighted with a star in figure 3. As this value is complex, this specific EP cannot be physically realized, although it would be possible to vary a second parameter (e.g. $\tau$) to find an EP at real-valued parameters (Mensah *et al.* Reference Mensah, Magri, Silva, Buschmann and Moeck2018; Orchini *et al.* Reference Orchini, Silva, Mensah and Moeck2020). Nonetheless, the imaginary part of $n_{sng}$ is small: the smaller is the imaginary part of the closest EP, the stronger is the eigenvalue veering in its vicinity. The eigenvalue identified at the singular point is indicated in figure 2(*a*) with a star; the same figure also shows how the eigenvalue trajectory exhibits the strongest veering in the vicinity of the EP. If we calculate thermoacoustic eigenvalues for real values of $n$, we can never reach exactly the EP. However, when $n=2.19$, the eigenvalue we track is very close to being defective. This implies that another eigenvalue exists in its vicinity since these two eigenvalues must coalesce at the EP. Using the contour-integration method suggested by Buschmann *et al.* (Reference Buschmann, Mensah, Nicoud and Moeck2020*b*) for thermoacoustic eigenvalue problems, we hence search for all eigenvalues found for $n=2.19$ in the vicinity of the EP eigenfrequency. We identify two eigenvalues: one is already known, as it lies on the trajectory that was already discussed, but the second is a new eigenvalue. By applying again a continuation method, we track this newly found eigenvalue trajectory in the entire range $n\in [0,2.5]$. This trajectory is the one shown at the top left of figure 2(*a*). When $n=0$, this eigenvalue has zero growth rate, and corresponds to a purely acoustic mode, specifically, a Helmholtz mode of the plenum. Modeshape A of figure 1 shows the magnitude of this mode.

We have now all the ingredients to interpret the modeshape of the thermoacoustic mode found at $n=1.68$, Modeshape C in figure 1. Starting from small values of $n$, two thermoacoustic eigenvalues exist, with similar frequencies but different growth rates: the one with zero growth rate is of acoustic origin, the one with very negative growth rate is of intrinsic origin. Their eigenvectors are indicated as Modeshapes A and B, respectively, in figure 1. As we increase the flame gain, the two eigenvalues first approach each other towards the EP, but eventually are repelled away from it for $n>2.19$. At the EP, not only the eigenvalues, but also the two modeshapes coalesce. The eigenvalue we considered at $n=1.68$, marked with a diamond in figure 2(*a*), lies between the acoustic and the ITA eigenvalues, and it is relatively close to the veering region caused by the EP (figure 3). Thus, the modeshape of the thermoacoustic eigenvalue at $n=1.68$ contains features from both the acoustic and the ITA one. This is clearly the case for Modeshape C in figure 1: the thermoacoustic modeshape has a strong plenum component (inherited from the acoustic modeshape) but also a strong component around the flame zone (inherited from the ITA mode).

The case just discussed (i) validates high-order perturbation theory for simple eigenvalues and its convergence limits and (ii) demonstrates that knowledge of the existence of EPs is essential for understanding the structure of thermoacoustic modeshapes, as well as in the identification of other eigenvalues in the vicinity of trajectories that exhibit strong veering. We conclude the analysis with some remarks on the numerical cost of the calculations. On a quad-core Intel i7 processor, the Newton-like method employed for calculating eigenvalues without perturbation theory takes approximately 2 s to convergence for each value of the parameter $n$ considered, provided that the initial guess is reasonably accurate. Perturbation theory, on the other hand, takes approximately 30 s to calculate the expansion coefficients up to $30$th order, but can then be used to evaluate the eigenvalues accurately for any value of $n\in (n_0-R_c, n_0+R_c)$ at negligible computational cost.

### 4.2. Annular combustors: degenerate thermoacoustic modes

We now consider an annular combustor geometry, which is directly relevant for aeronautical and power generation gas turbines. Annular combustors are known to exhibit degenerate eigenvalues (Evesque, Polifke & Pankiewitz Reference Evesque, Polifke and Pankiewitz2003; Noiray & Schuermans Reference Noiray and Schuermans2013; Bothien, Noiray & Schuermans Reference Bothien, Noiray and Schuermans2015), which arise from spatial symmetries of the system (typically discrete rotational symmetry and reflection symmetry). These degenerate eigenvalues have algebraic and geometric multiplicity two, i.e. they are semi-simple. From a physical point of view, these degenerate modes can be thought of as representing two travelling waves, one spinning in the clockwise direction and one in the counterclockwise direction, at the same frequency. A pair of degenerate thermoacoustic modes interacts nonlinearly, and a mode-selection process takes place, which can lead to the stabilization of spinning, standing or mixed-type thermoacoustic oscillations (Noiray, Bothien & Schuermans Reference Noiray, Bothien and Schuermans2011; Ghirardo, Juniper & Moeck Reference Ghirardo, Juniper and Moeck2016; Laera *et al.* Reference Laera, Schuller, Prieur, Durox, Camporeale and Candel2017). All these types of oscillations have been observed experimentally (Noiray *et al.* Reference Noiray, Bothien and Schuermans2011; Worth & Dawson Reference Worth and Dawson2013; Bourgouin *et al.* Reference Bourgouin, Durox, Moeck, Schuller and Candel2014*b*). An annular combustor in which azimuthal instabilities have been investigated in detail is known as the MICCA combustor (Bourgouin *et al.* Reference Bourgouin, Durox, Moeck, Schuller and Candel2014*b*; Prieur *et al.* Reference Prieur, Durox, Schuller and Candel2017). In addition to standing and spinning oscillations, this combustor also exhibits a more complex oscillation pattern that has been labelled slanted mode, and is believed to be due to the synchronization of a longitudinal and an azimuthal thermoacoustic instability (Bourgouin *et al.* Reference Bourgouin, Durox, Moeck, Schuller and Candel2014*a*; Orchini, Mensah & Moeck Reference Orchini, Mensah and Moeck2018; Moeck *et al.* Reference Moeck, Durox, Schuller and Candel2019; Yang, Laera & Morgans Reference Yang, Laera and Morgans2019). Given the interest of the thermoacoustic community in this combustor, we demonstrate the application of perturbation theory to the MICCA combustor geometry. Our goal here is not that of accurately reproducing the dynamics observed in the MICCA combustor, and we therefore simplify the linear flame response to an $n$–$\tau$ model with constant coefficients.

The MICCA geometry consists of two annular sections, a plenum and a combustion chamber, connected by 16 ducts that are equispaced around the annulus. The geometry we model is shown in figure 4. The geometrical details used for modelling this configuration are given by Mensah *et al.* (Reference Mensah, Magri, Orchini and Moeck2019). We impose sound hard boundary conditions ($Z=\infty$) on all walls, and a sound soft boundary condition ($Z=0$) at the outlet of the combustion chamber. An axially varying temperature field is used (see figure 4), as in Laera *et al.* (Reference Laera, Schuller, Prieur, Durox, Camporeale and Candel2017), with the speed of sound fixed to $c=348\ \textrm {m}\,\textrm {s}^{-1}$ in the cold plenum and ranging from $c=784\ \textrm {m}\,\textrm {s}^{-1}$ at the flame zone to $c=690\ \textrm {m}\,\textrm {s}^{-1}$ at the chamber exit. An acoustically compact flame zone is located at the exit of each duct. We fix the gain of the flame response to $n=1$, and we consider the flame time delay as the perturbation parameter, starting from the unperturbed value $\tau _0=3\ \textrm {ms}$. Variations in the flame's time delay response are known to have a strong impact on the thermoacoustic stability (Lord Rayleigh Reference Rayleigh1878; Dowling Reference Dowling1995), which can be achieved, for example, by means of fuel staging or the use of cylindrical burner outlets (Krebs *et al.* Reference Krebs, Flohr, Prade and Hoffmann2002; Noiray *et al.* Reference Noiray, Bothien and Schuermans2011).

We discretize the thermoacoustic equations on the MICCA geometry with finite elements, and we solve the resulting eigenproblem using PyHoltz. For the set of parameters we investigate, we identify an unstable thermoacoustic mode with frequency $f=527.78\ \textrm {Hz}$ and growth rate $\sigma =306.44\ \textrm {s}^{-1}$. This eigenvalue is degenerate, with algebraic multiplicity $a=2$, and it is associated with a mode of azimuthal order 1 of the plenum cavity. This degeneracy is due to the rotation and mirror symmetries of the combustor. Thus, the eigenvalue is semi-simple, and the geometric multiplicity also equals 2. There exist therefore two linearly independent eigenvectors, spanning the two-dimensional subspace associated with the degenerate eigenvalue, which we calculate together with their corresponding adjoint eigenvectors.

We consider two types of perturbation patterns: (I) we perturb the time delay $\tau$ of all flames; (II) we perturb the time delay of certain flames only, specifically flames 1, 4, 8 and 10, counting in the counterclockwise direction (see figure 5). While the former pattern preserves all the symmetries of the combustor geometry, the latter pattern is chosen such that both the mirror and discrete rotational symmetries are broken. For each of the two cases, we apply degenerate perturbation theory as discussed in § 4.2.

Since pattern I preserves the symmetry, the degeneracy has not unfolded. Regardless of the magnitude of the perturbation, the eigenvalue remains degenerate with algebraic and geometric multiplicity 2. The perturbation method correctly identifies that the eigenvalue does not split since the eigenvalues of the auxiliary eigenvalue problem (2.11) remain degenerate at all considered orders. We therefore have only one power series expansion, whose radius of convergence (shown as percentage variation from the unperturbed time delay $\tau _0$) is shown in figure 5 at various perturbation orders. The results show that high-order perturbation theory accurately predicts the evolution of the degenerate eigenvalue from the unperturbed state for variations in $\tau$ up to $22\,\%$, a relatively large value for a flame time delay. The subspace spanned by the two linearly independent eigenvectors varies as a function of the perturbation parameter, according to (2.4*a*,*b*), with the coefficients of the eigenvector expansion calculated according to (2.15). Note that, since the problem remains degenerate, any other linear combination of the eigenvectors identified by our method would be equally valid. Our method, however, always converges to the same solution, which is constrained by the bi-orthogonalization of the perturbed eigenvectors.

On the other hand, the perturbation pattern II unfolds the degeneracy as the chosen flame staging pattern breaks the spatial symmetries. The two eigenvalues therefore follow different trajectories when varying $\tau$ (shown in figure 6). We label these eigenvalue trajectories branch 1 and branch 2. Since their power series expansions have different coefficients, the estimated radius of convergence and closest EP, calculated via (3.2) and (3.3), will be different for the two branches. The radii of convergence estimated from the coefficients of the two eigenvalue expansions are shown in figure 5(*b*). Note that, for branch 2, there is a small mismatch between the radius of convergence estimated via the high-order perturbation theory expansion coefficients (dotted line) and the actual distance to the closest branch point obtained using the iterative procedure outlined in § 3.2 (shaded line). This deviation is justified in that (3.2*a*) is only an estimate of the location of the EP, even at high orders: the closer is the expansion point to the EP, the better is the estimate of the radius of convergence.

Despite this small mismatch, figure 5 indicates that perturbation theory should yield correct results in the prediction of the evolution of the eigenvalues in a significant range of values of $\tau$, which can vary up to $15\,\%$ and $20\,\%$ for the two eigenvalues, respectively. This is verified by comparing the reconstruction of the eigenvalue trajectories from the high-order perturbation theory with brute-force solutions (direct numerical solutions of the thermoacoustic eigenproblem (2.3)) at various values of $\tau \in [2.3,3.6]\ \textrm {ms}$, shown in figure 6. The shaded backgrounds in the figure indicate the range of convergence of each eigenvalue. It is evident that, within each convergence region, the brute-force and high-order perturbation results are almost indistinguishable, for both frequency and growth rate of the eigenvalues. On the other hand, as soon as perturbation theory is applied outside of the radius of convergence of its corresponding eigenvalue, the eigenvalue power series approximation rapidly diverges from the actual solution. For comparison, we have also reported in figure 6 the eigenvalue approximation that one would obtain by using only first-order sensitivity. For both growth rate and frequency, first-order perturbation theory correctly predicts the tangential direction along which the eigenvalues split, but fails in capturing the more complex, nonlinear behaviour of the eigenvalues at moderate changes of the perturbation parameter. This is true particularly for the growth rates shown in figure 6. First-order theory predicts that increasing (decreasing) the time delay will linearly increase (decrease) the growth rate of both eigenvalues. High-order theory, instead, shows that the growth rate of the unperturbed eigenvalue is close to a maximum, so that either increasing or decreasing the flame time lag response will lead to a reduction in the growth rate of both modes.

Convergence results analogous to those of the eigenvalues are also obtained for the eigenvectors. Since perturbation pattern II breaks the symmetries, there exists a specific basis in the unperturbed degenerate subspace that perturbation theory must identify, which consists of the pair of eigenvectors $\hat {\boldsymbol {p}}_{0,1}$ and $\hat {\boldsymbol {p}}_{0,2}$ into which the degeneracy linearly unfolds. This is also correctly captured by the presented high-order perturbation theory, which is able to accurately reconstruct the three-dimensional modeshape of the thermoacoustic modes at hand when varying $\tau$. As was shown in Mensah *et al.* (Reference Mensah, Orchini and Moeck2020), the accuracy of the estimated eigenvectors within the radius of convergence increases at higher perturbation orders. Figure 7 shows a plenum cut view of the magnitude of the pressure modeshapes of the split eigenvalues for $\tau =2.65\ \textrm {ms}$. Although the asymmetry is moderate – 4 flames have a time delay that differs by 0.35 ms from that of the other 12 flames – the modeshapes deviate markedly from the symmetric case. This is particularly evident for the lower-frequency mode (on the top in figure 7), for which the pressure nodes are visibly not aligned, and the two pressure maxima have different intensities. The modeshapes reconstructed using perturbation theory and the relative error between the exact and approximated eigenvectors are also shown. We recall that eigenvectors can always be arbitrarily scaled by a complex-valued coefficient. Thus, when comparing the exact and approximated eigenvectors, one needs to ensure that the same normalization has been applied to the eigenvectors calculated with the exact and approximated method. This was discussed in Mensah *et al.* (Reference Mensah, Orchini and Moeck2020).

We conclude the analysis with some numerical remarks. When two closely spaced eigenvalues exist, as those generated by symmetry breaking, it may be difficult to identify both of them using standard Newton-like algorithms. This is because these iterative algorithms need an initial guess and will eventually converge to an eigenvalue solution. We observe that, often, the basin of attraction of one of these two solutions is much larger than that of the other. Unless rather accurate initial guesses are provided to the iterative algorithm, only one of the two solutions will be identified. High-order perturbation methods are useful also in this respect. Within the radius of convergence, the eigenvalue estimates they provide are suitable initial guesses for the actual eigenvalues of the system, as shown above, and Newton-like methods are able to identify both closely spaced eigenvalues in a few iterations. The brute-force results shown in figure 6 have been calculated using this effective strategy. An alternative is to use Beyn's contour-integration method, which is able to identify all eigenvalues inside a circle in the complex-frequency space (Beyn Reference Beyn2012; Buschmann *et al.* Reference Buschmann, Mensah, Nicoud and Moeck2020*b*). Both methods have their advantages and disadvantages. High-order perturbation methods require only one longer calculation, needed to obtain the eigenvalue expansion coefficients, and then only short calculations to (i) obtain good initial guesses and (ii) converge to the eigenvalue(s) of interest via Newton methods, for any value of the perturbation parameter within the radius of convergence. Beyn's method, instead, requires a moderately long calculation for each of the perturbation parameter values in which one is interested, but does not suffer from radius of convergence limitations – although for very large perturbations, the split eigenvalues may be far from each other, which requires the integration over a large circle, with increasing computational effort needed. We are of the opinion that no general conclusion can be made regarding which of the two methods is preferable. The trade-off depends on the size of the eigenproblem at hand and on the range of parameters one wants to investigate.

## 5. Expansion of defective eigenvalues at EPs

At defective points, the theory described in § 2 breaks down. This is because the eigenvectors of a defective eigenvalue lose the bi-orthogonality property (2.8). We have demonstrated the existence of EPs, which are defective eigenvalues (§ 3), in the spectrum of thermoacoustic systems. To make our technique applicable to all types of eigenvalues, in this section, we present the perturbation theory at defective eigenvalues. For defective eigenvalues, the following relation between direct and adjoint eigenvectors holds in place of (2.8):

This is known as self-orthogonality in non-Hermitian quantum mechanics (Moiseyev Reference Moiseyev2011; Heiss Reference Heiss2012) and invalidates the derivation of the expansion equations of § 2.2. The latter relies on expressions of the form analogous to those shown in (2.7), in which the scalar product in the denominator vanishes in the defective eigenvalue scenario due to the self-orthogonality condition (5.1). Consequently, already the first-order sensitivity of the eigenvalues diverges to infinity, and it is not possible to expand the dependence of the eigenvalues on a parameter into a power series.

However, eigenvalues at an EP can be expanded into a fractional power series, also known as Puiseux series. In particular, for a defective eigenvalue with algebraic multiplicity $a$, it is possible (under mild assumptions) to expand the eigenvalue as follows (Leung Reference Leung1990; Lancaster *et al.* Reference Lancaster, Markus and Zhou2004):

By using this ansatz, and introducing the concept of generalized eigenvectors for defective eigenvalues (see supplementary material § 1), it is possible to follow the steps of § 2.2 to derive arbitrary high-order equations for the calculation of the coefficients of the Puiseux series (5.2). It is not the aim of this contribution to provide the entire (and lengthy) derivation of these expressions. Our goal is to show that expansions at EPs are possible, and that the Puiseux coefficients can be evaluated by means of adjoint methods. Additional details of the derivation of the Puiseux expansions at EPs are provided in the supplementary material.

Here, we focus only on the special case of defective eigenvalues with algebraic multiplicity $a=2$ and (consequently) geometric multiplicity 1. In this case, it can be shown (see (S1.8) in the supplementary material) that the first coefficient of the Puiseux series expansion is given by

where $\hat {\tilde {\boldsymbol {p}}}_{gen}$ is the generalized eigenvector, defined by $\boldsymbol{\mathsf{L}}_{0,0}\hat {\tilde {\boldsymbol {p}}}_{gen} \equiv -\boldsymbol{\mathsf{L}}_{1,0}\hat {\boldsymbol {p}}_{{def}}$. Note the differences between (5.3) and the first-order ordinary sensitivity equation, which reads $s_1 = - \left . \left \langle \left .\hat {\tilde {\boldsymbol {p}}}^{\dagger }_0\right |\boldsymbol{\mathsf{L}}_{0,1} \hat {\tilde {\boldsymbol {p}}}_0\right \rangle \right /\left \langle \left .\hat {\tilde {\boldsymbol {p}}}^{\dagger }_0\right |\boldsymbol{\mathsf{L}}_{1,0} \hat {\tilde {\boldsymbol {p}}}_0\right \rangle$. First, two branches (the + and $-$ roots) stem from the sensitivity equation for defective eigenvalues (5.3) due to the $a=2$ algebraic multiplicity. Second, a square root appears in the equation, highlighting the fact that the defective eigenvalue is a branch point. This is due to the fact that the first coefficient in the Puiseux series expansion is determined at second order (whereas first-order power series expansion coefficients are determined at first order). Lastly, a new term appears in the denominator, which involves the second derivative of the operator $\boldsymbol{\mathsf{L}}$ with respect to the eigenvalue. This results from the fact that the defective eigenvalue sensitivity is determined at second order and that bi-orthogonality conditions do not hold.

The above sensitivity equation represents an exception to the standard sensitivity equations for semi-simple eigenvalues outlined in § 2.2, and should be used if and only if the eigenvalue of interest is a defective eigenvalue. Note that, in practice, it is unlikely to converge exactly to a defective eigenvalue with numerical methods. Both the series (3.2*a*) – which converges to the closest EP in the complex-parameter space – and the method outlined in Orchini *et al.* (Reference Orchini, Silva, Mensah and Moeck2020) – which converges to an EP in the real parameter space – are iterative numerical methods, whose accuracy is limited by machine precision. Due to the infinite sensitivity of eigenvalues at EPs, even a very small deviation in the estimation of the defective eigenvalue is sufficient to strongly affect the eigenvalue sensitivity, causing large errors in the results of (5.3). The latter should therefore be applied only to problems for which EPs can be identified analytically.

### 5.1. Application to a one-dimensional thermoacoustic model

We now provide an application of Puiseux series expansion to numerically show that the sensitivity to small perturbations of eigenvalues at EPs is not polynomial, but scales with powers of $\varepsilon ^{1/2}$. The model we use is a single-mode Galerkin expansion of the thermoacoustic equations in an acoustically open tube (Juniper & Sujith Reference Juniper and Sujith2018), expressed in non-dimensional units, for which the scalar operator $L$ is given by

We wish to find defective eigenvalues of the above equation, which are obtained when

*a*,

*b*)\begin{equation} L(s) = 0 \quad \text{and} \quad \frac{\partial L}{\partial s} \equiv L_{1,0} = 0. \end{equation}

These are generally only necessary conditions for a defective point, as they would also be satisfied for a semi-simple degenerate eigenvalue with algebraic multiplicity of (at least) two. However, as shown in Seyranian *et al.* (Reference Seyranian, Kirillov and Mailybaev2005), when degenerate eigenvalues arise the occurrence of a defective eigenvalue is almost certain in a system without symmetries – semi-simple degenerate eigenvalues exist, but have a negligible measure compared to that of the defective ones. No symmetries are present here, so (5.5*a*,*b*) effectively identify EPs.

From (5.4) and (5.5*a*,*b*) it can be shown that choosing $\beta \tau e = {\pm }1$ for ${\rm \pi} \tau > 1$, and with the additional constraints that $\tan {[({\rm \pi} \tau )^2-1]^{1/2}}=[({\rm \pi} \tau )^2-1]^{1/2}$, yields a defective eigenvalue with algebraic multiplicity 2 and geometric multiplicity 1 of the form

The non-dimensional parameters at which we identify an EP are $\tau =1.4653$ and $\beta =0.2511$; the corresponding defective eigenvalue is $s_{def}=-0.6825 + 3.0666\mathrm {i}$.

Let us now consider a perturbation expansion of the defective eigenvalue around $\tau$. Since we have imposed that $L_{1,0} = 0$, it is not correct to use (2.7) to try to identify polynomial coefficients, as the latter equation diverges. Instead, we can calculate the first coefficient of the Puiseux series using (5.3), which yields the following expansion for the eigenvalues around the EP:

Because this example is one-dimensional, the direct and adjoint eigenvectors associated with the EP equal 1, and the generalized eigenvector vanishes.

Figure 8(*a*) shows the evolution of the eigenvalues around the EP while varying $\tau$. The solid lines are obtained by repeatedly applying Beyn's contour-integral method while varying the flame's time delay. The colour indicates the evolution of each tracked eigenvalue: at the EP (black dot) the eigenvalues coalesce, and the trajectories have a discontinuous first-order derivative, with cusp angles of $90^{\circ }$. The eigenvalue trajectories approximated with the first-order Puiseux series (5.7) are shown with dashed lines. The Puiseux series expansions accurately capture the angle along which the two cusps formed by the eigenvalue trajectories intersect at the EP, and are locally tangent to the exact trajectories. Figure 8(*b*) shows the relative error between the actual eigenvalue and that estimated by the Puiseux series for the red branch only; analogous results are obtained for the other branch. At expansion order $N=0$, the error equals the distance between the eigenvalue for a given $\tau$ and the eigenvalue found at the EP. For very small variations of $\tau$, machine precision limits the error to a lower bound. However, when increasing the variation in the time delay, $\Delta \tau$, the relative error follows a linear trend in log–log scale. Linear regression identifies the slope of this line to be $\alpha =0.501$, implying that the error scales with $\Delta \tau ^{1/2}$. This is the leading term in the Puiseux series (5.2) for algebraic multiplicity $a=2$. When calculating the error including the first term of the Puiseux expansion, instead, the relative error generally decreases and its slope becomes $\alpha \approx 1$, which is the power of the leading error in (5.7). A further increase in accuracy by a factor of $\Delta \tau ^{1/2}$ is observed when including also the second term of the Puiseux series, $N=2$, whose expression for the model investigated here is provided in (S2.4) of the supplementary material. The trajectories of the eigenvalues reconstructed including also the second term of the Puiseux series are shown in figure 8(*a*) with dotted lines. They provide better approximations of the actual eigenvalues for larger deviations from the EP, and they are able to capture both the local angle and also the curvature of the cusps that intersect at the EP.

## 6. Conclusions

In this study, we have established interconnections between several topics that are essential for the investigation of thermoacoustic stability, namely symmetry-induced degeneracies, high-order adjoint perturbation theory, the origin of thermoacoustic modes and EPs. We have first extended adjoint-based high-order perturbation formulae for thermoacoustic eigenproblems, currently available only for simple eigenvalues, to semi-simple degenerate eigenvalues, typical of annular configurations. Although degenerate perturbation theory could be extended to arbitrary level of degeneracies, the twofold-degenerate case, which arises due to spatial symmetries of combustor configurations, is the most relevant for applications in thermoacoustics, and is the one that has been considered in this study.

Well-defined boundaries of validity of eigenvalues reconstructed with perturbation methods can be defined by evaluating the radius of convergence of the power series, directly from their coefficients. The radius of convergence is generally finite but not necessarily small. In fact, it can be sufficiently large to use perturbation theory for investigating the effects of parametric variations on the thermoacoustic spectrum in a broad parameter range of physical relevance. The existence of a finite limit of convergence is due to singularities in the thermoacoustic eigenproblems that arise at EPs. Exceptional points appear in the spectrum of thermoacoustic systems when eigenvalues of acoustic and/or ITA origin coalesce, and we have shown how they can be identified by means of the coefficients of high-order perturbation methods. At EPs, eigenvalues have infinite sensitivities; moreover, the eigenvalue trajectories that form when a parameter is varied exhibit strong veering in the vicinity of an EP. This explains the high sensitivity of some thermoacoustic modes that has recently been discussed in the literature.

We applied the presented theory to two thermoacoustic configurations that model existing experiments. We demonstrated that the identification of EPs facilitates the prediction and understanding of the trajectories and modeshapes of thermoacoustic eigenvalues as the system parameters are varied. We have shown that degenerate perturbation theory is capable of correctly predicting whether a given perturbation pattern breaks or retains the degeneracy of semi-simple eigenvalues found in annular combustors. When splitting occurs, the method correctly predicts the two trajectories that are followed by both eigenvalues, as well as the eigenvector basis into which the unperturbed problem unfolds. All these results are valid within a radius of convergence that was accurately calculated from the power series coefficients. The framework enables one to predict the strong veering that thermoacoustic eigenvalues may experience, and explain the latter via the identification of the closest EP.

Lastly, we have shown how adjoint-based methods allow also for the calculation of the coefficients of Puiseux series, which can be used to expand the eigenvalues found at EPs. Combined together, our findings provide a comprehensive and efficient theory that can be applied to investigate the thermoacoustic spectrum in the vicinity of any eigenvalue (simple, semi-simple or defective) while varying a parameter.

Our results show that perturbation theory is a powerful tool that enables the efficient and accurate assessment of the effects that parametric variations have on the thermoacoustic spectrum. Its results are valid in a finite but broad range of design parameters, and can be exploited in the execution of tasks relevant to the thermoacoustic community such as design optimization and uncertainty quantification. In fact, although the trajectories followed by thermoacoustic eigenvalues are nonlinear, in particular in the vicinity of an EP, often only the linear first-order sensitivity is considered for the prediction of eigenvalue variations. This can lead to a misinterpretation of the behaviour of the thermoacoustic spectrum. The high-order adjoint-based method proposed here constitutes a significant step forward, as it enables the accurate estimation of the nonlinear behaviour of thermoacoustic eigenvalues in a wide parametric range. This helps in understanding and predicting the effect that real-world finite perturbations have on thermoacoustic stability. It can also be leveraged in advanced optimization methods, by providing nonlinear descent information that can lead to a quicker identification of local optima.

## Acknowledgements

A.O. acknowledges support from the Alexander von Humboldt Foundation and the German Research Foundation (DFG project no. 422037803). L.M. acknowledges support from the Royal Academy of Engineering Research Fellowship Scheme, and support from the Technical University of Munich Institute for Advanced Study, funded by the German Excellence Initiative and the European Union Seventh Framework Programme under grant agreement no. 291763.

## Declaration of interest

The authors report no conflict of interest.

## Supplementary material

Supplementary material is available at https://doi.org/10.1017/jfm.2020.586.

## Appendix A. Explicit expressions of perturbation equations

A detailed derivation of (2.6) has been provided in Mensah *et al.* (Reference Mensah, Orchini and Moeck2020), to which we refer the interested reader. For completeness, we report here without demonstration the explicit expression that can be used to calculate the vectors $\boldsymbol {r}_j$ at any order:

where we have introduced the multi-index notation $\boldsymbol{\mu} \equiv [\mu _1,\mu _2,\ldots ,\mu _N]$, where $\mu _n \in \mathbb {N}$, and we have defined the multi-index $\boldsymbol{1}_j$ as having index 1 at position $j$, and 0 otherwise. The following notation for multi-index properties holds:

*a*–

*d*)\begin{equation} |\boldsymbol{\mu}|\equiv\sum_{n=1}^N \mu_n, \quad |\boldsymbol{\mu}|_{w}\equiv\sum_{n=1}^N n\mu_n, \quad \binom{|\boldsymbol{\mu}|}{\boldsymbol{\mu} }\equiv\frac{|\boldsymbol{\mu}|!}{\displaystyle{\prod_{n=1}^N\left(\mu_n!\right)}}, \quad \boldsymbol{s}^{\boldsymbol{\mu}} \equiv \prod_{n=1}^N s_n^{\mu_n}. \end{equation}

From (A 1), the explicit expressions of $\boldsymbol {r}_j$ for the first two orders are

*a*)\begin{gather} \varepsilon: \quad \boldsymbol{r}_1 \equiv \boldsymbol{\mathsf{L}}_{0,1}\hat{\boldsymbol{p}}_0, \end{gather}

*b*)\begin{gather} \varepsilon^2: \quad \boldsymbol{r}_2 \equiv \left(\boldsymbol{\mathsf{L}}_{0,1} + s_1 \boldsymbol{\mathsf{L}}_{1,0}\right)\hat{\boldsymbol{p}}_1 + \left(\boldsymbol{\mathsf{L}}_{0,2} + s_1\boldsymbol{\mathsf{L}}_{1,1} + s_1^2\boldsymbol{\mathsf{L}}_{2,0}\right)\hat{\boldsymbol{p}}_0. \end{gather}

We lastly report the explicit expression for $\boldsymbol{\mathsf{L}}_{1,0}$, which is extensively used to define the bi-orthogonal conditions (2.8). From (2.5), in the case in which the operator $\boldsymbol{\mathsf{L}}$ is defined via (2.3), we have

where $s_0$ is an eigenvalue of $L$. Note that $\boldsymbol{\mathsf{L}}_{1,0}$ is actually a discretization of the expression on the right-hand side of (A 4). We again emphasize that the theory discussed in this study is valid not only for the thermoacoustic case, for which we have reported the explicit expressions (2.3) and (A 4), but also holds for arbitrary finite-dimensional operators $\boldsymbol{\mathsf{L}}$.

## Appendix B. Derivation of eigenvector coefficient equations

In this appendix we prove (2.18) for the specific case in which the degeneracy unfolds at first order, $d=1$, which is the relevant scenario for the applications shown in the present work. The equation is, however, more general and holds also for $d>1$, a proof of which is contained in § 3 of the supplementary material.

If the degeneracy unfolds at first order, the auxiliary eigenvalue problem $\boldsymbol{\mathsf{X}}_1\boldsymbol {\alpha } = s_1\boldsymbol {\alpha }$ has two different simple eigenvalues, $s_{1,\zeta }$ and $s_{1,\eta }$ – see (2.11). Furthermore, the coefficients $\boldsymbol {\alpha }$ identify the basis into which the degeneracy unfolds. In this basis, the auxiliary eigenvalue problem is diagonal, so that

where we used (A 3*a*) for the definition of $\boldsymbol {r}_{1,\zeta }$.

Then the eigenvector correction at first order is

analogous to (2.15).

Using the expression for $\boldsymbol {r}_2$, (A 3*b*), on branch $\zeta$, the perturbation equation at second order is

To be solvable, this equation needs to satisfy the conditions (2.17). We now focus only on the second of these solvability conditions, (2.17*b*), for which $\eta \neq \zeta$, and we use it to derive an equation for the coefficients $c_{1,\zeta ,\eta }$, (2.18). This condition explicitly reads

We immediately notice that the last term vanishes since $\left \langle \left .\hat {\boldsymbol {p}}_{0,\eta }^{\dagger }\right |\boldsymbol{\mathsf{L}}_{1,0} \hat {\boldsymbol {p}}_{0,\zeta } \right \rangle = 0$ for $\eta \neq \zeta$, from the bi-orthogonality conditions (2.8). We then expand the term containing $\hat {\boldsymbol {p}}_{1,\zeta }$ using its explicit expression, from (B 2). Equation (B 4) then becomes

We note that the right-hand-side argument of the inner product in the first line of this equation is formally identical to the definition of $\boldsymbol {r}_2$ from (A 3*b*), with the difference that $\boldsymbol {p}_1$ is replaced by $\boldsymbol {p}_1^{\bot }$. This suggests the definition of

We then tackle the inner product multiplying $c_{1,\zeta ,\zeta }$ in (B 5). Using the linearity of the inner product, we can split it into two parts, both of which vanish:

The first term vanishes because of the diagonalization of the basis (B 1), and the second term vanishes because of the bi-orthogonality condition (2.8). This proves that the coefficients $c_{1,\zeta ,\zeta }$ do not affect the solution of the perturbation equations at the next orders. These coefficients can be uniquely determined if a normalization condition is imposed on the eigenvector; however, this is not discussed here. Finally, we consider the inner product multiplying $c_{1,\zeta ,\eta }$ in (B 5), which simplifies to

where we have used again (B 1) for the first term, and the bi-orthogonality condition (2.8) on branch $\eta$ for the second term.

Substituting the simplified expressions (B 6)–(B 8) into (B 5), we finally have

from which an explicit expression for the coefficients $c_{1,\zeta ,\eta }$ follows:

which is (2.18) for $d=1$ and $j=1$. This proof can be extended to higher orders $n$ with appropriate definitions of the vectors $\boldsymbol {r}_{n,\zeta }^{\bot }$, which, analogously to (B 6), are obtained from the definition of $\boldsymbol {r}_{n,\zeta }$ by replacing $\hat {\boldsymbol {p}}_{n-1,\zeta }$ with $\hat {\boldsymbol {p}}_{n-1,\zeta }^{\bot }$.