Localization of higher order exceptional points from finite element model and their applications to duct acoustics.

Benoit Nennig 1

Institut supérieur de mécanique de Paris (SUPMECA), Laboratoire Quartz EA 7393, 3 rue Fernand Hainaut, 93407 Saint-Ouen, France.

Emmanuel Perrey-Debain 2

Sorbonne universités, Université de Technologie de Compiègne, Laboratoire Roberval CNRS FRE 2012, CS 60319, 60203 Compiègne cedex, France.

Martin Ghienne 3

Institut supérieur de mécanique de Paris (SUPMECA), Laboratoire Quartz EA 7393, 3 rue Fernand Hainaut, 93407 Saint-Ouen, France.

ABSTRACT This work reviews the state of the art for high order perturbation method for parametric eigenvalue problems and propose some extensions for the multiparameters case. This approach allows to locate high order exceptional points (EP) arising in eigenvalue problems. EP correspond to a particular tuning of some complex-valued parameters which render the problem degenerate. These non-Hermitian degeneracies have raised considerable attention in the scientific community as these can have a great impact in a variety of physical problems (PT-symmetry, thermo-acoustic or fluid- structure instability, etc.) and their numerical solution. For applications dealing with dissipative acoustic waveguides, strong modal attenuation can be achieved close to EP and a maximum of attenuation occurs at EP of high order corresponding to the coalescence of more than two modes. The method is based on the automatic computation of the successive derivatives of some selected eigenpairs with respect to the parameters so that, after recombination, regular functions can be constructed. This algebraic manipulations permit to build a reduced order model allowing i) to quickly solve the eigenvalue problem for other parameters values, ii) to follow modal branches, iii) to locate higher order EPs. The method is applied to the particular case of a circular duct with a locally reacting liner at its wall which admittance varies with azimuthal position. This particular application allows to gradually test the proposed method with EP of increasing order.

1 benoit.nennig@isea-supmeca.fr

2 emmanuel.perrey-debain@utc.fr

3 martin.ghienne@isae-supmeca.fr

ee ante. noize 21-24 AUGUST SCOTTISH EVENT CAMPUS O ? ? GLASGOW

1. INTRODUCTION

Acoustic treatments are often used to attenuate sound propagation in guiding structures generally encountered in HVAC, exhaust devices, and aircraft engines. In many applications, the walls of the waveguide are acoustically treated with absorbing materials, and it is usual for the analysis to distinguish between two classes of liners depending whether the latter is locally like perforated plate bounded with honeycomb or non-locally reacting like porous materials. In both cases, the propagation of sound waves is best described in terms of duct acoustic modes and the existence of exceptional point (EP) is known to be related to very strong modal attenuation [1–5]. At the EP, eigenvalues and associated eigenvectors coalesce at a branch point singularity in the parametric space. Moreover, in the vicinity of the EP, perturbation analysis show that eigenvalues take the form of a Puiseux series [6] i.e. a series in fractional power of the parameter. From practical point of view, this strong attenuation can be explained by a transmission dominated by surface waves [2,3]. Eigenvalue problems are ubiquous in all areas of physics, parametric eigenvalue problems have thus been the subject of intensive research. In acoustic waveguides, the usual parameters can be the frequency or the wall impedance. Most approaches are based on first or second order perturbation methods [7]. Recently, high order perturbation methods have been proposed to optimize the numerical e ffi ciency when the parameter appears explicitly in the numerical model. In [8, 9] the bordered matrix [10] is used to compute eigenvalue derivatives and the adjoint vector is used in [11, 12]. It is noteworthy that the bordered matrix approach is limited (at least in this form) to simple eigenvalue whereas [11] extend the adjoint method to semi-simple eigenvalue, ie a repeated eigenvalue with di ff erent eigenvectors. However, high order Taylor expansion of the eigenvalue is still limited by the presence of branch point singularity. In [8], it has been proposed to recombine the Taylor expansion of the eigenvalue to build analytic functions which are not a ff ected by the singularity at the EP. This algebraic manipulation allows to extend the radius of convergence of the eigenvalue approximation [9], the localization of EP using standard root-finding algorithms and the computation of the associated Puiseux series up to an arbitrary order. EP of order N (EP N ) require more parameters and their localization is more challenging than EP2. In addition, the eigenvalue sensitivity increases with N and the eigenvalue estimation could strongly depend on the rounding error and split [13]. Despite this, practical realization have been proposed [14], although, to the best authors knowledge, no generic algorithm exists to find EP N . The main strategy to obtained high order EP is generally to work on simple (or reduced) system and to combine them [15]. The aim of this work is to propose i) a new perturbation method for eigenvalue problem and ii) a generic algorithm to locate EP N . For this purpose, a new framework, based on analytic functions combining more than two eigenvalues to improve the radius of convergence and to take into account several parameters, is proposed. This can be done through the new concept of partial characteristic polynomial (PCP), which consists in reconstructing a part of the characteristic polynomial 4 from the Taylor expansion of some selected eigenpairs obtained at an initial value of the parameters. To compute the successive derivatives, the bordered matrix [10] is used here as in [8] (but other approaches can be used). For new parameters values, the eigenvalues can be obtained quickly as the roots of a low order polynomial. The EP of order N can then be found by imposing that the PCP and its N − 1 successive derivatives with respect to the eigenvalue vanish. The main ingredient of the method are given in sec. 2, applications to a 2D and 3D waveguides illustrating the accuracy and the interest of the proposed method is presented in sec. 3 and sec. 4.

4 The characteristic polynomial cannot be computed for big matrices obtained after standard discretization methods like FEM.

2. PARTIAL CHARACTERISTIC POLYNOMIAL

2.1. General presentation Let us consider a pair of eigenvalue, denoted by λ − ( ν ) and λ + ( ν ), depending on a scalar parameter ν . For the applications we have in mind, the ν parameter is tipically a wall admittance. In [9], it has been shown that the auxiliary functions

g ( ν ) = λ + + λ − , (1)

h ( ν ) = ( λ + − λ − ) 2 , (2)

are analytic even if the two eigenvalues present a branch point singularity and that eigenvalues’s Taylor series (denoted by T • ) is not convergent. In the same way, the product of the eigenvalues is also an analytic function of ν . It follows that analytic auxiliary functions (AAF) scheme,

2 ≈ T g ( ν ) ± √ T h ( ν )

λ ± ( ν ) = g ( ν ) ± √ h ( ν )

2 , (3)

behaves generally better than Puiseux or Taylor series to reconstruct eigenvalue loci as a function of ν . In the vicinity of its EP, the Puiseux series provides the good branch structure of the eigenvalue Riemann surfaces. Nevertheless, the implicit function theorem shows that Puiseux series radius of convergence is limited by the presence of the next EP even if this EP involved the same pair of eigenvalue, while in this specific case AAF scheme remains analytic and is able to accurately reconstruct the eigenvalues (but not directly the branches). Getting the eigenvalue from the analytic auxiliary function is equivalent to find the roots of the polynomial

P ( λ ) = ( λ − λ + )( λ − λ − )

= λ 2 − ( λ + + λ − ) λ + λ + λ − , (4)

whom coe ffi cients are analytic functions of the parameter ν (the explicit dependency of ν has been dropped for clarity). In particular, in Equation 3, the link between h and the polynomial discriminant is obvious. This approach can be generalized with the concept of partial characteristic polynomial and to vectorial parameters ν . For an analytic matrix function L ( λ ( ν ) , ν ), its characteristic polynomial ˆ P ( λ, ν ) is also an analytic function. Assuming that ˆ P ( λ, ν ) can be split into

ˆ P ( λ ( ν )) = λ − λ 1 ( ν )
. . . λ − λ ℓ ( ν )
. . . λ − λ L ( ν )
· f ( λ ( ν )) (5)

where the λ ℓ stand for all the eigenvalues (counting their multiplicities) contained in a disk in the λ -complex plane and f is an analytic function. We call this set Λ and its cardinal is L = | Λ | . We call

P ( λ, ν ) = λ − λ 1 ( ν )
. . . λ − λ ℓ ( ν )
. . . λ − λ L ( ν )
(6)

the partial characteristic polynomial of L ( λ ( ν ) , ν ). Using the Vieta’s formulas, the coe ffi cient of this polynomial can be expressed as analytic function of the eigenvalue.

P ( λ, ν ) = a L ( ν ) λ L + a L − 1 ( ν ) λ L − 1 + · · · + a 0 ( ν ) , (7)

where

a k = X

Y

i ∈ c λ i ( ν ) , (8)

c ∈ C k

and the set C k = { all ways to choose ( k −| Λ | ) di ff erent eigenvalues in Λ } . For instance, for a quintic polynomial build using a set of 5 eigenvalues, we get

a 5 = 1 , (9a) a 4 = − λ 0 − λ 1 − λ 2 − λ 3 − λ 4 (9b) a 3 = λ 0 λ 1 + λ 0 λ 2 + λ 0 λ 3 + λ 0 λ 4 + λ 1 λ 2 + λ 1 λ 3 (9c) + λ 1 λ 4 + λ 2 λ 3 + λ 2 λ 4 + λ 3 λ 4 , a 2 = − λ 0 λ 1 λ 2 − λ 0 λ 1 λ 3 − λ 0 λ 1 λ 4 − λ 0 λ 2 λ 3 (9d) − λ 0 λ 2 λ 4 − λ 0 λ 3 λ 4 − λ 1 λ 2 λ 3 − λ 1 λ 2 λ 4 − λ 1 λ 3 λ 4 − λ 2 λ 3 λ 4 , a 1 = λ 0 λ 1 λ 2 λ 3 + λ 0 λ 1 λ 2 λ 4 + λ 0 λ 1 λ 3 λ 4 + λ 0 λ 2 λ 3 λ 4 (9e) + λ 1 λ 2 λ 3 λ 4 , a 0 = − λ 0 λ 1 λ 2 λ 3 λ 4 . (9f)

Each polynomial coe ffi cient a k ( k = 0 , . . . , L − 1) is analytic, their Taylor series can thus be obtained thanks to multivariate Leibniz’ rule applied to Equation 8 from the Taylor series of each eigenvalue. The final step is to find the roots of the partial polynomial P ( λ, ν ) for a given range of ν to recover the λ i ( ν ) ( i = 0 , . . . , L − 1). Depending on the number of eigenvalue in Λ , the roots can be obtained numerically or analytically if the polynomial degree is lower than 4. Practicaly, the numerical burden due to the numerical root-finding is negligible when compared to eigenvalue computation of matrices obtained after standard discretization methods like finite element method. This representation also allows to locate EP of order two or higher.

2.2. Higher order EP location EP corresponds generically to multiple roots of the PCP for complex symmetric problems [16]. If the first derivative ∂ λ P ( λ ; ν ) vanishes, it corresponds to an EP2. If the second derivatives also vanishes it corresponds to an EP3 and so one. It is noteworthy that PCP derivative with respect to the eigenvalue also vanishes for semi-simple eigenvalue. Further checks could be required to distinguish EP from semi-simple, for instance by inspecting Puiseux series coe ffi cients [17]. To locate the finite set of EP3, we need to solve the multivariate polynomial system 

P ( λ ; ν ) = 0 , ∂ λ P ( λ ; ν ) = 0 ,

∂ 2 λ P ( λ ; ν ) = 0 , (10)

with ν ∈ C 2 . The Bezout theorem (see [18] for instance) ensures that this system has only discrete solutions. The system can be solved using

– Newton solver . Since such approach need an initial guess, a mesh of the parametric space is required and the solver need to be called for all starting points. This approach is generally fast but may miss some roots.

– Homotopy solver [18,19]. Methods of that kind are able to find all the solution of Equation 10 numerically. The basic idea is to exploit analyticity to move slowly from the solution of a trivial similar problem to the solution of the target problem.

Both approaches are e ffi cient, but due to the truncated Taylor series, the number of spurious solutions may be huge. An estimation can be obtained from the Bezout Number [18]. To filter the spurious roots, comparing solutions for two truncation orders in the Taylor approximation of a k is generally su ffi cient. It could also be combined with an estimation of the radius of convergence of the PCP terms.

Figure 1: Lined duct with two admittances.

3. 2D DUCT APPLICATION

3.1. Problem statement We consider a two-dimensional acoustic waveguide of infinite length and unit width, lined with two admittances µ and ν , as described in Fig. 1. This reference problem have been already studied in [5] and can lead to a finite set of EP3 or a continuum of EP2. Their values are useful to validate the proposed approach. In the duct, the acoustic pressure satisfies the Helmholtz equation (e − i ω t

convention is adopted here) ∆ p + k 2 a p = 0 , (11)

where k a = ω/ c a is the wavenumber, c a is the sound speed and ω is the angular frequency. On both walls, the liner is assumed to be locally reacting which implies that the pressure must satisfy the Robin boundary condition ∂ y p = − µ p , at y = 0 and, ∂ y p = ν p , at y = 1 (12)

Using invariance along the waveguide x -axis, the modal analysis is performed by assuming that the pressure field can be written in the separable form p = φ ( y )e i β x . Here, function φ ( y ) is the mode shape and β the axial wavenumber. The weak formulation associated with the Helmholtz equation gives

− Z 1

0 ∂ y ψ · ∂ y φ d y + ( k 2 a − β 2 ) Z 1

0 ψφ d y + µψ (0) φ (0) + νψ (1) φ (1) = 0 , (13)

where ψ stands for the test function. Once the variational formulation is discretized with linear Lagrangian finite element, we obtain a generalized eigenvalue problem of the form

L λ ( ν ) , ν
x ( ν ) = − K + ( k 2 a − λ ) M + ν Γ 1 + µ Γ 2
x ( ν ) = 0 . (14)

To be consistent with our notation, we put λ = β 2 and ν = ( µ, ν ). Here, the vector x contains the finite element nodal values of the acoustic pressure, the matrix Γ i ( i = 1 , 2) comes from the admittance boundary condition and K and M are the standard sti ff ness and mass matrices respectively. We can easily calculate formally the partial derivative of Equation 14 with respect to ν components and λ required to get the eigenvalue derivative as described in [8].

3.2. Results This approach is implemented in EasterEig [20], an open source framework dedicated to perturbation of eigenvalue problem. We use 200 linear Lagrange element in the cross section of the duct. To illustrate the numerical e ffi ciency of eigenvalue reconstruction based on PCP, the absolute error E between direct and PCP computation is given in Fig. 2 for increasing values of the perturbation parameter δ such ν = ν 0 + e 0 . 3i δ . The calculation are performed by combining the Taylor series of 2, 4, 6, 8 or 10 eigenvalues (sorted by modulus) computed at ν 0 = (7 . 01265 − 4 . 76715i , 2 . 89872 − 2 . 47i). It can be observed that increasing the number of eigenvalues in the PCP strongly enhances the

5

0

log 10 E

− 5

T λ N=2 N=4 N=6 N=8 N=10

− 10

− 15

0 2 4 6 8 10 | δ |

Figure 2: Evolution of error on eigenvalue reconstruction according to the number of eigenvalues in the PCP and the modulus of the perturbation parameter δ .

reconstruction quality although the error slightly increases when δ = 0 due to recombination round- o ff . This improvement can be explained by a mutual regularization of the branch point singularity in eigenvalue Taylor series. For this 1D problem, the prediction is valid up to 2 digits for perturbation parameter equal to 10, spanning most of interesting admittance values with a single eigenvalue direct computation. For comparison, the direct Taylor expansion of the eigenvalue T λ n has a smaller validity range and is limited to δ < 0 . 8. The second aspect of the proposed approach concerns the EP location. Again, we start from ν 0 = (7 . 01265 − 4 . 76715i , 2 . 89872 − 2 . 47i) and combine the 6 first eigenvalues. The homotopy method is used to solve Equation 10. The solution is a set containing the triplets ( λ , µ , ν ) and contains many spurious roots due to the Taylor series truncation. To filter-out these spurious roots we exploit the fact that spurious roots are extremely sensitive and we solve the system with two di ff erent truncation order with respect to µ an ν (here four and five). The value obtained for ν that yields to EP3 are presented in Fig. 3. Note that because of the symmetry of the problem, the role of ν and µ is interchangeable thus the solution are found by pairs. It can be shown that several EP 3 are accurately identified among the spurious roots when two truncation orders are compared.

4. CIRCULAR DUCT

It is well known that EP2 can be obtained from circular ducts treated with a homogeneous liners [1, 21, 22]. The degeneracy in this case concerns radial modes at fixed azimuthal order. In order to get higher order EP, more tuning parameters are required. Here, we consider a 3D duct with a circular cross section containing two lined parts separated by a rigid walls, as shown in Fig. 4. The cross-section is discretized with quadratic Lagrange finite element and yield 1100 degree of freedom. We consider only the no flow case. This configuration is geometrically rather similar to the 2D duct treated in the previous section, but in higher dimension the eigenvalues structure drastically di ff er from 1D uniform case since there are several family of modes like whispering gallery, azimuthal, radial and sometimes trapped modes [23]. The algorithm found several EP3 values ( ν 0 = (7 . 01 − 4 . 76i , 2 . 89 − 2 . 47i)) and one more iteration is done to improve accuracy. The eigenvalues (Fig. 5) and the eigenvector (Fig. 6) are presented for µ = 3 . 3189 + 3 . 1685i and ν = 3 . 1765 + 4 . 1483i at k = 1 and the numerical value for the axial wavenumbers are given in Fig. 6. For this specific case, the EP3 is not the least attenuated mode in the system (see Fig. 6) with a first mode with β 0 = 0 . 944 + 3 . 696i. Work is on going to better understand the e ff ect of the rigid splitters on the merging of modes.

12

10

8

Im ν

6

4

Homotopy (5, 5)

2

Homotopy (4, 4) ref. ν 0

0

0 1 2 3 4 5 Re ν

Figure 3: EP3 solutions obtained from the initial guess ν 0 = (7 . 01265 − 4 . 76715i , 2 . 89872 − 2 . 47i)) and combining the 6 first eigenvalues.

20

18

16

14

12

Im

10

8

6

4

0.0 0.5 1.0 1.5 2.0 2.5 Re

Figure 4: Circular duct cross-section of radius r = 0 . 5.

Figure 5: Axial wavenumber obtained at an EP3 with µ = 3 . 3189 + 3 . 1685i and ν = 3 . 1765 + 4 . 1483i at k = 1.

5. CONCLUSION

A new eigenvalue perturbation method has been proposed. The method is based on the partial characteristic polynomial concept. This approach allows to improve the radius of convergence of eigenvalue reconstruction for parametric problems depending on several parameters and to locate higher order EP. This generic approach is applied on circular ducts. For instance EP3 have been found for wall presenting two lined parts separated by rigid wall portions. Work is on going to better understand the merging mechanism and the relation with strong attenuation as observed in 2D ducts [5]. The method can also be applied to other parametric eigenvalue problems like those arising in thermoacoustic [11,12], in topological acoustic waveguides [24] or in structural dynamics.

0.944+3.696j

2.852+4.001j

2.856+4.003j

0.0 0.4 0.8 1.1 |p| mode #0

0.8 1.9 3.0 4.1 |p| mode #1

0.8 1.9 3.0 4.1 |p| mode #2

2.852+4.005j

1.913+5.816j

0.216+7.249j

0.8 1.9 3.0 4.1 |p| mode #3

0.0 0.5 1.0 1.5 |p| mode #4

0.0 0.8 1.6 2.4 |p| mode #5

Figure 6: Modulus of the pressure field for the 6 first eigenvectors obtained at an EP3 with µ = 3 . 3189 + 3 . 1685i and ν = 3 . 1765 + 4 . 1483i. The indicated values correspond to the axial wavenumbers β .

REFERENCES

[1] B. J. Tester. The optimization of modal sound attenuation in duct, in the absence of mean flow. J. Sound Vib. , 27:477–513, January 1973. [2] W. Bi and W. Pagneux. New insights into mode behaviours in waveguides with impedance boundary conditions. arXiv:1511.05508 , 2015. [3] L. Xiong, B. Nennig, Y. Aurégan, and W. Bi. Sound attenuation optimization using metaporous materials tuned on exceptional points. J. Acoust. Soc. Am. , 142(4):2288 – 2297, 2017. [4] X. Qiu, L. Du, X. Jing, X. Sun, M. Åbom, and H. Bodén. Optimality analysis of bulk-reacting liners based on mode-merging design method. J. Sound Vib. , 485:115581, 2020. [5] E. Perrey-Debain, B. Nennig, and J. B. Lawrie. Mode coalescence and the green’s function in a two-dimensional waveguide with arbitrary admittance boundary conditions. Journal of Sound and Vibration , page 116510, 2021. [6] T. Kato. Perturbation Theory for Linear Operators, 2nd edition , page 623pp. Springer-Verlag, Berlin, Heidelberg, 1980. [7] R. M. Lin, J. E. Mottershead, and T. Y. Ng. A state-of-the-art review on theory and engineering applications of eigenvalue and eigenvector derivatives. Mechanical Systems and Signal Processing , 138:106536, 2020. [8] B. Nennig and E. Perrey-Debain. A high order continuation method to locate exceptional points and to compute puiseux series with applications to acoustic waveguides. J. Comp. Phys. , page 109425, 2020. [9] M. Ghienne and B. Nennig. Beyond the limitations of perturbation methods for real random eigenvalue problems using exceptional points and analytic continuation. J. Sound Vib. , page

115398, 2020. [10] A. L. Andrew, K.-W. E. Chu, and P. Lancaster. Derivatives of eigenvalues and eigenvectors of matrix functions. SIAM J. Matrix Anal. Appl. , 14(4):903–926, 1993. [11] A. Orchini, G. A. Mensah, and J. P. Moeck. Perturbation theory of nonlinear, non-self-adjoint eigenvalue problems: Semisimple eigenvalues. Journal of Sound and Vibration , 507:116150, 2021. [12] G. A. Mensah, A. Orchini, and J. P. Moeck. Perturbation theory of nonlinear, non-self-adjoint eigenvalue problems: simple eigenvalues. Journal of Sound and Vibration , 473:115200, 2020. [13] J.-W. Ryu, J.-H. Han, and C.-H. Yi. Classification of multiple arbitrary-order non-hermitian singularities. arXiv preprint arXiv:2112.02547 , 2021. [14] K. Ding, G. Ma, M. Xiao, Z. Q. Zhang, and C. T. Chan. Emergence, coalescence, and topological properties of multiple exceptional points and their experimental realization. Physical Review X , 6(2):021007, 2016. [15] Q. Zhong, J. Kou, ¸SK. Özdemir, and R. El-Ganainy. Hierarchical construction of higher-order exceptional points. Physical review letters , 125(20):203602, 2020. [16] A. P. Seyranian, O. N. Kirillov, and A. A. Mailybaev. Coupling of eigenvalues of complex matrices at diabolic and exceptional points. J. Phys. A , 38(8):1723, 2005. [17] A. Welters. On explicit recursive formulas in the spectral perturbation analysis of a Jordan block. SIAM J. Matrix Anal. Appl. , 32(1):1–22, 2011. [18] S. M Wise, A. J. Sommese, and L. T. Watson. Algorithm 801: Polsys_plp: A partitioned linear product homotopy code for solving polynomial systems of equations. ACM Transactions on Mathematical Software (TOMS) , 26(1):176–200, 2000. [19] B. Nennig. A python wrapper to the fortran package polsys_plp that solve polynomial systems with homotopy method. https://github.com/nennigb/pypolsys , 2020. [20] B. Nennig and M. Ghienne. A library to locate exceptional points and to reconstruct eigenvalues loci. https://github.com/nennigb/EasterEig , 2020. [21] W. Koch. Attenuation of sound in multi-element acoustically lined rectangular ducts in the absence of mean flow. J. Sound Vib. , 52(4):459–496, 1977. [22] W. E. Zorumski and J. P. Mason. Multiple eigenvalues of sound-absorbing circular and annular ducts. J. Acoust. Soc. Am. , 55(6):1158–1165, 1974. [23] E. J. Brambley, A. M. J. Davis, and N. Peake. Eigenmodes of lined flow ducts with rigid splices. Journal of Fluid Mechanics , 690:399–425, 2012. [24] A. Coutant, V. Achilleos, O. Richoux, G. Theocharis, and V. Pagneux. Subwavelength Su- Schrie ff er-Heeger topological modes in acoustic waveguides. arXiv preprint arXiv:2112.13347 , 2021.