# A general spectral collocation method for computing the dispersion relations of guided acoustic waves in multilayer dissipative structures

Mathieu Maréchal, Alan Geslain, Jean-Philippe Groby, Vicente

Romero-García, Olivier Dazel

► To cite this version:

Mathieu Maréchal, Alan Geslain, Jean-Philippe Groby, Vicente Romero-García, Olivier Dazel. A general spectral collocation method for computing the dispersion relations of guided acoustic waves in multilayer dissipative structures. Journal of Applied Physics, 2025, 137 (10), pp.104902. 10.1063/5.0242405 . hal-05406069

HAL Id: hal-05406069

https://hal.science/hal-05406069v1

Submitted on 9 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A general spectral collocation method for computing the dispersion relations of guided acoustic waves in multilayer dissipative structures

Mathieu Maréchal*,<fnref n="1" /> Alan Geslain,<fnref n="1" /> Jean-Philippe Groby,<fnref n="1" /> Vicente Romero-García,<fnref n="2" /> and Olivier Dazel<fnref n="1" />

A spectral collocation method is proposed to compute the complex wavenumber-real frequency dispersion relations of guided acoustic waves in multilayer structures involving dissipative materials. The nature of these dissipative materials is initially considered to be arbitrary, i.e., poroelastic, viscoelastic, or viscoacoustic. For a given frequency, the complex wavenumbers as well as the physical fields, which are further used to evaluate the Poynting vectors and analyze the energy flux, are obtained by solving a generalized eigenvalue problem. The latter arises from a set of discretized equations of motion and appropriate boundary (coupling) conditions. These equations of motion and boundary (coupling) conditions are imposed by the nature of the material composing each layer of the structure. A focus is made on poroelastic layers. The dispersion relation of a two-layer elastic-poroelastic structure is analyzed, as well as the energy flows in the structure. The results as calculated with the present spectral collocation method are validated against those obtained with a classical complex root-finding (Müller) method and experiments.

Keywords: Spectral collocation method, dissipative multilayer systems, complex wavenumber-real frequency dispersion relation, guided waves.

## I. INTRODUCTION

Dealing with guided acoustic waves in multilayer structures of dissipative (viscoelastic, viscothermal, or poroelastic) material requires the evaluation of complex wavenumber-real frequency dispersion relations. Solving for these dispersion relations consists in finding complex roots of complex characteristic equations, which often arise from complex matrix determinants and can rarely be solved analytically. Various methods<fnref n="1" />,<fnref n="2" />, among which the Müller algorithm³, are commonly used to achieve this challenging task. These methods usually rely on preliminary semi-analytical derivations, require multiple initial guesses, and are iterative, so they can get stuck in local minima. Computing the complex wavenumber-real frequency dispersion relation in a reliable and robust way is thus of primary importance. Commercial software such as Disperse⁴ or open-source software such as Dispersion Calculator⁵,⁶ are available for calculating dispersion curves. They use mode tracing or root-finding methods and focus mainly on multilayer elastic structures. Alternatively, several numerical methods can be used, such as the Finite Element Method⁷, the Semi-Analytical Finite Element method (SAFE)⁸⁻¹⁰, or spectral methods¹¹. Spectral methods are especially relevant in this framework. They are numerical approximation techniques that seek the solution of a differential equation using a finite series of infinitely differentiable basis functions with spectral accuracy¹².

The propagation and dispersion of waves in multilayer systems have long been investigated¹³,¹⁴. Spectral methods have been developed to cope with wave radiation problems, but have also emerged in the last decades to solve dispersion relations¹⁵. Effectively, the numerical scheme boils down to discretizing the geometry along the layer thickness and solving an eigenvalue problem. Here, a collocation method is used to discretize the geometry. A finite-dimension polynomial basis is chosen and the collocation points are distributed according to the roots of this polynomial basis. Leaky modes appear in practice when multilayer structures are surrounded by one or two half-spaces. Accounting for this phenomenon brings us closer to practical considerations regarding these systems, as they would radiate part of the acoustic energy back to the surrounding medium, although it adds hardships to the modeling of the system; the waves modeled in these half-spaces should satisfy some radiation condition. These can be imposed by implementing perfectly matched layers¹⁶, by using a bounded mapping, and applying a Spectral Collocation Method (SCM) to all layers of the extended system¹⁷, or by using a suitable change of variable and solving a non-linear eigenvalue problem¹⁸. This last method inspires the current work. We propose a SCM suitable for solving complex wavenumber-real frequency dispersion relations in multilayer systems involving dissipative media. As a way of example, we will focus on multilayer poroelastic systems. Mechanical wave propagation in bulk poroelastic materials are well understood and modeled thanks to the seminal contributions of M.A. Biot¹⁹⁻²² and subsequent contributions, notably by D.L. Johnson et al. ²³ and J.-F. Allard et al. ²⁴,²⁵ concerning the modeling of viscothermal losses. Although poroelastic materials are bounded in practice and are usually encapsulated in

<footnote n="1">¹) Laboratoire d'Acoustique de l'Université du Mans (LAUM), UMR 6613, Institut d'Acoustique - Graduate School (IA-GS), CNRS, Le Mans Université, France</footnote>

<footnote n="2">²) Instituto Universitario de Matemática Pura y Aplicada (IUMPA), Universitat Politècnica de Valencia, Camino de Vera s/n, 46022 Valencia, Spain</footnote>

<footnote>(*Corresponding author address: mathieu.marechal@univ-lemans.fr)</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: a) Sketch of the geometry of the multilayer system and b) zoom on the n-th layer.</figcaption></figure>

more complex multilayer structures combining layers of materials of different nature in geophysics and acoustics, there are only a few studies on the dispersion relation of guided acoustic waves in poroelastic$^{26–28}$, porous$^{29,30}$ or anisotropic poroelastic layers$^{31}$. The SCM has already been used to solve the dispersion relation of guided acoustic waves in elastic layered rings$^{32,33}$ or anisotropic elastic multilayer systems$^{34}$, but has never been used to cope with multilayer systems involving dissipative materials of arbitrary nature.

This work aims at filling this gap and is divided into two sections. First, the SCM is derived for generic multilayer configurations possibly coupled to two surrounding identical fluid half-spaces. A two-layer elastic-poroelastic structure is then considered by way of example. The dispersion relation computed with the SCM method is compared to that obtained with a classical root-finding (Müller) method and that measured. The energy flow is further analyzed for this two-layer configuration, with very little additional computational cost. In the last section, we experimentally analyze the two-layer structure by retrieving its complex dispersion relation from displacement measurements using the SLaTCoW (Spatial LAplace Transform for COmplex Wavenumber recovery). The previously obtained numerical results can thus be correlated with the data of the dispersion relation.

## II. SOLVING THE DISPERSION RELATIONS FOR MULTILAYER STRUCTURES

### A. Spectral Collocation Method

Figure 1a depicts the general geometry of the multilayer structure considered in this article. It is composed of N homogeneous layers of arbitrary nature. The n-th layer has a thickness of $h^{(n)}$ and occupies the domain $\Omega^{(n)}$. The lower and upper interfaces are denoted $\Gamma^{(n-1)}$ and $\Gamma^{(n)}$. The multilayer structure is assumed to be two-dimensional, of infinite extent along $x_1$, and bounded along the $x_2$-axis between 0 and the total structure thickness $h = \sum_{n=1}^{N} h^{(n)}$. Outgoing waves radiate above and below the structure in two fluid half-spaces $\Omega^{(0^\pm)}$, occupied by the same fluid.

We assume an implicit time dependence $e^{-i\omega t}$. The physical fields associated to the guided waves in the n-th layer propagating along the positive $x_1$-axis take the form

$$ \Theta^{(n)}(\mathbf{x}) = s^{(n)}(k_1, x_2)e^{ik_1 x_1}, \quad (1) $$

where $s^{(n)}(k_1, x_2)$ is the spatial Laplace transform of $\Theta^{(n)}(\mathbf{x})$ and $k_1$ is the complex valued $x_1$-component of the wavenumber. Note that the spatial Laplace transform is used instead of the spatial Fourier transform because $k_1$ is complex-valued. because $k_1$, the $x_1$ component of the wavenumber is complex-valued. The physical fields that compose $\Theta^{(n)}(\mathbf{x})$ depend on the nature of the material occupying the n-th layer. Introducing $\Theta^{(n)}(\mathbf{x})$ in the corresponding wave equations (of second order in space in this work) reduces to a second-order differential operator $\mathcal{L}^{(n)}$ in $x_2$, that depends on the material properties and can further be developed as a second-order polynomial in
$k_1$

$$ \mathcal{L}^{(n)}(x_2)\mathbf{s}^{(n)}(k_1, x_2) = \left(k_1^2 \mathcal{L}_2^{(n)}(x_2) + k_1 \mathcal{L}_1^{(n)}(x_2) + \mathcal{L}_0^{(n)}(x_2)\right) \mathbf{s}^{(n)}(k_1, x_2) = \mathbf{0}, \qquad (2) $$

where $\mathcal{L}_j^{(n)}$ are the coefficients of the polynomial expansion $j = 0, 1, 2$. Please note that this operator is nothing but the spatial Laplace transform of the equations of motion.

The interface (boundary) conditions at $\Gamma^{(n)}$ between the $n+1$-th and the $n$-th layers again depend on the nature of these two layers, but can be formally written in the form

$$ C^{(n)-} \mathbf{s}^{(n)}(k_1)|_{\Gamma^{(n)}} - C^{(n)+} \mathbf{s}^{(n+1)}(k_1)|_{\Gamma^{(n)}} = \mathbf{0}, \qquad (3) $$

where $C^{(n)\pm}$ are interfaces operators, which involve differential operators of maximum second-order in $x_2$. Following this procedure and introducing $\mathbf{S} = (\mathbf{s}^{(0^-)} \dots \mathbf{s}^{(n)} \dots \mathbf{s}^{(0^+)})^\mathrm{T}$, with $\mathrm{T}$ denoting transposition, the problem can be cast in matrix form, whose coefficients are differential operators in $x_2$,

$$ \mathbf{K} \mathbf{S}(k_1, x_2) = \begin{pmatrix} C^{(0^-)} & -C^{(0^+)} & 0 & & & \\ 0 & \mathcal{L}^{(1)}(x_2) & 0 & & & \\ 0 & C^{(1^-)} & -C^{(1^+)} & 0 & & \\ & 0 & \mathcal{L}^{(2)}(x_2) & 0 & & \\ & & \ddots & \ddots & \ddots & \\ & & 0 & \mathcal{L}^{(N)}(x_2) & 0 & \\ & & 0 & C^{(N^-)} & -C^{(N^+)} & \end{pmatrix} \mathbf{S}(k_1, x_2) = \mathbf{0}. \qquad (4) $$

Note that the equations of motion in both upper and lower half-spaces are not required in Eq.(4). This system only depends on the values of the two spatial Laplace transforms $s^{(0^\pm)}$ at the upper and lower interfaces. Instead, we specify the nature of the two half-spaces, i.e. a fluid medium, and make use of the Sommerfeld condition, to explicitly give the form of the spatial Laplace transforms $\tilde{p}^{(0^\pm)}(k_1, x_2)$ of the pressure fields $p^{(0^\pm)}(\mathbf{x})$

$$ \begin{aligned} \tilde{p}^{(0^-)}(k_1, x_2) &= A^{(0^-)} e^{-ik_2^{(0)} x_2}, \quad \in \Omega^{(0^-)}, \\ \tilde{p}^{(0^+)}(k_1, x_2) &= A^{(0^+)} e^{ik_2^{(0)} (x_2 - h)}, \quad \in \Omega^{(0^+)}, \end{aligned} \qquad (5) $$

with $k_2^{(0)} = \sqrt{(k^{(0)})^2 - k_1^2}$, such that $\mathrm{Re}(k_2^{(0)}) \ge 0$, where $k^{(0)} = \omega/c_f$ is the wavenumber in the fluid and $c_f$ the wave velocity. Thus, $s^{(0^\pm)}$ appearing in Eq.(4) reduces to $s^{(0^\pm)} = A^{(0^\pm)}$. $C^{(0^-)}$ and $C^{(N^+)}$ are therefore no longer operators and reduce to $C^{(0^-)} = C^{(0^-)}$ and $C^{(N^+)} = C^{(N^+)}$. In our fluid case, the interface conditions at $\Gamma^{(0)}$ and $\Gamma^{(N)}$ involve the spatial Laplace transform of the pressure and its first spatial derivative with respect to $x_2$ whatever the nature of the 1-st and N-th layers. Both matrices $C^{(j)}$, $j = (0^-, (N^+)$, can thus be formally written as

$$ C^{(j)} = k_1 C_1^{(j)} + C_0^{(j)} + k_2^{(0)} C_{1'}^{(j)}, \qquad (6) $$

with the subscript denoting the submatrices rearranged according to the power order $k_1$. Index $1'$ corresponds to the term in $k_2^{(0)}$. Introducing Eqs.(2) and (6) in Eq. (4) leads to

$$ (k_1^2 \mathbf{K}_2 + k_1 \mathbf{K}_1 + \mathbf{K}_0 + i k_2^{(0)} \mathbf{K}_{1'}) \mathbf{S} = \mathbf{0}. \qquad (7) $$

This eigenvalue problem is non-linear because of the presence of the term $k_2^{(0)}$. Making use of the following changes of variable $k_1 = k^{(0)}(\gamma + \gamma^{-1})/2$ and $k_2^{(0)} = k^{(0)}(\gamma - \gamma^{-1})/2i^{35,36}$ and companion linearization$^{37}$, a generalized eigenvalue problem is formed and solved for $\gamma$, from which the dispersion relation can be calculated. More details on this procedure can be found in$^{18}$. Note that the first or the last column of Eq.(4) vanishes when radiation is absent on either side of the multilayer system. The situation is different when both fluid half-spaces disappear. The eigenvalue problem Eq.(7) is no longer non-linear and the problem can directly be solved after companion linearization. Note also that fluid half-spaces are only considered for the sake of simplicity and that this procedure
can be adapted to any kind of material occupying the half-spaces, as described recently $^{38}$. Note finally that this procedure avoids discretization of these two half-spaces.

To numerically solve the problem, the spatial Laplace transform of the physical fields in each n-th layer are
expanded as

$$
\mathbf{s}^{(n)}(x_2) \approx \underline{s}^{(n)}(x_2) = \sum_{m=0}^{M} \alpha_{m}^{(n)} \underline{\psi}_{m}(\xi^{(n)}), \qquad (8)
$$

where $\alpha_m$ are coefficients of the polynomial expansion, $\underline{\psi}_m$ are the Chebyshev polynomial of order $m$, and $\xi^{(n)} = 2(x_2 - \sum_{j=1}^{n-1} h^{(j)})/h^{(n)} - 1, \in [-1, 1]$. Underlined variables represent discrete vectors and double-underlined variables discrete matrices. The goal of this approximation is to discretize the differential operators in $x_2$ appearing in the problem. To do so, the SCM is employed. $\xi^{(n)}$ is discretized on the roots of the Chebyshev polynomials, $\underline{\xi}_j^{(n)} = \cos(\frac{j\pi}{M^{(n)}})$ with $j = 0...M^{(n)}$. Please note that $\underline{\xi}_j^{(n)}$ discretely runs from 1 to -1 with increasing $j$, while $\xi^{(n)}$ continuously runs over the same interval but in the opposite direction, i.e. from -1 to 1. These collocation points/nodes form a non-uniform grid along the thickness of each layer. This discretization is usually preferred because the collocation points are clustered at the interfaces$^{12}$. The locations of these nodes are presented in Fig. 1b for the n-th layer with an arbitrary value $M^{(n)} = 8$. The differential operators in the n-th layer are thus represented by differentiation matrices (DMs), such that

$$
\begin{align}
\mathbf{s}^{(n)}(x_2) &\rightarrow \underline{\mathcal{I}} \underline{s}^{(n)}(\underline{\xi}_j), \qquad \partial_2 \mathbf{s}^{(n)}(x_2) \rightarrow \left(\frac{2}{h^{(n)}}\right) \underline{\mathcal{D}}_2 \underline{s}^{(n)}(\underline{\xi}_j) = \underline{\mathcal{D}}_2^{(n)} \underline{s}^{(n)}(\underline{\xi}_j), \nonumber \\
\partial_{22} \mathbf{s}^{(n)}(x_2) &\rightarrow \left(\frac{2}{h^{(n)}}\right)^2 \underline{\mathcal{D}}_{22} \underline{s}^{(n)}(\underline{\xi}_j) = \underline{\mathcal{D}}_{22}^{(n)} \underline{s}^{(n)}(\underline{\xi}_j). \tag{9}
\end{align}
$$

where $\underline{\mathcal{I}}$ is the identity matrix and $\underline{\mathcal{D}}_2$ and $\underline{\mathcal{D}}_{22}$ are respectively the first-order and second-order normalized DMs.
Contrary to an interpolation with $\underline{\mathcal{C}}$hebyshev polynomials, DMs interpolate the functions and their spatial derivatives directly and thus do not involve the polynomial expansion coefficients per se.$^{39}$. Note that the normalized DMs
are denoted as $\underline{\mathcal{D}}_2^{(n)}$ and $\underline{\mathcal{D}}_{22}^{(n)}$ in the following. These matrices are square and have a size of $(M \times M)$. This size also
directly depends on the order of the considered Chebyshev polynomials. The formulation prevents the occurrence
of roundoff errors, which is of particular importance when SCM is applied to solve eigenvalue problems$^{40}$.

The linear operator $\mathcal{L}^{(n)}$ is discretized as a matrix $\underline{\mathcal{L}}^{(n)}$ with size $(P^{(n)}(M^{(n)} - 1), P^{(n)}(M^{(n)} + 1))$, where $P^{(n)}$ is the number of physical fields required to describe the wave propagation in the layer and again depends on the nature of the material this layer is composed of. The first and the last rows of the DMs are removed because they are used to apply the interface conditions at $\Gamma^{(n+1)}$ and $\Gamma^{(n)}$ respectively. The interface operators $\mathcal{C}^{(n)\pm}$ are effectively expressed in terms of this first row of the DMs of the n-th layer and this last row of the DMs of the $(n+1)$-th layer. These rows respectively correspond to the first, i.e. $\underline{\xi}_{M^{(n)}} = -1$ and $\xi^{(n)} = 1$, and last collocation nodes of the n-th and $(n+1)$-th layers, i.e. $\underline{\xi}_0 = 1$ and $\xi^{(n+1)} = -1$. These vectors are subsequently denoted $\underline{\mathcal{I}}^{(n+1)-}$, $\underline{\mathcal{I}}^{(n)+}$, $\underline{\mathcal{D}}_2^{(n+1)-}$, $\underline{\mathcal{D}}_2^{(n)+}$, $\underline{\mathcal{D}}_{22}^{(n+1)-}$, and $\underline{\mathcal{D}}_{22}^{(n)+}$. The interface matrix $\underline{\mathcal{C}}^{(n)\mp}$ is of size $(P^{(n)} + P^{(n+1)}, P^{(n)}(M^{(n)} + 1))$ and the interface matrix $\underline{\mathcal{C}}^{(n)\pm}$ is of size $(P^{(n)} + P^{(n+1)}, P^{(n+1)}(M^{(n+1)} + 1))$.

Finally, the full $(\sum_{n=0}^{N} M^{(n)} P^{(n)} + 2, \sum_{n=0}^{N} M^{(n)} P^{(n)} + 2)$-matrix system for the general multilayer system is

$$
\underline{\underline{K}} \underline{\underline{S}} = \begin{pmatrix}
\underline{\underline{C}}^{(0-)} & -\underline{\underline{C}}^{(0+)} & 0 \\
0 & \underline{\underline{L}}^{(1)} & 0 \\
0 & \underline{\underline{C}}^{(1-)} & -\underline{\underline{C}}^{(1+)} & 0 \\
& 0 & \underline{\underline{L}}^{(2)} & 0 \\
& \ddots & \ddots & \ddots \\
& 0 & \underline{\underline{L}}^{(N)} & 0 \\
& 0 & \underline{\underline{C}}^{(N-)} & -\underline{\underline{C}}^{(N+)}
\end{pmatrix} \begin{pmatrix}
\underline{\underline{A}}^{(0-)} \\
\underline{\underline{s}}^{(1)} \\
\underline{\underline{s}}^{(2)} \\
\vdots \\
\underline{\underline{s}}^{(N)} \\
\underline{\underline{A}}^{(0+)}
\end{pmatrix} = \underline{\underline{0}}. \quad (10)
$$

This system is cast in the form of a generalized eigenvalue problem that is solved by traditional eigenvalue solvers
at each frequency. The dispersion relation for the acoustic guided wave in the dissipative multilayer system is
calculated, repeating the procedure for each frequency $\omega$.
## B. Solving the dispersion relation with Müller algorithm

The dispersion relation of multi-layer dissipative structure are more commonly calculated with secant-based algorithm as the Müller algorithm. The latter relies on a semi-analytic description of the fields. These fields are best written in each layer, independently of the nature of the material it is composed of, using the potentials, i.e. decomposing the field $\theta = \nabla\varphi + \nabla \times \psi$ into an irrotational component $\varphi$ and a transverse component $\psi = \psi e_3$. The spatial Laplace transform of the potentials in the n-th layer are entirely described by

$$ \tilde{\theta}_n(k_1, x_2) = A_X^{(n)+} e^{ik_{2,X}^{(n)}(x_2-h^{(n)})} + A_X^{(n)-} e^{-ik_{2,X}^{(n)}(x_2-h^{(n)})}, \quad (11) $$

with $A_X^{(n)\pm}$ the up and down-going wave amplitudes, $k_{2,X}^{(n)} = \sqrt{(k_X^{(n)})^2 - k_1^2}$, such that $\mathrm{Re}(k_{2,X}^{(n)}) \ge 0$, and X refers to the type of wave, i.e. shear, compressional or acoustic. If necessary, these Laplace transform representations are complemented by those of the pressure field in the fluid half-spaces provided in Eq. (5).

The boundary conditions are applied at each interface and the set of equations is cast into the following matrix form $M(k_1)q = 0$, with $q$ a vector containing the amplitudes of each field/potential. The modes of the system correspond to the wavenumbers $k_1$ for which the matrix $M(k_1)$ is singular, i.e.

$$ \det(M(k_1)) = 0. \quad (12) $$

These wavenumbers are the complex roots of an implicit complex-valued equation. They are evaluated by running our own implementation of the Müller algorithm³ from 3 initial guesses. This equation is solved iteratively for each frequency $\omega$, giving the dispersion relation.

In addition, $q$ can be obtained by computing the nullspace of $M^{41}$ (This null space is computed using singular value decomposition as embedded in the `lINALg null_space` function from the Python package SciPy). The amplitudes of the potentials in the layers are thus evaluated, enabling all the physical fields that depend on these potentials to be calculated. However, the $M$ matrix becomes increasingly tedious to write as the number of layers increases, and even more so when it comes to calculating the roots of its determinant and nullspaces. Nevertheless, this approach remains a fairly effective way of providing a reference solution for validating numerical results when studying a two-layer system.

# III. APPLICATION TO A TWO-LAYER POROELASTIC-ELASTIC STRUCTURE

In this section, the calculation of the dispersion relation for guided waves in a two-layer structure, i.e. $N = 2$, is used as an example of the application of the present SCM. This two-layer structure consists of a $h^{(1)}$-thick poroelastic layer, coated on one side by a $h^{(2)}$-thick aluminum plate, that is radiating in a fluid half-space, and rigidly backed on the other side as depicted in Fig. 2a. The SCM results are analyzed and compared to those obtained with the Müller root-finding method.

## A. Numerical scheme

The pressure field $p^{(1)}$ and the two components of the frame displacement, $u_{1,s}^{(1)}$ and $u_{2,s}^{(1)}$, i.e. $P^{(1)} = 3$, as well as the two-components of the elastic displacement, $u_1^{(2)}$ and $u_2^{(2)}$, i.e. $P^{(2)} = 2$, are used to model the propagation in the poroelastic and in the elastic layers respectively. The $\{u_s, p\}$ formulation is used⁴² for the poroelastic medium. The more usual displacement formulations ($\{u_s, u_f\}^{19}$ or $\{u_s, w\}^{21}$) cannot be used in the present case, because the number of physical fields they rely on, i.e. 4, is different from that of the potentials in the layer, i.e. 3.

In total, 5 physical fields inside the two-layer structure and an additional amplitude for the acoustic wave radiated towards $x_2 \to \infty$ are needed to solve for the dispersion relation. Each domain is discretized on $M^{(j)}$, $j = 1, 2$, collocation points.

In more details, the $\{u_s, p\}$ formulation of the Biot theory⁴² is employed to model the propagation in the homogeneous isotropic poroelastic layer. The coupled equations of motion read as

$$ \left\{ \begin{array}{l} \nabla \cdot \hat{\sigma}_s^{(1)} + \hat{\rho} \omega^2 \mathbf{u}_s^{(1)} + \beta \nabla p^{(1)} = 0, \\ \frac{\nabla^2 p^{(1)}}{\rho_{22} \omega^2} - \frac{\beta}{\phi^2} \nabla \cdot \mathbf{u}_s^{(1)} + \frac{1}{R} p^{(1)} = 0, \end{array} \right. \quad (13) $$
where $\hat{\boldsymbol{\sigma}}_s^{(1)} = \hat{A}\nabla \cdot \boldsymbol{u}_s^{(1)}\boldsymbol{I} + 2N_s\boldsymbol{\epsilon}_s$ is the in vacuo solid stress tensor, with $\boldsymbol{\epsilon}_s$, the solid-phase strain tensor. The parameters are generally expressed in terms of the elastic coefficients $P, Q$ and $R$ and the effective densities $\rho_{11}, \rho_{12}$ and $\rho_{22}$ initially introduced by Biot$^{19,43}$ as

$$ \hat{A} = P - 2N_s - \frac{Q^2}{R}, \quad \hat{\rho} = \rho_{11} - \frac{\rho_{12}^2}{\rho_{22}}, \quad \beta = \phi \left( \frac{\rho_{12}}{\rho_{22}} - \frac{Q}{R} \right). \qquad (14) $$

These parameters depend on the frame properties and on the effective density $\rho_{eq}(\omega)$ and bulk modulus $K_{eq}(\omega)$ of the fluid phase as reminded in A. In short, a poroelastic material is entirely described by the porosity $\phi$, the shear modulus of the frame $N_s$, the Poisson ratio $\nu$, the density of the frame $\rho_1$, the viscous and thermal characteristic lengths $\Lambda$ and $\Lambda'$, tortuosity $\alpha_\infty$, and the flow resistivity $R_f$ when saturated by a light fluid$^{19,23,24}$. All these parameters are real valued except for the shear modulus, which incorporates a viscoelastic damping factor.

Introducing $\boldsymbol{s}^{(1)} = (\tilde{u}_1^s, \tilde{u}_2^s, \tilde{p})^\mathrm{T}$, the spatial Laplace transform of Eq. (13) is expanded in the form of Eq. (2) with

$$ \begin{aligned} \mathcal{L}_0^{(1)} &= \begin{pmatrix} N_s \partial_{22} + \hat{\rho} \omega^2 & 0 & 0 \\ 0 & \hat{P} \partial_{22} + \hat{\rho} \omega^2 & \beta \partial_2 \\ 0 & -\beta / \phi^2 \partial_2 & 1/R + \partial_{22} / (\rho_{22} \omega^2) \end{pmatrix}, \\ \mathcal{L}_1^{(1)} &= \begin{pmatrix} 0 & i \hat{P} \partial_2 & i \beta \\ i(\hat{A} + N_s) \partial_2 & 0 & 0 \\ -i \beta / \phi^2 & 0 & 0 \end{pmatrix}, \quad \text{and} \quad \mathcal{L}_2^{(1)} = \begin{pmatrix} \hat{P} & 0 & 0 \\ 0 & N_s & 0 \\ 0 & 0 & 1 / (\rho_{22} \omega^2) \end{pmatrix}. \end{aligned} \qquad (15) $$

Finally, the discretized equations of motion read as

$$ \underline{L}^{(1)} \underline{s}^{(1)} = \begin{pmatrix} (\omega^2 \hat{\rho} - k_1^2 \hat{P}) \underline{I} + N_s \underline{D}_{22} & i k_1 (\hat{P} - N_s) \underline{D}_{22} & -i k_1 \beta \underline{I} \\ i k_1 (\hat{P} - N_s) \underline{D}_2 & \hat{P} \underline{D}_{22} + (\omega^2 \hat{\rho} - k_1^2 N_s) \underline{I} & \beta \underline{D}_2 \\ \frac{-i \beta \omega^2}{\phi^2} \underline{I} & \frac{-\beta \omega^2}{\phi^2} \underline{D}_2 & \frac{k_1^2}{\rho_{22}} \underline{I} + \frac{\omega^2}{R} \underline{I} \end{pmatrix} \underline{s}^{(1)} = 0. \qquad (16) $$

The equation of motion in the elastic medium is derived from the fundamental elasticity equations$^{44}$ and reads as,

$$ \rho\omega^2\boldsymbol{u}^{(2)} + (\lambda + \mu)\nabla(\nabla\cdot\boldsymbol{u}^{(2)}) + \mu\nabla^2\boldsymbol{u}^{(2)} = 0, \qquad (17) $$

where $\rho$ is the density and $\lambda$ and $\mu$ are the Lamé coefficients. The corresponding discretized equation of motion is

$$ \underline{L}^{(2)} \underline{s}^{(2)} = \left( \frac{\rho^{(2)}}{\mu} \omega^2 \underline{I} + \begin{pmatrix} -k_1^2(\lambda/\mu + 2)\underline{I} + \underline{D}_{22} & i k_1(\lambda/\mu + 1)\underline{D}_2 \\ i k_1(\lambda/\mu + 1)\underline{D}_2 & (\lambda/\mu + 2)\underline{D}_{22} - k_1^2\underline{I} \end{pmatrix} \right) \underline{s}^{(2)} = 0. \qquad (18) $$

In addition, the same procedure can be followed for the equation of motion of a wave propagating in a layer of fluid material. We give its expression in order to provide an exhaustive list of the types of layers that can be modeled using the present method. The Helmholtz equation, $(\nabla^2 + k_0^2)p = 0$, with $k_0$ the bulk wavenumber and $p$ the pressure field in the fluid material, is discretized into the form,

$$ (\underline{D}_{22} - k_1^2 + k_0^2) \underline{p} = 0. \qquad (19) $$

The interface and boundary conditions depend on the vectors $\underline{I}^\pm, \underline{D}_2^\pm$ and $\underline{D}_{22}^\pm$ which corresponds to the last (subscript +) and first row (subscript −) of the DMs. From the bottom to the top of the two-layer system, they are

• Rigid boundary condition at $\Gamma^{(0)}$ which leads to

$$ \boldsymbol{u}^{(1)}(x_2 = 0) = 0; \quad \boldsymbol{u}_f^{(1)}(x_2 = 0) \cdot \boldsymbol{n} - \boldsymbol{u}_s^{(1)}(x_2 = 0) \cdot \boldsymbol{n} = 0, \qquad (20) $$

the discretized form of which is,

$$ \underline{C}^{(0)} \underline{s}^{(1)} = \begin{pmatrix} \underline{0} & \underline{I}^+ & \underline{0} \\ \underline{0} & \underline{0} & \underline{I}^+ \\ \frac{\phi}{\rho_{22}} \underline{D}_2^+ & \underline{0} & -\left(1 + \frac{\rho_{12}}{\rho_{22}}\right) \omega^2 \underline{I}^+ \end{pmatrix} \begin{pmatrix} \underline{p}^{(1)} \\ \underline{u}_1^{(1)} \\ \underline{u}_2^{(1)} \end{pmatrix} = \underline{0}. \qquad (21) $$
* Continuity between the poroelastic and the elastic media<sup>45</sup> at $\Gamma^{(1)}$ which leads to

$$ \begin{gathered} \sigma^{t(1)}(x_2 = h^{(1)}) - \sigma^{(2)}(x_2 = h^{(1)}) = 0; \quad \mathbf{u}^{(2)}(x_2 = h^{(1)}) = \mathbf{u}_s^{(1)}(x_2 = h^{(1)}); \\ \mathbf{u}_f^{(1)}(x_2 = h^{(1)}) \cdot \mathbf{n} - \mathbf{u}_s^{(1)}(x_2 = h^{(1)}) \cdot \mathbf{n} = 0, \end{gathered} \qquad (22) $$

where $\sigma^{t(1)} = \hat{\sigma}_s - \phi(1 + Q/R)$ is the total stress tensor in the poroelastic layer. The corresponding discretized interface matrices, arising from

$$ \underline{C}^{(1)+} \underline{s}^{(2)} - \underline{C}^{(1)-} \underline{s}^{(1)} = \underline{C}^{(1)} \left( \underline{p}^{(1)} \quad \underline{u}_1^{(1)} \quad \underline{u}_2^{(1)} \quad \underline{u}_1^{(2)} \quad \underline{u}_2^{(2)} \right)^{\mathrm{T}} = 0, \qquad (23) $$

are cast as

$$ \underline{C}^{(1)} = \begin{pmatrix} -N_s \underline{D}_2^+ & -ik_1 N_s \underline{I} & \underline{0} & \underline{\mu \underline{D}_2^-} & ik_1 \underline{\mu \underline{I}^-} \\ -ik_1 \hat{A} \underline{I}^+ & -\hat{P} \underline{D}_2^+ & \underline{\phi} \left(1 + \frac{Q}{R}\right) \underline{I}^+ & ik_1 \lambda \underline{I}^- & (\lambda + 2\mu) \underline{D}_2^- \\ \underline{0} & -\underline{I}^+ & \underline{0} & \underline{I}^- & \underline{0} \\ \underline{0} & \underline{0} & -\underline{I}^+ & \underline{0} & \underline{I}^- \\ -\frac{\phi}{\rho_{22}} \underline{D}_2^+ \left(1 + \frac{\rho_{12}}{\rho_{22}}\right) \omega^2 \underline{I}^+ & \underline{0} & \underline{0} & \underline{0} & \underline{0} \end{pmatrix}, \qquad (24) $$

with the left-hand side corresponding to $\underline{C}^{(1)-}$ and the right-hand side to $\underline{C}^{(1)+}$.

* Continuity between the elastic and the fluid media at $\Gamma^{(2)}$ which leads to

$$ \begin{gathered} \sigma_{12}^{(2)}(x_2 = h) = 0; \quad \sigma_{22}^{(2)}(x_2 = h) = -p^{(0+)}; \\ u_2^{(2)}(x_2 = h) = u_2^{(0+)}(x_2 = h). \end{gathered} \qquad (25) $$

The corresponding discretized interface matrix $\underline{C}^{(2)} \begin{pmatrix} \underline{u}_1^{(2)} & \underline{u}_2^{(2)} & A^{(0+)} \end{pmatrix}^{\mathrm{T}} = 0$ is

$$ \underline{C}^{(2)} = \begin{pmatrix} \underline{\underline{D}}_2^+ & ik_1 \underline{I}^+ & \underline{0} \\ ik_1 \lambda \underline{I}^+ & (\lambda + 2\mu) \underline{D}_2^+ & \underline{1} \\ \underline{0} & \omega^2 \underline{I}^+ & ik_2^{(0)} \\ \underline{0} & \underline{0} & \underline{\rho}_f \end{pmatrix}. \qquad (26) $$

Combining the discretized interface conditions, Eqs. (21), (24) and, (26), with the discretized equations of motion, Eqs. (16) and (18), leads to the following system

$$ \begin{pmatrix} \underline{C}^{(0)} & 0 & 0 \\ \underline{L}^{(1)} & 0 & 0 \\ \underline{C}^{(1)-} & -\underline{C}^{(1)+} & 0 \\ 0 & \underline{L}^{(2)} & 0 \\ 0 & \underline{C}^{(2)-} & -\underline{C}^{(2)+} \end{pmatrix} \begin{pmatrix} \underline{s}^{(1)} \\ \underline{s}^{(2)} \\ A^{(0+)} \end{pmatrix} = \underline{0}. \qquad (27) $$

This system is further expanded in the form of Eq. (7) and solved for a set of eigenvalues $\gamma$ and associated eigenvectors $\mathbf{S}$ at each frequency $\omega$. Wavenumbers $k_1$ are then obtained from $\gamma$. The dispersion curve is thus calculated over a given frequency range. However, the condition $\mathrm{Re}(k_2^0) \ge 0$ should be imposed to sort out the solution that do not meet it.

Although usual and simple, this two-layer structure is challenging from a numerical point of view because it involves materials with very large impedance contrasts and layers whose thicknesses differ by several orders of magnitude. These problems are partially solved by normalizing each line and rearranging the terms of the matrix Eq. (7), leading to a better numerical conditioning of the system.
<figure><img src="image_3.png" /><figcaption>Figure 2: a) Sketch of the two-layer structure geometry. b) 3D view of the complex wavenumber-real frequency dispersion relation and c) real and imaginary parts of the dispersion relation. Blue dots represents the results as calculated with the SCM and black open markers represents those as calculated with the Müller method. Bulk wavenumbers in the poroelastic and elastic materials are highlighted in yellow and dark red respectively.</figcaption></figure>

## B. Dispersion relation

The material properties of the $h^{(1)} = 52$ mm-thick poroelastic (melanime) layer are listed in Table I, while the properties of the $h^{(2)} = 1$ mm-thick aluminium plate are the density $\rho = 2700$ kg.m$^{-3}$ and the Lamé coefficients $\lambda = 60.75$ GPa and $\mu = 26.03$ GPa. The saturating fluid and that occupying the half-space is air with density $\rho_f = 1.213$ kg.m$^{-3}$, the heat capacity ratio $\kappa = 1.4$, kinematic compressibility $\mu_f = 1.839 \times 10^{-5}$ Pa.s, the Prandtl number $Pr = 0.71$, and adiabatic bulk modulus $\kappa P_0$, where the atmospheric pressure is $P_0 = 1.013 \times 10^5$ Pa. The total two-layer thickness is thus $h = h^{(1)} + h^{(2)} = 53$ mm. $M^{(1)} = 5$ and $M^{(2)} = 11$ collocation points are employed respectively in the poroelastic and elastic layers.

The 3D view of the complex dispersion diagram is depicted in Fig. 2b. To ease readability, real and imaginary parts are depicted separately in Fig. 2c. Please note that this last representation can be misleading because only the modes that lie in the complex wavenumber range depicted in Fig. 2b are represented in Fig. 2c. Cut-on frequencies are thus fictitious. The results calculated using the SCM method, plotted in blue markers, are compared with those obtained using the Müller method, plotted with black open markers. Both methods provide identical results, up to around four digits, therefore validating the present approach. Note that the SCM solutions are used as initial guesses of the Müller algorithm to shorten computational time. In this way, fewer iterations are needed for the
<table><caption>Table I: Parameters of the poroelastic material (melamine foam)</caption><thead><tr><th>Parameter</th><th>φ</th><th>ρ<sub>1</sub></th><th>R<sub>f</sub></th><th>Λ</th><th>Λ'</th><th>α<sub>∞</sub></th><th>ν</th><th>N<sub>s</sub></th></tr><tr><th>Unit</th><th>-</th><th>kg.m<sup>-3</sup></th><th>kPa.s.m<sup>-2</sup></th><th>μm</th><th>μm</th><th>-</th><th>-</th><th>kPa</th></tr></thead><tbody><tr><td>Melamine</td><td>0.98</td><td>6.5</td><td>5.6</td><td>214</td><td>214</td><td>1</td><td>0.24</td><td>11.96 (1 + i0.07)</td></tr></tbody></table>

results of the root-finding method to converge. The main advantages of SCM are that no initial guesses are needed, i.e. as long as the discretization is sufficient to model the full wave behavior in the structure, the dispersion relation is guaranteed to be calculated. In addition, the matrices used for the numerical model are very small compared to those required by other methods such as FEM, and are therefore much less demanding in terms of calculation.

Wavenumbers are perfectly symmetric with both Re($k_1h/2\pi$) = 0, and Im($k_1h/2\pi$) = 0 axis, which is a feature of reciprocal systems. Wavenumbers having a very large slope and a low imaginary part correspond to modes mostly propagating in the aluminium plate. They are very close to the bulk longitudinal $k_l$ and transverse $k_t$ wavenumbers in the aluminium. This is due to the very large contrast between the elastic properties of the aluminium and those of the poroelastic frame. The other branches asymptotically tend toward the bulk acoustic, $k_a$, compression, $k_p$, or shear, $k_s$, wavenumbers of the poroelastic medium at high frequency. At low frequencies, these modes are highly dispersive.

If we compare these results with the dispersion relation depicted in Bc) for the guided waves in a single poroelastic plate twice as thick and in the absence of the elastic plate, we find some essential differences. The branches associated with the elastic plate and the branch clearly associated with the coupling between the elastic and poroelastic layers (referenced as (1) in Fig.Bc)) do not exist. The other branches are similar to those of this single poroelastic layer, although slightly modified by the presence of the elastic plate. This proves a weak coupling between the elastic and poroelastic layers, particularly at high frequencies. At lower frequencies, the coupling between the two layers is stronger. Effectively, the wavelength is much larger than the thickness of at least one of the two layers in these frequency ranges, which favors coupling. The $A_0$ branch observed in the case of a single poroelastic layer disappear in the two-layer system rigidly backed, because this rigid-backing does not allow the existence of such a mode.

### C. Mode shapes and energy fluxes

Beyond the SCM numerical efficiency, a significant advantage of the method is that the eigenvectors associated to each eigenvalue directly provides the associated mode shape. Since these mode shapes, i.e. eigenvectors, are calculated from an eigenvalue problem, their amplitude are defined up to a constant. They are normalized to the value of the pressure field at the bottom of the poroelastic layer, i.e. at the location of the rigid backing, because this condition imposes $\tilde{u}^{(1)}(0) = 0$ (and in particular $\tilde{u}_2^{(1)}(0) = 0$), which is equivalent to a maximum pressure field. From these normalized eigenvectors, any physical fields can be calculated as detailed in C.

Let us consider two solutions $k_1^\pm = \pm|\text{Re}(k_1)|$ symmetric in the dispersion diagram with respect to the axis $\text{Re}(k_1) = 0$ as shown in Figs. 3a-b. The displacement fields and total energy fluxes in each direction are depicted in Fig. 3c-f. The results calculated using the SCM method, plotted in continuous lines, match those obtained using the Müller method, plotted with marker lines. These two modes are identical but propagate in opposite directions, i.e., the mode $k_1^+$ propagates towards positive $x_1$ while the mode $k_1^-$ propagates towards negative $x_1$. Displacement fields and energy fluxes are identical along the $x_2$ axis but are of equal modulus but opposite sign along the $x_1$ axis (see Fig. 3c-d and e-f). Both displacement field and energy flux amplitudes are two to three orders of magnitude smaller in the elastic plate (for $h_1 < x_2 < h$) than in the poroelastic plate. The aluminium coating is almost acting as a rigid boundary condition, as also testified by the fact that the $x_1$ components of the displacement field and energy fluxes are almost symmetric with respect to the axis $x_2 = h_1/2$, while their $x_2$ components are almost antisymmetric along the $x_2$ axis. Note that the fluid phase Poynting vectors depicted in Fig. 3f are of equal modulus but opposite sign, symmetric with respect to the axis $x_2 = h_1/2$, and more than a hundred times smaller than total energy flux in amplitude. The energy of these modes are thus mainly localized in the solid phase of the poroelastic plate.

To go a step further, we evaluated the energy transport velocity $\overline{V}_e = \overline{P}_t/\overline{U}$ as the ratio between the average energy flux in each layer $\overline{P}_t$ over the mean total energy $\overline{U}$. The energy transport velocity of the mode highlighted
<figure><img src="image_4.png" /><figcaption>Figure 3: a) and b) represent the real and imaginary parts of the dispersion relation of the two-layer structure. Two branches are highlighted in yellow and green, as well as two specific solutions, $k_1^-$ (blue) and $k_1^+$ (red), symmetric with respect to the Re($k_1$) = 0 axis. Displacement fields in the $x_1$ c) and $x_2$ direction d) associated with the $k_1^+$ and $k_1^-$ solutions displayed. Total Poynting e) vector and fluid Poynting vector f) with the length of the arrow indicating the amplitude and angle indicating the direction of the vector field. g) Group and energy transport velocities for the yellow branch. h) Energy transport velocity for the green branch. The gray area represents the frequency range over which this velocity is negative.</figcaption></figure>
in orange in Fig. 3a-b is compared to the group velocity $v_g = \partial\omega/\partial\text{Re}(k_1)^{46}$ in Fig. 3g. It is numerically computed by doing finite differences on a branch of the dispersion relation. Both velocities match for frequency ranges where the imaginary part of $k_1$ is low, i.e., for weak attenuation modes$^{47-50}$, but are different for strong attenuation. In particular, the group velocity diverges around the cut-on frequency. This is the reason why the energy transport velocity is usually preferred to study absorbing structures. The energy transport velocity $\bar{V}_e$ of the mode highlighted in green Fig. 3a-b is depicted in Fig. 3h. This energy transport velocity is always positive, while $k_1$ crosses the $\text{Re}(k_1) = 0$ axis. When $\text{Re}(k_1) < 0$ the phase velocity is negative, while the $\bar{V}_e > 0$. This branch (and the symmetric one) has also a negative group velocity as also testified by the fact $\text{Re}(k_1) < 0$ while $\text{Im}(k_1) > 0$ in the shaded area in Fig. 3h. This branch has also an exponentially growing amplitude.

D. Experimental validation

The experimental set-up is depicted in Fig. 4a. The sample is the two-layer structure considered in the previous subsection. The aluminum plate is glued on top of the poroelastic layer. The latter layer is glued on a 3 cm-thick aluminium plate, which mimic the rigid backing. The three elements are 45 cm wide and 85 cm long. The surrounding medium and the saturating fluid is air, at regular pressure and temperature conditions. The sample is positioned perpendicularly to the ground on its longest side as depicted in Fig. 4a to ease the mode excitation in the poroelastic layer. The excitation is provided by a shaker (Brüel and Kjaer type 4810), which is rigidly attached to the sample with a threaded steel rod fixed to the shaker on one side and glued on a 1mm thick aluminum plate of width 10 cm and height 1.5 cm. This plate is cut at the edge opposite to the threaded steel rod and glued to the porous sample, creating a line source at the edge of the sample. The resonance of this part was measured to be 4500 Hz, which is the upper limit of the measured frequency range. The excitation signal is a swept-sine function with 400 points ranging from 1 kHz to 3 kHz. The general layout of the experiments is similar to that used originally in$^{51}$. The key difference in our case is that the excitation is operated directly to the poroelastic layer and not to the top of the aluminum plate, because the modes that corresponds to propagating solutions in the poroelastic layer are encapsulated in between two almost rigid plates (see Section III C). The normal displacement is measured with a laser Doppler vibrometer (Polytec VibroFlex Neo) on a line of length $L = 40$ cm, along the side of the thickness of the porous layer. These measurements are averaged 100 times. The excitation signal and the measured field are interfaced through a computer via a Zürich Instruments acquisition card.

The SLaTCoW method is used to analyze the space-frequency displacement measurements and recover the complex wavenumber-real frequency dispersion relation. This method relies on the Laplace spatial transform of these measurements, which gives access to the real and imaginary components of the complex wavenumber. The difference between this Laplace transform and a correctly chosen ansatz function is finally minimized for each frequency, to recover the dispersion relation. More details on this method are given in Appendix E as well as in the original article$^{51}$. The complex wavenumber-real frequency dispersion relation in Fig. 4b-c depicts the wavenumbers recovered experimentally $k_s$ and simulated $k_n$. Four modes were sought during the minimization with the SLaTCow method: three of them correspond to the three bulk waves due to excitation and the fourth is depicted. The experimentally recovered mode jumps from one branch to the other because of the large modal density. The cut-on frequencies of the modes are clearly visible, especially on the imaginary part of the wavenumber. The lower real wavenumber parts of the numerical dispersion relation are not retrieved with the experimental results. This is because they correspond to a high attenuation of the mode, a part that is inherently difficult to measure in this experimental conditions. Generally, the recovered branches are those around the bulk solid compression wavenumber, corresponding to the wave that we most probably excite in the structure. The difficulties in efficiently exciting all the guided waves supported by this configuration largely explain the missing curves in the dispersion relation. The location and orientation of the excitation line limit possible symmetries and the generation of potential shear waves. This choice was made to avoid the main excitation of guided waves in the elastic plate when excitation and displacement measurement are performed only on this element. In addition, the modes mainly supported by the fluid phase should also be better excited by a loudspeaker and therefore better measured by a microphone. However, we have chosen to focus on the modes mainly supported by the solid phase. Nevertheless, the agreement between experimental results and simulation is quite good, thus validating the method.
<figure><img src="image_5.png" /><figcaption>Figure 4: a) Experimental set-up used for the measurement of the normal displacement field of the two-layer structure, with a close-up of the sample excitation. b) Real and c) imaginary parts of the dispersion relations recovered experimentally from the SLaTCoW method $k_s$ (black crosses) and calculated from the SCM $k_s$ (colored lines and black-to-grey dots) and the bulk wavenumbers in the poroelastic material (magenta lines).</figcaption></figure>

## IV. CONCLUSION

A spectral collocation method (SCM) is proposed to calculate dispersion relations for guided acoustic wave propagation in multilayer structures. The formalism is general enough to be able to consider arbitrary arrangements of (visco-) elastic, poroelastic, and (viscothermal-) fluid layers. Outgoing fluid radiation, leading to Leaky modes, can be accounted for. The fields in each layer is approximated by Chebyshev polynomials and discretized on the corresponding nodes. Interface conditions are implemented to couple layers together. The solution provides the full complex wavenumber-real frequency dispersion relation of the considered structure. In addition, the physical fields in each layer are directly evaluated since they are the eigenvectors associated with the eigenvalue problem giving the dispersion relation. The dispersion diagram calculated with the present SCM is validated against experimental results and that calculated with a root-finding method (Müller method) in a simple two-layer structure consisting of a rigidly backed poroelastic layer covered by a thin aluminum plate radiating in a fluid half-space. The semi-analytical approach underlying the root-finding method is taken a step further to solve for the wave amplitudes as well, providing a reference for the mode shapes calculated with the SCM. The dispersion properties in such structures, as well as the direction of the energy density in each phase of the poroelastic layers using a decoupled expression for the Poynting vectors is analyzed. Results on energy transport velocity are also reported.

This method paves the way of complex wavenumber-real frequency dispersion relation calculation and analysis of more complex structures, like layers with embedded periodic inclusions (elastic or resonant), i.e. metaporolastic surfaces.

## CONFLICT OF INTEREST

The authors have no conflicts to disclose
## AUTHOR CONTRIBUTIONS

**Mathieu Maréchal** Formal analysis (lead); Investigation (equal); Methodology (equal); Software (lead); Visualization (lead); Writing original draft (lead). **Alan Geslain** Writing review & editing (supporting); Investigation (equal). **Jean-Philippe Groby** Formal analysis (equal); Methodology (equal); Supervision (equal); Writing review & editing (lead). **Vicent Romero-Garcia** Formal analysis (supporting); Methodology (supporting); Supervision (equal); Writing review & editing (equal). **Olivier Dazel** Formal analysis (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing review & editing (equal).

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## APPENDICES

### Appendix A: Expression of the Biot elastic coefficients and effective parameters in poroelastic materials

The complex and frequency dependent density and bulk modulus of the fluid phase, which accounts respectively for the viscous and thermal losses are$^{23,24}$ as,

$$ \rho_{\text{eq}}(\omega) = \frac{\rho_f \alpha_{\infty}}{\phi} \left( 1 + \frac{i R_f \phi}{\rho_f \alpha_{\infty} \omega} \sqrt{1 - i \omega \rho_f \mu_f \left( \frac{2 \alpha_{\infty}}{R_f \phi \Lambda} \right)^2} \right), \\ K_{\text{eq}} = \frac{\kappa P_0}{\phi} \left( \frac{\kappa - (\kappa - 1)}{1 + \frac{i R_f \phi}{\rho_f \alpha_{\infty} \text{Pr}\omega} \sqrt{1 - i \text{Pr}\omega \rho_f \mu_f \left( \frac{2 \alpha_{\infty}}{R_f \phi \Lambda'} \right)^2}} \right)^{-1}, \quad (\text{A1}) $$

where $\phi$ is the porosity, $\alpha_{\infty}$ is the tortuosity, $\Lambda$ and $\Lambda'$ are respectively the viscous and thermal characteristic lengths, $R_f$ is the flow resistivity, $\mu_f$ is the dynamic viscosity, $\kappa$ is the heat capacity ratio, $\rho_f$ the density of the saturating fluid, and $P_0$ is the ambient pressure.

The Biot elastic coefficients, $P$, $Q$, and $R$, are more commonly expressed in terms of the frame, $K_b$, effective fluid, $K_{eq}$, and shear, $N_s$, moduli when saturated by a light fluid as

$$ P = K_b + \frac{4}{3} N_s + (1 - \phi)^2 \frac{K_{\text{eq}}}{\phi}; \quad Q = K_f(1 - \phi), \quad R = \phi K_{\text{eq}}, \quad (\text{A2}) $$

while the apparent densities are

$$ \rho_{22} = \phi^2 \rho_{eq}, \quad \rho_{12} = \phi \rho_f - \rho_{22}, \quad \rho_{11} = (1 - \phi) \rho_s - \rho_{12}, \quad (\text{A3}) $$

with $\rho_s$ the density of the solid phase. The coefficients in Eq. (13) are thus

$$ \hat{A} = K_b - \frac{2}{3} N_s + \frac{K_{\text{eq}}}{\phi}((1 + (1 - \phi))), \quad \hat{\rho} = \rho_1 + \phi \rho_f - 2 \phi^2 \rho_{eq} - \frac{\rho_f^2}{\rho_{eq}}, \quad \beta = \frac{\rho_f}{\rho_{eq}} - 2\phi - 1, \quad (\text{A4}) $$

where $\rho_1 = \phi \rho_f + (1 - \phi) \rho_s$ is the apparent density of the frame. Finally, the frame bulk modulus and shear modulus are linked by the Poisson ratio $\nu$ with $K_b = 2N_s(1 + \nu)/3(1 - 2\nu)$.

### Appendix B: Application to a single poroelastic layer surrounded by two fluid half-spaces

The dispersion relation of guided waves in a single layer of poroelastic material surrounded by two identical fluid half-spaces is considered (See Fig. 5a). The poroelastic material properties are those considered in Section III B and
<figure><img src="image_6.png" /><figcaption>Figure 5: a) Sketch of the single layer poroelastic structure geometry. b) 3D view of the complex wavenumber-real frequency dispersion relation and c) real and imaginary parts of the dispersion relation. The branches of the symmetric and antisymmetric modes are represented by blue and red curves respectively. The dispersion relation of the two-layer elatsic-poroelastic system studied in Section III are reminded with black lines.</figcaption></figure>

the layer thickness is $2h^{(1)} = 104$ mm, i.e. twice the thickness of the poroelastic layer considered in the two-layer system.

The discretized equations of motion is that provided in Eq. (16). $M^{(1)} = 11$ collocation points are considered. Only the boundary conditions at the lower, i.e. $\Gamma_0$ at $x_2 = 0$, and upper, i.e. $\Gamma_1$ at $x_2 = h$, interfaces are modified. They read as$^{45}$

$$ \left(1 + \phi + \phi \frac{\rho_{12}}{\rho_{22}}\right) u_2^{(1)} + \frac{\phi^2}{\rho_{22}\omega^2} \partial_2 p^{(1)} = u_2^{(0\pm)}, \quad p^{(1)} = p^{(0\pm)}, \tag{B1} \\ \hat{\sigma}_{22}^s = \left(1 - \phi \left(1 + \frac{Q}{R}\right)\right) p^{(0\pm)}, \quad \hat{\sigma}_{12}^s = 0. $$
The discretized form of these boundary conditions is

$$ \underline{C}_1 = \begin{pmatrix} \frac{\phi^2}{\rho_{22}} \underline{D}_2^+ & \underline{0} & \omega^2 \left(1 + \phi \left(1 + \frac{\rho_{12}}{\rho_{22}}\right)\right) \underline{I}^+ & \frac{ik_2^{(f)}}{\rho_f} \\ \underline{I}^+ & \underline{0} & \underline{0} & -1 \\ \underline{0} & ik_1(\hat{P} - 2N_s)\underline{I}^+ & \hat{P}\underline{D}_2^+ & 1 - \phi \left(1 + \frac{Q}{R}\right) \\ \underline{0} & N_s\underline{D}_2^+ & ik_1 N_s \underline{I}^+ & 0 \end{pmatrix}. \quad (B2) $$

These matrices are then arranged in the same manner as Eq. (10) and the wavenumbers are calculated. Fig. 5b-c depicts the results.

In a similar way as Lamb modes, the different modes can be sorted in symmetric and antisymmetric modes. This is done by comparing the sign of the radiated amplitudes at each edge of the layer. Results are shown in Fig. 5b. All modes necessarily have an imaginary part, both because of the radiation condition, i.e. leakage, and because of viscothermal dissipation as can be noticed in Fig. 5c.

## Appendix C: Evaluation of the physical fields from the SCM solutions

The discrete sets of eigenvectors resulting from the SCM are the spatial Laplace transforms of some physical fields at collocation points, i.e. the locations where the discrete equations of motion (residue) are exactly satisfied. These discrete fields are interpolated into a continuous form using

$$ s^{(n)}(x_2) = \left[ \psi_m^{-1}(\xi) \underline{s}_m^{(n)}(\xi) \right] \psi(\xi^{(n)}) \quad (C1) $$

where the terms in brackets correspond to the coefficients $\alpha_m^{(n)}$ given in Eq. (8). The Chebyshev polynomials $\psi$ is the basis that maps from the discrete grid $\xi$ back to the physical space.

The stress tensor components depend on the direct expressions of the fields and their first order derivative in $x_2$, denoted with a prime symbol in the following. The various fields, differentiated up to any order, can be calculated by replacing $\psi$ by its derivative in Eq. (C1). Explicit expressions are given in the following,

* in the elastic layer, the spatial Laplace transforms of stress tensor components are computed from $\tilde{u}$, with

$$ \begin{aligned} \tilde{\sigma}_{11} &= (\lambda + 2\mu)ik_1\tilde{u}_1 + \lambda\tilde{u}_2', & \tilde{\sigma}_{12} &= \mu(\tilde{u}_1' + ik_1\tilde{u}_2), \\ \tilde{\sigma}_{22} &= \lambda ik_1\tilde{u}_1 + (\lambda + 2\mu)\tilde{u}_2'. \end{aligned} \quad (C2) $$

* in the poroelastic layer, the full solid stress tensor is expressed as $\tilde{\sigma}_s = \tilde{\sigma}_s - \phi(Q/R)\tilde{p}$, with $\tilde{\sigma}_s$ the spatial Laplace transform of the in vacuo solid stress tensor introduced in Section III. The components of this full solid stress tensor are,

$$ \begin{aligned} \tilde{\sigma}_{s,11} &= P k_1 \tilde{u}_{s,1} + \hat{A} \tilde{u}_{s,2}' - \phi \frac{Q}{R} \tilde{p}, & \tilde{\sigma}_{s,12} &= N_s (\tilde{u}_{s,1}' + ik_1 \tilde{u}_{s,2}), \\ \tilde{\sigma}_{s,22} &= \hat{A} k_1 \tilde{u}_{s,2} + P \tilde{u}_{s,2}' - \phi \frac{Q}{R} \tilde{p}. \end{aligned} \quad (C3) $$

The Laplace transform of the fluid stress tensor is calculated as

$$ \tilde{\sigma}_f = -\phi \tilde{p} \delta_{ij}, \quad (C4) $$

and the fluid displacement components are

$$ \tilde{u}_{f,1} = \frac{\phi}{\rho_{22}\omega^2} ik_1 \tilde{p} - \frac{\rho_{12}}{\rho_{22}} \tilde{u}_{s,1}, \quad \tilde{u}_{f,2} = \frac{\phi}{\rho_{22}\omega^2} \tilde{p}' - \frac{\rho_{12}}{\rho_{22}} \tilde{u}_{s,2}. \quad (C5) $$

All physical fields $\theta(\mathbf{x})$ can finally be calculated by the inverse spatial Laplace transform with $\theta(\mathbf{x}) = \tilde{\theta}(k_1, x_2)e^{ik_1 x_1}$.
## Appendix D: Energy flux in the multilayer structure - Poynting theorem in a poroelastic medium

Energy conservation is commonly conducted in the time domain. Let us consider a volume $\Omega$, bounded by $\partial\Omega$, with $n$ the outgoing normal vector. The energy balance or Poynting theorem in a poroelastic materials reads as¹⁹,⁵²

$$ -\int_{\partial\Omega} (-\boldsymbol{\sigma}_s \cdot \boldsymbol{v}_s^*) \cdot \boldsymbol{n} + (-\boldsymbol{\sigma}_f \cdot \boldsymbol{v}_f^*) \cdot \boldsymbol{n} \, dS = \frac{1}{2} \int_{\Omega} \rho_1 \partial_t |\boldsymbol{v}_s|^2 - \rho_{12} \partial_t (\boldsymbol{v}_f - \boldsymbol{v}_s)^2 + \rho_2 \partial_t |\boldsymbol{v}_f|^2 \, d\Omega \\ + \int_{\Omega} \boldsymbol{\sigma}_s : \partial_t \boldsymbol{\varepsilon}_s^* + \boldsymbol{\sigma}_f : \partial_t \boldsymbol{\varepsilon}_f^* - b (\boldsymbol{v}_s - \boldsymbol{v}_f)^2 \, d\Omega, \quad (D1) $$

where $\boldsymbol{v}_s = \partial_t \boldsymbol{u}_s$ and $\boldsymbol{v}_f = \partial_t \boldsymbol{u}_f$ are the velocity fields, $\boldsymbol{\sigma}_s$ and $\boldsymbol{\sigma}_f$ are the stress tensors, $\boldsymbol{\varepsilon}_s = 1/2 (\boldsymbol{\nabla} \boldsymbol{u}_s + \boldsymbol{\nabla}^\mathrm{T} \boldsymbol{u}_s)$ and $\boldsymbol{\varepsilon}_f = 1/2 (\boldsymbol{\nabla} \boldsymbol{u}_f + \boldsymbol{\nabla}^\mathrm{T} \boldsymbol{u}_f)$ are the strain tensors in the solid frame and in the fluid phase respectively, and ⋆ denotes the complex conjugate. Note that a similar expression is provided in⁵⁰, where the alternative 1962 Biot formulation is used.

The left-hand side contains the energy density flux $\boldsymbol{P}_s = -\boldsymbol{\sigma}_s \cdot \boldsymbol{v}_s^*$ and $\boldsymbol{P}_f = -\boldsymbol{\sigma}_f \cdot \boldsymbol{v}_f^*$ in the elastic frame and fluid phase. The right-hand side contains the kinetic energy $K$⁴³ and the internal forces $P_{\mathrm{int}} = \partial_t W + P_{\mathrm{dis}}$. The latter is the sum of the time derivative of the strain energy $W$ and the dissipated power $P_{\mathrm{dis}}$ by the elastic damping and the viscous and thermal losses⁵². In the frequency domain, the kinetic and strain energies read as

$$ K = \frac{1}{2} \rho_1 \partial_t |\boldsymbol{v}_s|^2 - \frac{1}{2} \rho_{12} \partial_t (\boldsymbol{v}_f - \boldsymbol{v}_s)^2 + \frac{1}{2} \rho_2 \partial_t |\boldsymbol{v}_f|^2, \quad (D2) \\ W = \frac{1}{2} \hat{\boldsymbol{\sigma}}_s : \boldsymbol{\varepsilon}_s^* + \frac{1}{2} R \nabla \cdot (Q \boldsymbol{u}_s + R \boldsymbol{u}_f), $$

Finally, the time average total Poynting vector is $\langle \boldsymbol{P}_t \rangle = \langle \boldsymbol{P}_s \rangle + \langle \boldsymbol{P}_f \rangle$, with $\langle \boldsymbol{P}_s \rangle = -\mathrm{Re}(\boldsymbol{\sigma}_s \cdot \boldsymbol{v}_s^*)/2$ and $\langle \boldsymbol{P}_f \rangle = -\mathrm{Re}(\boldsymbol{\sigma}_f \cdot \boldsymbol{v}_f^*)/2$ and the total energy density is $\langle U \rangle = \langle K \rangle + \langle W \rangle = \mathrm{Re}(K)/2 + \mathrm{Re}(W)/2$. Note that a similar procedure can be followed for isotropic elastic or fluid medium⁵³. In particular, the elastic Poynting vector is $\langle \boldsymbol{P} \rangle = \mathrm{Re}(\boldsymbol{\sigma} \cdot \boldsymbol{v}^*)/2$, where $\boldsymbol{\sigma}$ is the stress tensor and $\boldsymbol{v}$ is the velocity in the elastic medium.

The energy transport velocity $V_e$ is evaluated from the latter expressions

$$ V_e = \frac{\langle \boldsymbol{P}_t \rangle \cdot \boldsymbol{n}}{\langle U \rangle}. \quad (D3) $$

When a multi-layer configuration is considered, the latter velocity is averaged over the structure thickness $\bar{V}_e = \bar{P}_t / \bar{U}$, with the average Poynting vector and the average total energy defined as

$$ \bar{P}_t = \sum_n \frac{1}{h_n} \int_{h_{n-1}}^{h_n} \langle \boldsymbol{P}_t^{(n)} \rangle \cdot e_{x_2} \, dx_2, \quad \bar{U} = \sum_n \frac{1}{h_n} \int_{h_{n-1}}^{h_n} \langle U^{(n)} \rangle \, dx_2. \quad (D4) $$

## Appendix E: Reminder of the SLaTCoW method for obtaining complex-valued wavenumbers.

The key point of the SLaTCoW method is the Laplace transform, $\mathcal{L}_{\mathrm{mes}}(z)$ with $z = z_r + iz_t$ the Laplace space coordinate, of the data measured in space at a fixed frequency $\omega$. In doing so, we move from the spatial domain to the complex wavenumber domain. The difference between $\mathcal{L}_{\mathrm{mes}}(z)$ and the spatial Laplace transform of the displacement ansatz along the same line, $\mathcal{L}_{\mathrm{rec}}(z)$, is then minimized in the least square sense to recover both the complex amplitude (magnitude and phase) and the complex wavenumber of each mode assumed to be part of the displacement ansatz.

The experimental Laplace spatial transform $\mathcal{L}_{\mathrm{mes}}(z) = \int_0^L r(x)e^{-izx} dx$ is computed at each frequency $\omega$ from the transfer function $r(x)$ between the excitation and the normal displacement measured along a line $x \in [0, L]$, with $L = 0.3$ m and spacing $d_x = 0.001$ m. The displacement ansatz is assumed to be of the form $w(x) = \sum_{\beta} A^{\beta} e^{ik^{\beta}x}$, with $\beta$ the number of guided-waves, i.e., the number of modes, and where $A^{\beta} = |A^{\beta}|e^{i\phi^{\beta}}$ is the complex amplitude
<figure><img src="image_7.png" /><figcaption>Figure 6: Description of the procedure of the SLaTCoW method: a) Spatial Laplace transform from the experimental data, b) Spatial Laplace transform of the reconstructed wave fields. Real c) and imaginary d) parts of the experimentally retrieved wavenumbers (crosses) compared to the SCM (black-to-grey dots). The frequency at which the Laplace transforms are represented corresponds to the magenta dashed line in the graph.</figcaption></figure>

and $k^β$ is the complex wavenumber of the $β$-th mode.. In practice, $β$ is the smallest integer that gives a good fit between the two Laplace transforms and therefore depends on the frequency under consideration. Its spatial Laplace transform reads as $\mathcal{L}_{\text{rec}}(z) = \sum_{\beta} A^{\beta} Le^{i(k^{\beta}z - \frac{L}{2})} \text{sinc}((k^{\beta} - z)\frac{L}{2})$.

The cost function $\mathcal{C} = (\|\mathcal{L}_{\text{rec}}(z) - \mathcal{L}_{\text{mes}}(z)\|) / \|\mathcal{L}_{\text{mes}}(z)\|$ is finally minimized using the Matlab function `fminsearchbnd` for `{|A^β|, φ^β, Re(k^β), Im(k^β)}`, $\forall \beta$. The initial guess are the recovered values of the minimization parameters at the previous frequency, with the exception of the first frequency. The method is graphically exposed in Figure 6. Figures 6 c) and d) depict the complex dispersion relation as already depicted in Fig. 4. Figure 6a) depicts the spatial Laplace transform of the experimental data and Figure 6b) depicts $\mathcal{L}_{\text{rec}}(z)$ with the parameters derived from minimization of the cost function at the frequency $f = 3.56$ kHz. This frequency is highlighted by the purple dashed line in Figures 6 c) and d). The two Laplace transforms match almost perfectly.
REFERENCES

<footnote n="1">¹P. Kowalczyk, W. Marynowski, Efficient Complex Root Tracing Algorithm for Propagation and Radiation Problems, IEEE Trans. Antennas Propagat. 65 (5) (2017) 2540-2546. doi:10.1109/TAP.2017.2684198.</footnote>

<footnote n="2">²P. S. Dubbelday, Application of a new complex root-finding technique to the dispersion relations for elastic waves in a fluid-loaded plate, SIAM Journal on Applied Mathematics 43 (5) (1983) 1127-1139.</footnote>

<footnote n="3">³D. E. Muller, A Method for Solving Algebraic Equations Using an Automatic Computer, Mathematical Tables and Other Aids to Computation 10 (56) (1956) 208. doi:10.2307/2001916.</footnote>

<footnote n="4">⁴B. Pavlakovic, M. Lowe, D. Alleyne, P. Cawley, Disperse: A General Purpose Program for Creating Dispersion Curves, Springer US, Boston, MA, 1997, pp. 185-192. doi:10.1007/978-1-4615-5947-4_24.
URL https://doi.org/10.1007/978-1-4615-5947-4_24</footnote>

<footnote n="5">⁵A. M. A. Huber, M. G. R. Sause, Classification of solutions for guided waves in anisotropic composites with large numbers of layers, The Journal of the Acoustical Society of America 144 (6) (2018) 3236-3251. arXiv:https://pubs.aip.org/asa/jasa/article-pdf/144/6/3236/14738983/3236\_1\_online.pdf, doi:10.1121/1.5082299.
URL https://doi.org/10.1121/1.5082299</footnote>

<footnote n="6">⁶A. M. A. Huber, Classification of solutions for guided waves in fluid-loaded viscoelastic composites with large numbers of layers, The Journal of the Acoustical Society of America 154 (2) (2023) 1073-1094. arXiv:https://pubs.aip.org/asa/jasa/article-pdf/154/2/1073/18094891/1073\_1\_10.0020584.pdf, doi:10.1121/10.0020584.
URL https://doi.org/10.1121/10.0020584</footnote>

<footnote n="7">⁷F. Treyssède, K. Nguyen, A.-S. Bonnet-BenDhia, C. Hazard, Finite element computation of trapped and leaky elastic waves in open stratified waveguides, Wave Motion 51 (7) (2014) 1093-1107. doi:10.1016/j.wavemoti.2014.05.003.</footnote>

<footnote n="8">⁸I. Bartoli, A. Marzani, F. Lanza Di Scalea, E. Viola, Modeling wave propagation in damped waveguides of arbitrary cross-section, Journal of Sound and Vibration 295 (3-5) (2006) 685-707. doi:10.1016/j.jsv.2006.01.021.</footnote>

<footnote n="9">⁹A. Marzani, E. Viola, I. Bartoli, F. Lanza Di Scalea, P. Rizzo, A semi-analytical finite element formulation for modeling stress wave propagation in axisymmetric damped waveguides, Journal of Sound and Vibration 318 (3) (2008) 488-505. doi:10.1016/j.jsv.2008.04.028.</footnote>

<footnote n="10">¹⁰J. L. Rose, Ultrasonic Guided Waves in Solid Media, 1st Edition, Cambridge University Press, 2014. doi:10.1017/CB09781107273610.</footnote>

<footnote n="11">¹¹C. Canuto (Ed.), Spectral Methods: Fundamentals in Single Domains; with 19 Tables, Scientific Computation, Springer, Berlin Heidelberg New York, 2006.</footnote>

<footnote n="12">¹²L. N. Trefethen, Spectral Methods in MATLAB, Society for Industrial and Applied Mathematics, 2000. doi:10.1137/1.9780898719598.</footnote>

<footnote n="13">¹³M. Lowe, Matrix techniques for modeling ultrasonic waves in multilayered media, IEEE Trans. Ultrason., Ferroelect., Freq. Contr. 42 (4) (1995) 525-542. doi:10.1109/58.393096.</footnote>

<footnote n="14">¹⁴N. A. Haskell, The dispersion of surface waves on multilayered media*, Bulletin of the Seismological Society of America 43 (1) (1953) 17-34. doi:10.1785/BSSA0430010017.</footnote>

<footnote n="15">¹⁵V. Pagneux, A. Maurel, Determination of Lamb mode eigenvalues, J. Acoust. Soc. Am. 110 (3) (2001) 1307-1314. doi:10.1121/1.1391248.</footnote>

<footnote n="16">¹⁶M. Gallezot, F. Treyssède, L. Laguerre, Contribution of leaky modes in the modal analysis of unbounded problems with perfectly matched layers, The Journal of the Acoustical Society of America 141 (1) (2017) EL16-EL21. doi:10.1121/1.4973313.</footnote>

<footnote n="17">¹⁷E. Georgiades, M. J. S. Lowe, R. V. Craster, Leaky wave characterisation using spectral methods, J. Acoust. Soc. Am. 152 (3) (2022) 1487-1497. doi:10.1121/10.0013897.</footnote>

<footnote n="18">¹⁸D. A. Kiefer, M. Ponschab, S. J. Rupitsch, M. Mayle, Calculating the full leaky Lamb wave spectrum with exact fluid interaction, J. Acoust. Soc. Am. 145 (6) (2019) 3341-3350. doi:10.1121/1.5109399.</footnote>

<footnote n="19">¹⁹M. A. Biot, Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. I. Low-Frequency Range, J. Acoust. Soc. Am. 28 (2) (1956) 168-178. doi:10.1121/1.1908239.</footnote>

<footnote n="20">²⁰M. A. Biot, Theory of Propagation of Elastic Waves in a Fluid-Saturated Porous Solid. II. Higher Frequency Range, J. Acoust. Soc. Am. 28 (2) (1956) 179-191. doi:10.1121/1.1908241.</footnote>

<footnote n="21">²¹M. A. Biot, Mechanics of Deformation and Acoustic Propagation in Porous Media, J. Appl. Phys. 33 (4) (1962) 1482-1498. doi:10.1063/1.1728759.</footnote>

<footnote n="22">²²M. A. Biot, Generalized Theory of Acoustic Propagation in Porous Dissipative Media, J. Acoust. Soc. Am. 34 (9A) (1962) 1254-1264. doi:10.1121/1.1918315.</footnote>

<footnote n="23">²³D. L. Johnson, J. Koplik, R. Dashen, Theory of dynamic permeability and tortuosity in fluid-saturated porous media, J. Fluid Mech. 176 (-1) (1987) 379. doi:10.1017/S0022112087000727.</footnote>

<footnote n="24">²⁴Y. Champoux, J.-F. Allard, Dynamic tortuosity and bulk modulus in air-saturated porous media, J. Appl. Phys. 70 (4) (1991) 1975-1979. doi:10.1063/1.349482.</footnote>

<footnote n="25">²⁵D. Lafarge, P. Lemarinier, J. F. Allard, V. Tarnow, Dynamic compressibility of air in porous structures at audible frequencies, J. Acoust. Soc. Am. 102 (4) (1997) 1995-2006. doi:10.1121/1.419690.</footnote>

<footnote n="26">²⁶T. Weisser, J.-P. Groby, O. Dazel, F. Gaultier, E. Deckers, S. Futatsugi, L. Monteiro, Acoustic behavior of a rigidly backed poroelastic layer with periodic resonant inclusions by a multiple scattering approach, J. Acoust. Soc. Am. 139 (2) (2016) 617-629. doi:10.1121/1.4940669.</footnote>

<footnote n="27">²⁷L. Boeckx, P. Leclaire, P. Khurana, C. Glorieux, W. Lauriks, J. F. Allard, Investigation of the phase velocities of guided acoustic waves in soft porous layers, J. Acoust. Soc. Am. 117 (2) (2005) 545-554. doi:10.1121/1.1847848.</footnote>

<footnote n="28">²⁸G. Belloncle, H. Franklin, F. Luppé, J. Conoir, Normal modes of a poroelastic plate and their relation to the reflection and transmission coefficients, Ultrasonics 41 (3) (2003) 207-216. doi:10.1016/S0041-624X(03)00095-7.</footnote>

<footnote n="29">²⁹J. F. Allard, M. Henry, V. Gareton, Pseudo-surface waves above thin porous layers (L), The Journal of the Acoustical Society of America 116 (3) (2004) 1345-1347. doi:10.1121/1.1778904.</footnote>

<footnote n="30">³⁰J.-P. Groby, A. Wirgin, E. Ogam, Acoustic response of a periodic distribution of macroscopic inclusions within a rigid frame porous plate, Waves in Random and Complex Media 18 (3) (2008) 409-433. doi:10.1080/17455030802061300.</footnote>
<footnote n="31">³¹F. Seyfaddini, H. Nguyen-Xuan, V.-H. Nguyen, Semi-analytical IGA-based computation of wave dispersion in fluid-coupled anisotropic poroelastic plates, International Journal of Mechanical Sciences 212 (2021) 106830. doi:10.1016/j.ijmecsci.2021.106830.</footnote>

<footnote n="32">³²A. T. I. Adamou, R. V. Craster, Spectral methods for modelling guided waves in elastic media, J. Acoust. Soc. Am. 116 (3) (2004) 1524–1535. doi:10.1121/1.1777871.</footnote>

<footnote n="33">³³F. Karpfinger, B. Gurevich, A. Bakulin, Modeling of wave dispersion along cylindrical structures using the spectral method, J. Acoust. Soc. Am. 124 (2) (2008) 859–865. doi:10.1121/1.2940577.</footnote>

<footnote n="34">³⁴F. H. Quintanilla, M. J. S. Lowe, R. V. Craster, Modeling guided elastic waves in generally anisotropic media using a spectral collocation method, J. Acoust. Soc. Am. 137 (3) (2015) 1180–1194. doi:10.1121/1.4913777.</footnote>

<footnote n="35">³⁵A. Hood, Localizing the eigenvalues of matrix-valued functions: Analysis and applications, Ph.D. thesis, Cornell University Library (2017).</footnote>

<footnote n="36">³⁶T. Betcke, N. J. Higham, V. Mehrmann, C. Schröder, F. Tisseur, NLEVP: A Collection of Nonlinear Eigenvalue Problems, ACM Trans. Math. Softw. 39 (2) (2013) 1–28. doi:10.1145/2427023.2427024.</footnote>

<footnote n="37">³⁷I. Gohberg, P. Lancaster, L. Rodman, Matrix Polynomials, Society for Industrial and Applied Mathematics, 2009. doi:10.1137/1.9780898719024.</footnote>

<footnote n="38">³⁸H. Gravenkamp, B. Plestenjak, D. A. Kiefer, E. Jarlebring, Computation of leaky waves in layered structures coupled to unbounded media by exploiting multiparameter eigenvalue problems, arXiv preprint arXiv:2404.15277 (2024).</footnote>

<footnote n="39">³⁹J. A. Weideman, S. C. Reddy, A MATLAB differentiation matrix suite, ACM Trans. Math. Softw. 26 (4) (2000) 465–519. doi:10.1145/365723.365727.</footnote>

<footnote n="40">⁴⁰R. Baltensperger, M. R. Trummer, Spectral Differencing with a Twist, SIAM J. Sci. Comput. 24 (5) (2003) 1465–1487. doi:10.1137/S1064827501388182.</footnote>

<footnote n="41">⁴¹G. Strang, Introduction to Linear Algebra, SIAM, 2022.</footnote>

<footnote n="42">⁴²N. Atalla, R. Panneton, P. Debergue, A mixed displacement-pressure formulation for poroelastic materials, J. Acoust. Soc. Am. 104 (3) (1998) 1444–1452. doi:10.1121/1.424355.</footnote>

<footnote n="43">⁴³J.-F. Allard, N. Atalla, Propagation of Sound in Porous Media: Modelling Sound Absorbing Materials, 2nd Edition, Wiley, Hoboken, N.J, 2009.</footnote>

<footnote n="44">⁴⁴J. D. Achenbach, H. A. Lauwerier, W. T. Koiter, Wave Propagation in Elastic Solids: North-Holland Series in Applied Mathematics and Mechanics, Elsevier Science, Amsterdam, 2014.</footnote>

<footnote n="45">⁴⁵P. Debergue, R. Panneton, N. Atalla, Boundary conditions for the weak formulation of the mixed (u,p) poroelasticity problem, J. Acoust. Soc. Am. 106 (5) (1999) 2383–2390. doi:10.1121/1.428075.</footnote>

<footnote n="46">⁴⁶L. Brillouin, Wave Propagation in Periodic Structures, MGH, 1946.</footnote>

<footnote n="47">⁴⁷A. Bernard, M. J. S. Lowe, M. Deschamps, Guided waves energy velocity in absorbing and non-absorbing plates, J. Acoust. Soc. Am. 110 (1) (2001) 186–196. doi:10.1121/1.1375845.</footnote>

<footnote n="48">⁴⁸F. Simonetti, M. J. S. Lowe, On the meaning of Lamb mode nonpropagating branches, J. Acoust. Soc. Am. 118 (1) (2005) 186–192. doi:10.1121/1.1938528.</footnote>

<footnote n="49">⁴⁹M. Castaings, B. Hosten, Guided waves propagating in sandwich structures made of anisotropic, viscoelastic, composite materials, The Journal of the Acoustical Society of America 113 (5) (2003) 2622–2634. doi:10.1121/1.1562913.</footnote>

<footnote n="50">⁵⁰J. M. Carcione, Wave propagation in anisotropic, saturated porous media: Plane-wave theory and numerical simulation, J. Acoust. Soc. Am. 99 (5) (1996) 2655–2666. doi:10.1121/1.414809.</footnote>

<footnote n="51">⁵¹A. Geslain, S. Raetz, M. Hiraiwa, M. Abi Ghanem, S. P. Wallen, A. Khanolkar, N. Boechler, J. Laurent, C. Prada, A. Duclos, P. Leclaire, J.-P. Groby, Spatial Laplace transform for complex wavenumber recovery and its application to the analysis of attenuation in acoustic systems, J. Appl. Phys. 120 (13) (2016) 135107. doi:10.1063/1.4963827.</footnote>

<footnote n="52">⁵²O. Dazel, F. Sgard, F.-X. Becot, N. Atalla, Expressions of dissipated powers and stored energies in poroelastic media modeled by (u,U) and (u,P) formulations, J. Acoust. Soc. Am. 123 (4) (2008). doi:10.1121/1.2874520.</footnote>

<footnote n="53">⁵³D. Royer, T. Valier-Brasier, Ondes élastiques dans les solides, Collecion Ondes, ISTE éditions, London, 2021.</footnote>