Degenerate perturbation theory in thermoacoustics: high-order sensitivities and exceptional points

Abstract

In this study, we connect concepts that have been recently developed in thermoacoustics, specifically (i) high-order spectral perturbation theory, (ii) symmetry-induced degenerate thermoacoustic modes, (iii) intrinsic thermoacoustic modes and (iv) exceptional points. Their connection helps gain physical insight into the behaviour of the thermoacoustic spectrum when parameters of the system are varied. First, we extend high-order adjoint-based perturbation theory of thermoacoustic modes to the degenerate case. We provide explicit formulae for the calculation of the eigenvalue corrections to any order. These formulae are valid for self-adjoint, non-self-adjoint or even non-normal systems; therefore, they can be applied to a large range of problems, including fluid dynamics. Second, by analysing the expansion coefficients of the eigenvalue corrections as a function of a parameter of interest, we accurately estimate the radius of convergence of the power series. Third, we connect the existence of a finite radius of convergence to the existence of singularities in parameter space. We identify these singularities as exceptional points, which correspond to defective thermoacoustic eigenvalues, with infinite sensitivity to infinitesimal changes in the parameters. At an exceptional point, two eigenvalues and their associated eigenvectors coalesce. Close to an exceptional point, strong veering of the eigenvalue trajectories is observed. As demonstrated in recent work, exceptional points naturally arise in thermoacoustic systems due to the interaction between modes of acoustic and intrinsic origin. The role of exceptional points in thermoacoustic systems sheds new light on the physics and sensitivity of thermoacoustic stability, which can be leveraged for passive control by small design modifications.

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 2003; Dowling & Morgans 2005; Lieuwen & Yang 2005; Culick 2006; Poinsot 2017). 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 2005). 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. 2007). 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 2020).

1.1. Thermoacoustic eigenvalues: classification and origin

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

(1.1)\begin{equation} \boldsymbol{\mathsf{L}}(s) \hat{\boldsymbol{p}} = \boldsymbol{0}, \end{equation}

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 1968; Mehrmann & Voss 2004; Betcke et al. 2013; Güttel & Tisseur 2017). 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 2004). As recently shown, these spectral singularities are general features of thermoacoustic systems (Mensah et al. 2018; Orchini et al. 2020). 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. 2014; Bomberg, Emmert & Polifke 2015). 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. 2015; Hoeijmakers et al. 2016) and reflective environments (Emmert et al. 2017; Mukherjee & Shrira 2017; Silva et al. 2017; Buschmann, Mensah & Moeck 2020a; Orchini et al. 2020). 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 2018; Magri 2019). For the accurate calculation of these sensitivities, adjoint methods proved to be efficient mathematical and computational tools, as reviewed by Magri (2019). Adjoint methods for thermoacoustic eigenvalue sensitivity analysis were developed for design parameter and passive control by Magri & Juniper (2013), and subsequently applied to more complex flames in Magri & Juniper (2014); Orchini & Juniper (2016). Rigas et al. (2015) 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 2017).

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. 2007). 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. 2020). 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 2012; Miri & 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 2019; Sogaro, Schmid & Morgans 2019). 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 1980). 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. 2018). Orchini et al. (2020) 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 (2019) 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

(2.1)\begin{equation} \boldsymbol{\mathsf{L}}(s,\varepsilon)\hat{\boldsymbol{p}} = \boldsymbol{0} \end{equation}

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

(2.2)\begin{equation} \left\langle\left.\boldsymbol{g}\right|\boldsymbol{\mathsf{L}} \boldsymbol{f}\right\rangle \equiv \left\langle\left.\boldsymbol{\mathsf{L}}^{\dagger} \boldsymbol{g}\right|\boldsymbol{f}\right\rangle, \end{equation}

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 2007; Güttel & Tisseur 2017) 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 2003; Culick 2006; Nicoud et al. 2007)

(2.3)\begin{equation} \left[\boldsymbol{\nabla}\boldsymbol{\cdot}(c^2\boldsymbol{\nabla}) - s^2 -\frac{(\gamma - 1)}{\bar{\rho}} \frac{\bar{Q}}{\bar{U}}n \textrm{e}^{-s\tau}\hat{\boldsymbol{n}}_{{ref}} \boldsymbol{\cdot} \boldsymbol{\nabla}_{{ref}}\right]\hat{p}=0, \end{equation}

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 1965). 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. 2007), 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. (2019) 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

(2.4a,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.4a,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. (2020). 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

(2.5)\begin{equation} \boldsymbol{\mathsf{L}}_{n,m} \equiv \frac{1}{n!m!}\left.\frac{\partial^{n+m} \boldsymbol{\mathsf{L}}}{\partial s^n\partial \varepsilon^m}\right|_{\substack{s=s_0\\\varepsilon=0}}. \end{equation}

The power series approximations (2.4a,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:

(2.6)\begin{equation} \boldsymbol{\mathsf{L}}_{0,0}\hat{\boldsymbol{p}}_j = -\boldsymbol{r}_j -s_j\boldsymbol{\mathsf{L}}_{1,0}\hat{\boldsymbol{p}}_0, \quad \text{for } j=1,\ldots,N. \end{equation}

We refer to Mensah et al. (2020) 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. 2020). 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$:

(2.7)\begin{equation} s_j = - \frac{\left\langle\left.\hat{\boldsymbol{p}}_0^{\dagger}\right|\boldsymbol{r}_j\right\rangle}{\left\langle\left.\hat{\boldsymbol{p}}_0^{\dagger}\right|\boldsymbol{\mathsf{L}}_{1,0}\hat{\boldsymbol{p}}_0\right\rangle}. \end{equation}

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 2016). 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. 2016; Mensah et al. 2019). 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 2017)

(2.8)\begin{equation} \left\langle\left.\hat{\tilde{\boldsymbol{p}}}^{\dagger}_{0,\zeta}\right|\boldsymbol{\mathsf{L}}_{1,0} \hat{\tilde{\boldsymbol{p}}}_{0,\eta}\right\rangle = \delta_{\zeta,\eta}, \end{equation}

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

(2.9)\begin{equation} \hat{\boldsymbol{p}}_k = \alpha_1 \hat{\tilde{\boldsymbol{p}}}_{k,1} + \alpha_2 \hat{\tilde{\boldsymbol{p}}}_{k,2}, \end{equation}

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

(2.10a)\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}

(2.10b)\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

(2.11)\begin{equation} \boldsymbol{X}_j \boldsymbol{\alpha} = s_j \boldsymbol{\alpha}, \end{equation}

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:

1. (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.

2. (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.4a,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

(2.15)\begin{equation} \hat{\boldsymbol{p}}_{j,\zeta} = \hat{\boldsymbol{p}}^{\bot}_{j,\zeta} + c_{j,\zeta,1} \hat{\boldsymbol{p}}_{0,1} + c_{j,\zeta,2} \hat{\boldsymbol{p}}_{0,2}, \end{equation}

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

(2.16)\begin{equation} \left\langle\left.\hat{\boldsymbol{p}}_{0,\eta}^{\dagger}\right|\boldsymbol{r}_{j,\zeta}\right\rangle + \left\langle\left.\hat{\boldsymbol{p}}_{0,\eta}^{\dagger}\right|s_{j,\zeta}\boldsymbol{\mathsf{L}}_{1,0}\hat{\boldsymbol{p}}_{0,\zeta}\right\rangle = 0, \quad \text{for } \zeta=1,2 \text{ and } \eta = 1,2. \end{equation}

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

(2.17a)\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}

(2.17b)\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.17a) 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.17b) 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. 2016) 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. (2019) 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.17b) 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. (2020), 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

(2.18)\begin{equation} c_{j-d,\zeta,\eta} = \frac{\left\langle\left.\hat{\boldsymbol{p}}^{\dagger}_{0,\eta}\right|\boldsymbol{r}_{j,\zeta}^{\bot}\right\rangle}{s_{d,\eta}-s_{d,\zeta}} \quad \text{for } \eta\neq \zeta, \end{equation}

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 2004). 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. (2016), 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 1980). 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. 2004). However, regardless of the degeneracy of the eigenvalue of interest, power series expansions generally have a finite radius of convergence (Fisher 1999). 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 1980). 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 2000), 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

(3.1)\begin{equation} s(\varepsilon)\sim(\varepsilon-\varepsilon_{sng})^k, \end{equation}

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

(3.2a)\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}

(3.2b)\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 (2000, 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

(3.3)\begin{equation} R_c \equiv |\varepsilon_{{sng}} - \varepsilon_0|, \end{equation}

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.2b). 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.2b) 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 2005). A strategy to locate EPs while varying two parameters was suggested by Orchini et al. (2020). 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. (2020) and observed in, for example, Bauerheim et al. (2014) and Sogaro et al. (2019). This phenomenon is known as mode veering (Seyranian et al. 2005), 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 2010) and the MICCA annular combustor (Bourgouin et al. 2013). 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. (2020) 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 (2010). 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.

Figure 1. Modeshapes of the lowest frequency acoustic mode for $\tau =4\ \textrm {ms}$ and different values of $n$. Modeshapes A and B appear for vanishing values of $n$, and correspond to an acoustic and an ITA mode, respectively. Thermoacosutic modes for finite values of $n$ will generally inherit features from both acoustic and ITA modes. This is particularly evident when the eigenvalue is close to an EP, as for Modeshape C shown here.

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. 2012). 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. (2017), 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. (2012) 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. 2017). 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. (2017) and Orchini et al. (2020). 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 2015).

Figure 2. $(a)$ Eigenvalue trajectories estimated with perturbation theory at 1st (dashed black line) and 30th (solid black line) order, compared with exact solutions (thick shaded line). Within the radius of convergence (black markers), comparison with high-order perturbation theory is excellent. $(b)$ Convergence of the estimated radius of convergence (and therefore location of the closest EP) with the order of the perturbation expansion, from (3.3).

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. 2005), 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 2a) 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 2a), 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.

Figure 3. $(a)$ Trajectories of eigenvalues in the vicinity of the EP. The closer the trajectory is to the EP, the stronger is the veering. The thicker lines are those on which $n$ is real-valued, thus physically realizable. $(b)$ The distance of the EP (red star) from the chosen expansion point defines the radius of convergence of the power series. The convergence region in frequency space, between the two black lines in $(a)$, is not small, which means that the perturbation method is robust.

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. 2018; Orchini et al. 2020). 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. (2020b) 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 2003; Noiray & Schuermans 2013; Bothien, Noiray & Schuermans 2015), 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 2011; Ghirardo, Juniper & Moeck 2016; Laera et al. 2017). All these types of oscillations have been observed experimentally (Noiray et al. 2011; Worth & Dawson 2013; Bourgouin et al. 2014b). An annular combustor in which azimuthal instabilities have been investigated in detail is known as the MICCA combustor (Bourgouin et al. 2014b; Prieur et al. 2017). 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. 2014a; Orchini, Mensah & Moeck 2018; Moeck et al. 2019; Yang, Laera & Morgans 2019). 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. (2019). 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. (2017), 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 1878; Dowling 1995), which can be achieved, for example, by means of fuel staging or the use of cylindrical burner outlets (Krebs et al. 2002; Noiray et al. 2011).

Figure 4. MICCA combustor model employed in this study. The geometrical details and speed of sound field are the same as in Mensah et al. (2019). The mesh maintains the reflection symmetry of each burner segment and the discrete rotational symmetry of the whole geometry; it contains approximately 60 000 tetrahedra.

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.

Figure 5. Radius of convergence of the perturbed eigenvalues estimated from the expansion coefficients (markers) and obtained by identifying the closest branch point (shaded lines). The perturbation patterns are shown as insets. $(a)$ Perturbation pattern I retains the symmetry and the eigenvalue remains degenerate regardless of the perturbation strength. $(b)$ Perturbation pattern II breaks the geometrical symmetries, and each eigenvalue branch exhibits a different convergence behaviour.

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.4a,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.2a) 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.

Figure 6. Comparison between the growth rate and frequency of the eigenvalues close to a degenerate state ($\tau _0=3\ \textrm {ms}$) for the symmetry breaking perturbation pattern II. The coloured, shaded regions indicate the radius of convergence of the respective branch. The exact solutions (solid lines) are poorly approximated by the first-order sensitivity (dashed lines) but well approximated by high-order power series expansions (markers).

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. (2020), 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. (2020).

Figure 7. Cut view of the absolute value of the pressure eigenvector in the plenum of the MICCA combustor when the non-symmetric flame staging pattern II is considered. The degeneracy is resolved, and the two modeshapes that were spanning the twofold-degenerate subspace in the unperturbed case are now well defined. The symmetries in the system's solution are lost, as correctly captured by perturbation theory. (a,d) Eigenvalues and eigenvectors obtained by solving the non-symmetric configuration directly with the Helmholtz solver. (b,e) Eigenvalues and eigenvectors obtained from perturbation theory at 10th order, applied to the nominally symmetric, degenerate case. (c, f) Error between the exact and approximated eigenvectors.

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 2012; Buschmann et al. 2020b). 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):

(5.1)\begin{equation} \left\langle\left.\hat{\tilde{\boldsymbol{p}}}^{\dagger}_{{def}}\right|\boldsymbol{\mathsf{L}}_{1,0} \hat{\tilde{\boldsymbol{p}}}_{def}\right\rangle = 0. \end{equation}

This is known as self-orthogonality in non-Hermitian quantum mechanics (Moiseyev 2011; Heiss 2012) 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 1990; Lancaster et al. 2004):

(5.2)\begin{equation} s(\varepsilon) \approx s_0 + \sum_{j=1}^N s_j \varepsilon^{\,j/a}. \end{equation}

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

(5.3)\begin{equation} s_1 = {\pm} \sqrt{-\frac{\left\langle\left.\hat{\boldsymbol{p}}^{\dagger}_{{def}}\right|\boldsymbol{\mathsf{L}}_{0,1} \hat{\boldsymbol{p}}_{def}\right\rangle}{\left\langle\left.\hat{\boldsymbol{p}}^{\dagger}_{{def}}\right|\boldsymbol{\mathsf{L}}_{1,0} \hat{\tilde{\boldsymbol{p}}}_{gen}\right\rangle+\left\langle\left.\hat{\boldsymbol{p}}^{\dagger}_{{def}}\right|\boldsymbol{\mathsf{L}}_{2,0} \hat{\boldsymbol{p}}_{def}\right\rangle}}, \end{equation}

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.2a) – which converges to the closest EP in the complex-parameter space – and the method outlined in Orchini et al. (2020) – 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 2018), expressed in non-dimensional units, for which the scalar operator $L$ is given by

(5.4)\begin{equation} L(s) \equiv s^2 + 2 {\rm \pi}\beta \textrm{e}^{-s\tau} + {\rm \pi}^2. \end{equation}

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

(5.5a,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. (2005), 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.5a,b) effectively identify EPs.

From (5.4) and (5.5a,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

(5.6)\begin{equation} s_{{def}} = -\frac{1}{\tau} \pm \mathrm{i} \frac{\sqrt{\left({\rm \pi}\tau\right)^2-1}}{\tau}. \end{equation}

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:

(5.7)\begin{equation} s = s_{{def}} \pm \sqrt{\frac{2{\rm \pi}\beta s_{{def}} \textrm{e}^{-s_{{def}}\tau}}{1 + {\rm \pi}\beta\tau^2 \textrm{e}^{-s_{{def}}\tau}}} (\Delta\tau)^{1/2} +O\left(\Delta\tau\right). \end{equation}

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.

Figure 8. $(a)$ Trajectories of two eigenvalues coalescing at an EP (solid lines), and first-order (dashed) and second-order (dotted) Puiseux approximations of these trajectories. The discontinuous branch point behaviour at the EP is correctly captured by the Puiseux approximations. $(b)$ Relative error between the exact eigenvalues and their approximations by Puiseux series, at various orders in log–log scale. The error scales as $(\Delta \tau )^{1/2}$, which is consistent at an EP with algebraic multiplicity 2.

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.