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.

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.

This section briefly recalls some theoretical results1,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 . The wave field solution to the problem satisfies the Helmholtz equation in both layers and is denoted in the angular frequency-wavenumber domain by ψ(x,β,ω), where β and ω are the axial wavenumber and the angular frequency. A convention in ej(βzωt) is adopted for the harmonic dependence of wave fields. The transverse wavenumbers in the two layers are given by α0,2 = k0,2β2, with k0, = ω/c0,. 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)=(h,0) through the Neumann boundary condition dψ(h,z,ω)/dx=s(ω)δ(z).

Fig. 1.

Example of a two-dimensional bilayer open waveguide.

Fig. 1.

Example of a two-dimensional bilayer open waveguide.

Close modal

The solution must also satisfy continuity conditions at the interface (x = 0) and must be bounded and outwards at infinity2 (x). The spatial inverse Fourier transform of the solution leads to an expression of the following form:

ψ(x,z,ω)={ s(ω)2πα0cos(α0x)+jCmαsin(α0x)α0[α0sin(α0h)+jCmαcos(α0h)]ejβzdβifhx0,s(ω)2π1α0sin(α0h)+jCmαcos(α0h)ejαxejβzdβifx>0,
(1)

where Cm 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 > 0). The key point is the value of α, which depends on the sign of the square root of k2β2 and hence defines a two-sheeted Riemann surface.2 Let us consider the two sheets defined by α > 0 and α < 0, respectively. The corresponding branch cut is given by α = 0 and is a portion of hyperbola1 in the (β,β) plane beginning at branch point (k,k), as shown in Fig. 2(a). The remaining portion of the hyperbola satisfies the equation α = 0.

Fig. 2.

Representation of the proper sheet α > 0 in the upper complex plane β > 0: (a) lossy medium, (b) lossless. Branch cut α = 0 (thick line), curve α = 0 (dotted line), trapped modes (), and integration contour (dashed line). Hatched region: α > 0. White region: α < 0.

Fig. 2.

Representation of the proper sheet α > 0 in the upper complex plane β > 0: (a) lossy medium, (b) lossless. Branch cut α = 0 (thick line), curve α = 0 (dotted line), trapped modes (), and integration contour (dashed line). Hatched region: α > 0. White region: α < 0.

Close modal

The integrals in Eq. (1) can be evaluated thanks to Cauchy's residue theorem. The path is conveniently closed by a semi-infinite circle along which the integrals vanish. This implies that the exponentials ejαx and ejβ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 (β > 0). In this case, the integration must be carried out on the upper semi-circle C+, defined by β > 0 and located on the proper sheet, which is α > 0 (the improper sheet being α < 0). Furthermore, the branch cut is a discontinuity that cannot be crossed. It must be circumvented by the path Γ+ [see Fig. 2(a)]. Finally, on the contour +Γ++C+, the application of the residue theorem leads to =2πj residuesΓ+. The residues are calculated from the poles [zeros of denominators of Eq. (1)] that lie on the proper Riemann sheet α > 0, that is to say, the trapped modes. Thus, two parts contribute to ψ(x,z,ω). 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 

ψ(x,z,ω)= trapped+Γ+radiationmodes.
(2)

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 method8). 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 sheet2,α < 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 = 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 wavenumbers2 and exist only if k < β < k0.

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.

A PML is now applied in the infinite top layer. This amounts to perform an analytical extension10 for x > 0 with the change of variable x̃ = 0xγ(ζ)dζ, where γ(x) is a complex function such that argγ[0;π/2[. For simplicity, let us consider a constant function inside the whole top layer [γ(x) = γ]. In this case, x̃ = γ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αx is now replaced with ejα̃x, where α̃= γα. 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 α̃ > 0. Consequently, the branch cut is now given by α̃ = 0 and defines a new integration path denoted by Γ̃+. 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 Γ+ to Γ̃+.

The branch cut α̃=0=±ImReα̃2+jImα̃2 in the (β,β) plane can be unequivocally determined by the following conditions: Imα̃2=0 and Reα̃2>0 (see Ref. 1, pp. 180 and 181 for the case without PML). Then, from the expansion α̃2=(γ+jγ)(k2k2β2+β2+2j(kkββ)), the equation Imα̃2=2α̃α̃=0 becomes

Aβ2+Bββ+Cβ2+D=0,
(3)

with C = −A = γγ,B = γ2γ2,D = γγ(k2k2)+(γ2γ2)kk. Equation (3) represents a hyperbola in the frame (β,β). The principal axes of this hyperbola make an angle θ with the axes of the frame given by the formula cot2θ = (AC)/B, which leads to

θ=π4argγ.
(4)

From the study of the sign of Reα̃2, it can be checked that the branch cut α̃ = 0 is on the upper part of the hyperbola and that the lowest part corresponds to the curve α̃ = 0. Both parts are delimited by the branch point (k,k), analogously to the case without PML.

Without PML, argγ = 0 and θ = π/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 (π/4)argγ: the asymptotes are rotated by argγ. Hence, the branch cut transformation from Γ+ to Γ̃+ is a rotation of an angle equal to argγ [see Fig. 3(a)]. In the lossless case [Fig. 3(b)], note that the PML branch cut is still a hyperbola.

Fig. 3.

Representation of the proper sheet with PML α̃ > 0 in the upper complex plane β > 0: (a) lossy medium, (b) lossless. Same legend as Fig. 2 and PML branch cut α̃ = 0 (thick gray line) with its asymptotes (dashed dotted), leaky modes (), revealed part of the improper sheet α < 0 (gray shaded region, satisfying α > 0).

Fig. 3.

Representation of the proper sheet with PML α̃ > 0 in the upper complex plane β > 0: (a) lossy medium, (b) lossless. Same legend as Fig. 2 and PML branch cut α̃ = 0 (thick gray line) with its asymptotes (dashed dotted), leaky modes (), revealed part of the improper sheet α < 0 (gray shaded region, satisfying α > 0).

Close modal

As shown hereafter, this branch cut rotation has an important consequence. With PML, the proper Riemann sheet satisfies α̃ > 0 (i.e., 0 < argα̃ < π) or, equivalently, by using the equality argα̃ = argγ+argα:

argγ<argα<πargγ.
(5)

In particular, the PML gives access to the region argγ < argα < 0 (shaded region in Fig. 3), which belongs to the initial improper sheet α < 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α > argγ. 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  + C+ + Γ̃+ as shown in Fig. 3. Applying Cauchy's residue theorem on this contour yields

ψ(x,z,ω)= trapped+ revealedleaky+Γ̃+PMLmodes.
(6)

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 α < 0 (i.e., the shaded region in Fig. 3) the sign of α 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, αα > 0 (hatched region) and αα < 0 (white region)—see Fig. 2. The former and the latter are hence such that α > 0 and α < 0 in the proper Riemann sheet (i.e., α > 0), and conversely in the improper one (α < 0). Consequently, the revealed part of the sheet α < 0 is such that α > 0. As expected leaky waves are outgoing in the transverse direction.

In the elastodynamic case, two transverse wavenumbers must be considered: one for the longitudinal transverse wavenumber, αl, and one for the shear transverse wavenumber, αs. These are defined from αl2 = kl2β2 and αs2 = ks2β2, with kl= ω/cl and ks = ω/cs. 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, αl = 0 and αs = 0. This yields four Riemann sheets, labelled 1 for (αl > 0, αs > 0), 2 for (αl > 0, αs < 0), 3 for (αl < 0, αs < 0), and 4 for (αl < 0, αs > 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 argγ. 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.

Fig. 4.

Representation in the elastodynamic case of the proper sheet with PML (α̃l > 0, α̃s > 0): (a) lossy medium, (b) lossless. Same legend as in Fig. 3, except that sheet 1 (αl > 0, αs > 0) is now decomposed into hatched (αl > 0, αs > 0), white (αl < 0, αs < 0) and pea (αl < 0, αs > 0) regions. Shaded regions: revealed parts of the initial improper sheets 2, 3, and 4.

Fig. 4.

Representation in the elastodynamic case of the proper sheet with PML (α̃l > 0, α̃s > 0): (a) lossy medium, (b) lossless. Same legend as in Fig. 3, except that sheet 1 (αl > 0, αs > 0) is now decomposed into hatched (αl > 0, αs > 0), white (αl < 0, αs < 0) and pea (αl < 0, αs > 0) regions. Shaded regions: revealed parts of the initial improper sheets 2, 3, and 4.

Close modal

Similarly to Sec. 3, the sign of αl and αs can be determined in the accessible regions of sheets 2, 3, and 4. This yields the following results:

Region2:αl<0,αs>0(αl>0,αs<0),Region3:αl>0,αs>0(αl<0,αs<0),Region4:αl>0,αs>0(αl<0,αs>0).
(7)

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

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.

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

1.
P.
Malischewsky
,
Surface Waves and Discontinuities
(
Elsevier
,
Amsterdam
,
1987
).
2.
R. E.
Collin
,
Field Theory of Guided Waves
(
IEEE Press
,
Wiley
, New York,
1991
).
3.
B. N.
Pavlakovic
,
M. J. S.
Lowe
, and
P.
Cawley
, “
High-frequency low-loss ultrasonic modes in imbedded bars
,”
J. Appl. Mech.
68
,
67
75
(
2001
).
4.
S.
Kim
and
J. E.
Pasciak
, “
The computation of resonances in open systems using a perfectly matched layer
,”
Math. Comput.
78
,
1375
1398
(
2009
).
5.
F.
Treyssède
,
K. L.
Nguyen
,
A.-S.
Bonnet-BenDhia
, and
C.
Hazard
, “
Finite element computation of trapped and leaky elastic waves in open stratified waveguides
,”
Wave Motion
51
,
1093
1107
(
2014
).
6.
H.
Derudder
,
F.
Olyslager
,
D.
De Zutter
, and
S.
Van den Berghe
, “
Efficient mode-matching analysis of discontinuities in finite planar substrates using perfectly matched layers
,”
IEEE Trans. Ant. Propag.
49
,
185
195
(
2001
).
7.
D.
Marcuse
,
Theory of Dielectric Optical Waveguides
(
Academic Press
,
New York
,
1991
).
8.
T.
Tamir
and
A. A.
Oliner
, “
Guided complex waves. Part 1: Fields at an interface
,”
Proc. Inst. Electr. Eng.
110
,
310
324
(
1963
).
9.
J.
Hu
and
C. R.
Menyuk
, “
Understanding leaky modes: Slab waveguide revisited
,”
Adv. Opt. Photon.
1
,
58
106
(
2009
).
10.
W. C.
Chew
and
W. H.
Weedon
, “
A 3D perfectly matched medium from modified Maxwell's equations with stretched coordinates
,”
Microwave Opt. Technol. Lett.
7
,
599
604
(
1994
).
11.
V.
Maupin
, “
The radiation modes of a vertically varying half-space: A new representation of the complete Green's function in terms of modes
,”
Geophys. J. Int.
126
,
762
870
(
1996
).
12.
M.
Mazzotti
,
I.
Bartoli
,
A.
Marzani
, and
E.
Viola
, “
A coupled SAFE-2.5D BEM approach for the dispersion analysis of damped leaky guided waves in embedded waveguides of arbitrary cross-section
,”
Ultrasonics
53
,
1227
1241
(
2013
).
13.
A. L.
Kurkjian
, “
Numerical computation of individual far-field arrivals excited by an acoustic source in a borehole
,”
Geophysics
50
,
852
866
(
1985
).
14.
F.
Olyslager
, “
Discretization of continuous spectra based on perfectly matched layers
,”
SIAM J. Appl. Math.
64
,
1408
1433
(
2004
).
15.
E.
Becache
,
S.
Fauqueux
, and
P.
Joly
, “
Stability of perfectly matched layers, group velocities and anisotropic waves
,”
J. Comput. Phys.
188
,
399
433
(
2003
).