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)

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

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

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

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. 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).

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. 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 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. 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 Moeck2020b) 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 Candel2014b). 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 Candel2014b; 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 Candel2014a; 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).

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. (Reference Mensah, Magri, Orchini and Moeck2019). 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. (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).

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 Reference Beyn2012; Buschmann et al. Reference Buschmann, Mensah, Nicoud and Moeck2020b). 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.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

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