The modal analysis of wave problems of unbounded type involves a continuous sum of radiation modes. This continuum is difficult to handle mathematically and physically. It can be approximated by a discrete set of leaky modes, corresponding to improper modes growing to infinity. Perfectly matched layers (PMLs) have been widely applied in numerical methods to efficiently simulate infinite media, most often without considering a modal approach. This letter aims to bring insight into the modal basis computed with PMLs. PMLs actually enable to reveal of the contribution of leaky modes by redefining the continua (two for elastodynamics), discretized after PML truncation.

## 1. Introduction

Modal analysis is a powerful tool for treating bounded wave problems, for which modes provide a discrete basis that allows in-depth physical analysis. For unbounded problems, the modal basis is yet more complicated because it involves a continuous sum of modes (radiation modes) in addition to a discrete set (trapped modes, if any).^{1} This contribution of the continuum is usually quite difficult to calculate and interpret from a physical point of view. Based on steepest descent techniques, for example, it can be approximated by a discrete set of another class of modes, the so-called leaky modes.^{2} These modes behave like attenuated modes and reveal some interesting physical features. For instance, they can give the leakage attenuation of wave packets in open waveguides, which is of particular interest in the field of nondestructive evaluation (NDE).^{3} However, leaky modes correspond to improper modes that spatially grow to infinity in the unbounded direction of the problem (they are not part of the modal basis). Their use in modal analysis hence requires special care.

As for numerical methods, perfectly matched layers (PMLs) have been widely employed to efficiently simulate infinite media. More particularly, they have been successfully applied to compute modes in open systems.^{4,5} It has been shown that PMLs allow getting both trapped and leaky modes as well as a third type of mode sometimes called PML modes.^{6} PML modes mainly resonate inside the PML region (their fields strongly depend on PML parameters). However, the contribution of these different kinds of modes in a modal analysis may still be unclear. This letter aims to bring insight into the modal basis computed with PMLs and to justify its use for the analysis of excitation and scattering problems using a modal approach. Without loss of generality, two-dimensional waveguides are considered. For simplicity, the analysis starts from the scalar wave equation and is then extended to elastodynamics, where two continua occur instead of one.

## 2. Background

This section briefly recalls some theoretical results^{1,2,7} taking the example of a two-dimensional bilayer waveguide, homogeneous and infinite in the axial direction *z* as shown in Fig. 1. In the transverse direction, the bottom layer is labelled 0 and extends from *x* = −*h* to *x* = 0. The top layer is infinite and labelled $\u221e$. The wave field solution to the problem satisfies the Helmholtz equation in both layers and is denoted in the angular frequency-wavenumber domain by $\psi (x,\beta ,\omega )$, where *β* and *ω* are the axial wavenumber and the angular frequency. A convention in $ej(\beta z\u2212\omega t)$ is adopted for the harmonic dependence of wave fields. The transverse wavenumbers in the two layers are given by $\alpha 0,\u221e2$ = $k0,\u221e2\u2212\beta 2$, with $k0,\u221e$ = $\omega /c0,\u221e$. In this paper, ′ or $Re$ and ″ or $Im$ will denote, respectively, the real part and the imaginary part of a complex number. As an example, let us assume that the waveguide in Fig. 1 is excited by an out-of-plane line source located at $(x,z)=(\u2212h,0)$ through the Neumann boundary condition $d\psi (\u2212h,z,\omega )/dx=s(\omega )\delta (z)$.

The solution must also satisfy continuity conditions at the interface (*x* = 0) and must be bounded and outwards at infinity^{2} ($x\u2192\u221e$). The spatial inverse Fourier transform of the solution leads to an expression of the following form:

where *C _{m}* is a ratio of material properties which depends on the problem type (i.e., dielectric, acoustic, or SH-elastic). First, a lossy outer medium is assumed (i.e., $k\u221e\u2033$ > 0). The key point is the value of $\alpha \u221e$, which depends on the sign of the square root of $k\u221e2\u2212\beta 2$ and hence defines a two-sheeted Riemann surface.

^{2}Let us consider the two sheets defined by $\alpha \u221e\u2033$ > 0 and $\alpha \u221e\u2033$ < 0, respectively. The corresponding branch cut is given by $\alpha \u221e\u2033$ = 0 and is a portion of hyperbola

^{1}in the ($\beta \u2032,\u2009\beta \u2033$) plane beginning at branch point ($k\u221e\u2032,\u2009k\u221e\u2033$), as shown in Fig. 2(a). The remaining portion of the hyperbola satisfies the equation $\alpha \u221e\u2032$ = 0.

The integrals in Eq. (1) can be evaluated thanks to Cauchy's residue theorem. The path $\mathbb{R}$ is conveniently closed by a semi-infinite circle along which the integrals vanish. This implies that the exponentials $ej\alpha \u221ex$ and $ej\beta z$ must tend to zero at infinity. For clarity, let us focus on the solution in the region *z* > 0, where waves are positive-going and decreasing ($\beta \u2033$ > 0). In this case, the integration must be carried out on the upper semi-circle $C+$, defined by $\beta \u2033$ > 0 and located on the proper sheet, which is $\alpha \u221e\u2033$ > 0 (the improper sheet being $\alpha \u221e\u2033$ < 0). Furthermore, the branch cut is a discontinuity that cannot be crossed. It must be circumvented by the path $\Gamma +$ [see Fig. 2(a)]. Finally, on the contour $\mathbb{R}+\Gamma ++C+$, the application of the residue theorem leads to $\u222b\mathbb{R}=2\pi j\u2211 residues\u2212\u222b\Gamma +$. The residues are calculated from the poles [zeros of denominators of Eq. (1)] that lie on the proper Riemann sheet $\alpha \u221e\u2033$ > 0, that is to say, the trapped modes. Thus, two parts contribute to $\psi (x,z,\omega )$. The first one is a discrete sum of positive going trapped modes. The second part is a continuous sum and stands for the continuum of radiation modes.

Together, trapped modes and radiation modes form the modal basis of an open waveguide:^{1}

As mentioned previously, the continuum is difficult to handle. Its contribution can be approximated by a discrete sum of a third class of modes, the so-called leaky modes (e.g., using the steepest descent method^{8}). It has to be emphasised that leaky modes are also poles of Eq. (1). However, they do not belong to the modal basis because they are improper poles that lie in the negative Riemann sheet^{2}^{,}$\alpha \u221e\u2033$ < 0. Hence, leaky modes grow to infinity in the transverse direction. They can approximate the wave field only in a strictly delimited spatial region.^{9} From a physical point of view, leaky modes might be seen in the bottom layer as imperfect guided waves that suffer from leakage attenuation, while trapped modes are guided waves that propagate without leakage attenuation. Radiation modes form standing waves in the transverse direction and can be propagative or evanescent in the axial direction.^{7}

The lossless case ($k\u221e\u2033$ = 0), not considered so far, can be readily deduced asymptotically and is shown in Fig. 2(b). If the medium is lossless, trapped modes have purely real wavenumbers^{2} and exist only if $k\u221e$ < *β* < *k*_{0}.

Note that the problem given by Eq. (1) has been introduced in this section for illustrative purposes only. The main results of this letter will apply to more complex open waveguides (e.g., multilayered or excited by other types of source), since the analysis in the remainder essentially lies in the definition of the branch cut of the square root, which is intrinsic to the halfspace.

## 3. Modal basis with a PML

A PML is now applied in the infinite top layer. This amounts to perform an analytical extension^{10} for *x* > 0 with the change of variable $x\u0303$ = $\u222b0x\gamma (\zeta )d\zeta $, where $\gamma (x)$ is a complex function such that $arg\u2009\gamma \u2208[0;\pi /2[$. For simplicity, let us consider a constant function inside the whole top layer [$\gamma (x)$ = *γ*]. In this case, $x\u0303$ = *γx* in the top layer. It can be readily checked that solving the Helmholtz equation in both layers yields the same solution as Eq. (1) except that the term $ej\alpha \u221ex$ is now replaced with $ej\alpha \u0303\u221ex$, where $\alpha \u0303\u221e$= $\gamma \alpha \u221e$. The poles of the denominators (trapped and leaky modes) are left unchanged by the PML.

Following the same approach as in Sec. 2, the application of the residue theorem requires considering a new proper Riemann sheet, defined by $\alpha \u0303\u221e\u2033$ > 0. Consequently, the branch cut is now given by $\alpha \u0303\u221e\u2033$ = 0 and defines a new integration path denoted by $\Gamma \u0303+$. To understand the effect of the PML on the modal basis, it is necessary to look further into the branch cut transformation induced by the PML, namely, from $\Gamma +$ to $\Gamma \u0303+$.

The branch cut $\alpha \u0303\u221e\u2033=0=\xb1Im\u2009Re\u2009\alpha \u0303\u221e2+jIm\u2009\alpha \u0303\u221e2$ in the $(\beta \u2032,\u2009\beta \u2033)$ plane can be unequivocally determined by the following conditions: $Im\u2009\alpha \u0303\u221e2=0$ and $Re\u2009\alpha \u0303\u221e2>0$ (see Ref. 1, pp. 180 and 181 for the case without PML). Then, from the expansion $\alpha \u0303\u221e2=(\gamma \u2032+j\gamma \u2033)(k\u221e\u20322\u2212k\u221e\u20332\u2212\beta \u20322+\beta \u20332+2j(k\u221e\u2032k\u221e\u2033\u2212\beta \u2032\beta \u2033))$, the equation $Im\u2009\alpha \u0303\u221e2=2\alpha \u0303\u221e\u2032\alpha \u0303\u221e\u2033=0$ becomes

with *C* = −*A* = $\gamma \u2032\gamma \u2033,\u2009B$ = $\gamma \u20332\u2212\gamma \u20322,\u2009D$ = $\gamma \u2032\gamma \u2033(k\u221e\u20322\u2212k\u221e\u20332)+(\gamma \u20322\u2212\gamma \u20332)k\u221e\u2032k\u221e\u2033$. Equation (3) represents a hyperbola in the frame ($\beta \u2032,\u2009\beta \u2033$). The principal axes of this hyperbola make an angle *θ* with the axes of the frame given by the formula $cot2\theta $ = $(A\u2212C)/B$, which leads to

From the study of the sign of $Re\u2009\alpha \u0303\u221e2$, it can be checked that the branch cut $\alpha \u0303\u221e\u2033$ = 0 is on the upper part of the hyperbola and that the lowest part corresponds to the curve $\alpha \u0303\u221e\u2032$ = 0. Both parts are delimited by the branch point $(k\u221e\u2032,\u2009k\u221e\u2033)$, analogously to the case without PML.

Without PML, $arg\u2009\gamma $ = 0 and *θ* = $\pi /4$: the asymptotes of the hyperbola coincide with the real and imaginary axes [as already shown in Fig. 2(a)]. With PML, the angle becomes $(\pi /4)\u2212arg\u2009\gamma $: the asymptotes are rotated by $\u2212arg\u2009\gamma $. Hence, the branch cut transformation from $\Gamma +$ to $\Gamma \u0303+$ is a rotation of an angle equal to $\u2212arg\u2009\gamma $ [see Fig. 3(a)]. In the lossless case [Fig. 3(b)], note that the PML branch cut is still a hyperbola.

As shown hereafter, this branch cut rotation has an important consequence. With PML, the proper Riemann sheet satisfies $\alpha \u0303\u221e\u2033$ > 0 (i.e., 0 < $arg\u2009\alpha \u0303\u221e$ < *π*) or, equivalently, by using the equality $arg\u2009\alpha \u0303\u221e$ = $arg\u2009\gamma +arg\u2009\alpha \u221e$:

In particular, the PML gives access to the region $\u2212arg\u2009\gamma $ < $arg\u2009\alpha \u221e$ < 0 (shaded region in Fig. 3), which belongs to the initial improper sheet $\alpha \u221e\u2033$ < 0 where leaky modes occur. In other words, the rotation induced by the PML enables to reveal of the contribution of the leaky modes satisfying $arg\u2009\alpha \u221e$ > $\u2212arg\u2009\gamma $. These revealed leaky modes are proper poles of the problem with PML. Note that Eq. (5) coincides with the condition found for a PML to cancel the transverse exponential growth of a leaky mode.^{4,5}

To summarize, the integration contour with PML is $\mathbb{R}$ + $C+$ + $\Gamma \u0303+$ as shown in Fig. 3. Applying Cauchy's residue theorem on this contour yields

The PML hence defines a new modal basis which properly includes some leaky modes. This modal basis still involves a continuum of modes, which mainly resonate inside the PML (PML modes). One must yet keep in mind that the field inside the PML is not physical: Eq. (6) coincides with the physical solution (2) only in the region out of the PML.

Inside the accessible region of the sheet $\alpha \u221e\u2033$ < 0 (i.e., the shaded region in Fig. 3) the sign of $\alpha \u221e\u2032$ is single-valued and can be determined as follows. In the case without PML, the hyperbola given by Eq. (3) delimits two regions in the complex plane, $\alpha \u221e\u2032\alpha \u221e\u2033$ > 0 (hatched region) and $\alpha \u221e\u2032\alpha \u221e\u2033$ < 0 (white region)—see Fig. 2. The former and the latter are hence such that $\alpha \u221e\u2032$ > 0 and $\alpha \u221e\u2032$ < 0 in the proper Riemann sheet (i.e., $\alpha \u221e\u2033$ > 0), and conversely in the improper one ($\alpha \u221e\u2033$ < 0). Consequently, the revealed part of the sheet $\alpha \u221e\u2033$ < 0 is such that $\alpha \u221e\u2032$ > 0. As expected leaky waves are outgoing in the transverse direction.

## 4. The elastodynamic problem

In the elastodynamic case, two transverse wavenumbers must be considered: one for the longitudinal transverse wavenumber, $\alpha l\u221e$, and one for the shear transverse wavenumber, $\alpha s\u221e$. These are defined from $\alpha l\u221e2$ = $kl\u221e2\u2212\beta 2$ and $\alpha s\u221e2$ = $ks\u221e2\u2212\beta 2$, with $kl\u221e$= $\omega /cl\u221e$ and $ks\u221e$ = $\omega /cs\u221e$. As a consequence, two branch cuts now occur (one for each square root).^{11} Following Sec. 3, the branch cuts to consider correspond to zero imaginary parts, $\alpha l\u221e\u2033$ = 0 and $\alpha s\u221e\u2033$ = 0. This yields four Riemann sheets, labelled **1** for ($\alpha l\u221e\u2033$ > 0, $\alpha s\u221e\u2033$ > 0), **2** for ($\alpha l\u221e\u2033$ > 0, $\alpha s\u221e\u2033$ < 0), **3** for ($\alpha l\u221e\u2033$ < 0, $\alpha s\u221e\u2033$ < 0), and **4** for ($\alpha l\u221e\u2033$ < 0, $\alpha s\u221e\u2033$ > 0).

Applying the residue theorem, the elastodynamic solution without PML can then be expressed as a discrete sum of trapped modes and two continuous sums of radiation modes, instead of one in the scalar case. Trapped modes lie on the proper sheet **1** (on the real axis in the lossless case). Extending the results of Sec. 3, the PML rotates both branch cuts by $\u2212arg\u2009\gamma $. As shown in Fig. 4(a) (lossy case), this rotation reveals parts of the initial improper sheets **2**, **3**, and **4**, containing leaky modes. The solution then becomes expressed as a discrete sum of trapped modes, a discrete sum of revealed leaky modes and two continua of PML modes.

Similarly to Sec. 3, the sign of $\alpha l\u221e\u2032$ and $\alpha s\u221e\u2032$ can be determined in the accessible regions of sheets **2**, **3**, and **4**. This yields the following results:

Furthermore, sheet **1** can be decomposed into three regions, as shown in Fig. 4(a). Note that the transverse longitudinal wavenumbers of leaky modes revealed in sheet **2** have a positive imaginary part together with a negative real part. This partial wave hence behaves like a transverse backward wave decaying at infinity. This does not change the nature of the modes in region **2**, which are of leaky type owing to the transverse behavior of the partial shear wave (forward and growing at infinity).

In the lossless case [Fig. 4(b)], the sheet **4** cannot be accessible by the PML. In Fig. 4(a), the revealed region of sheet **4** is located between the two hyperbolas. In Fig. 4(b), these hyperbolas collapse with the axes of the frame, and so does the revealed part of **4**. Nevertheless, it is pointed out that the contribution of this particular sheet is usually neglected with other approaches regardless PMLs.^{12,13}

## 5. Discussion

In addition to trapped modes, a second class of proper modes is likely to occur depending on material properties.^{8} These modes lie on the proper sheet but are complex even in the lossless case, as opposed to trapped ones. Hence, the discrete sum of trapped modes in Eqs. (2) and (6) should also include these complex proper modes (when they exist). Furthermore, trapped modes of backwards type may occur. Such modes have opposite sign of phase and group velocity and have not been represented in the figures of this paper for simplicity (backwards modes lie on the negative real axis while propagating in the positive direction, and conversely).

Only infinite PMLs have been considered here. Yet for numerical implementation, PMLs must be truncated at some finite distance. Both trapped and leaky modes can be accurately computed, even after PML truncation, provided that the PML is properly parametrized.^{4,5} However, the PML truncation transforms the continuous sum of PML modes into a discrete one.^{6} At first sight, the role of this discrete sum in the convergence of solutions might be questionable owing to the fact that discrete PML modes strongly depend on the boundary conditions arbitrarily set at the end of the truncated PML. Nevertheless, it has been proved for homogeneous scalar problems that the continuum of PML modes can be truly approximated by a sum over discrete PML modes.^{14} The convergence of the discrete sum towards the continuous sum can be achieved depending on the thickness of the PML, its attenuation, and the number of PML modes. Such a convergence cannot be necessarily achieved with absorbing layers of non-perfectly matched type.^{9}

To conclude, the PML technique may significantly simplify the modal analysis of unbounded problems: this tool, easy to implement from a computational point of view, enables to naturally reveal of the contribution of leaky modes. Furthermore, the discretized sum of PML modes can be included in the modal basis in order to improve the accuracy of numerical solutions, which may be necessary when the contribution of leaky modes is not sufficient for a desired accuracy. This fully justifies from a numerical point of view the use of the PML technique for the modal analysis of unbounded problems. Numerical results will be left for a separate publication. It must be yet recalled that the PML technique is known to fail in elastodynamics for some specific properties of anisotropy of the PML medium.^{15} Such situations are beyond the scope of this paper, in which isotropic materials have been implicitly assumed.

## Acknowledgments

The authors wish to thank Région Pays de la Loire for their financial support.