Electroelastic guided wave dispersion in piezoelectric plates: Spectral methods and laser-ultrasound experiments

In a waveguide, waves are constrained to travel along a defined path due to boundaries that channel the energy.^1,2 Guided elastic waves have applications in non-destructive testing,³ material characterization,^4–7 and sensing.^8,9 If the waveguide’s material is piezoelectric, the mechanical and electrical fields are coupled, resulting in electroelastic waves.^1,10–13 These have found widespread applications in micro-electromechanical systems (MEMS), e.g., as surface acoustic wave (SAW) devices for signal processing and sensing.¹⁴ In other devices, such as bulk acoustic wave filters, spurious guided modes need to be avoided.¹⁵ Due to technological advances in materials, thin film technology, and further development of MEMS devices, application-driven research on guided waves in piezoelectric plates is still very active today.^16–19 

Electroelastic waves in piezoelectric media are well studied. The theoretical fundamentals are covered in classic textbooks.^1,2,10 The initial work describes the propagation of bulk waves.^12,13 Work on modeling guided waves in piezoelectric plates was pioneered by Tiersten,²⁰ who considered metallized surfaces while restricting propagation to two specific directions. Later, fundamental modes propagating along principal directions with open (unmetallized) and shorted (metallized) surfaces were studied by Joshi and Jin²¹ and validated experimentally with interdigital transducers. Syngellakis and Lee²² extended the modeling to arbitrary anisotropy and propagation direction and computed the dispersion of electrically shorted higher-order modes. Electrically open modes at different propagation directions were computed by Kuznetsova et al.²³ Sun et al.²⁴ presented solutions for open–open, open–shorted, and shorted–shorted surfaces. A periodic grating of open and shorted regions can be used to design phononic crystals with tunable bandgaps.²⁵ Guo et al.²⁶ analyzed multilayered piezoelectric plates based on a state-space formalism. Last, Zhu et al.²⁷ accounted for the electrical losses in the metallized surfaces of plates.

Root-searching of the characteristic equations was used in all previously mentioned studies (except the phononic crystal) to compute dispersion curves. Although this method is known to be slow and to miss solutions, it is still prevalent today. A software based on root-searching was made available in 1990 by Adler²⁸ and is still actively being used.^29,30 An alternative treatment of the transcendental characteristic equations is by analytic approximations for limiting cases.³¹ 

In contrast to root-searching of the characteristic equations, semi-analytical methods discretize the waveguide cross section and then solve a standard algebraic eigenvalue problem to reliably and efficiently obtain solutions. Various discretization methods have been used for semi-analytical models of piezoelectric waveguides. Finite elements were used very early by Lagasse³² for ridge and wedge guides, while Cortes et al.³³ treated multilayered plates. Spectral methods are based on a global function basis to expand the field and are particularly well-suited for layered structures. In this context, Laguerre polynomials were used for surface waves^34,35 and Legendre polynomials³⁶ for plates. Laude et al.³⁷ used a Fourier series representation but encountered difficulties due to the non-vanishing electrical field outside the plate. Despite the advantages of semi-analytical methods, their widespread adoption is hindered by the lack of publicly available implementations and the relatively high level of knowledge on numerical computing required to implement them.

From an experimental point of view, laser-based ultrasound (LUS) has proven to be a valuable tool to precisely resolve the multi-modal dispersion curves of plates^5,38–41 and to access specific modes in the response spectrum.^4,7,42–45 Importantly, the contact-free operation of LUS allows for flexible scanning of the wave field under various environmental and boundary conditions. This is particularly valuable for model validation and for inverse characterization of material parameters. We are aware of two studies^46,47 that present LUS measurements of piezoelectric plates, but the effect of electrical BCs was not examined therein.

In this work, we present a semi-analytical model for guided electroelastic waves in piezoelectric plates based on spectral methods. Our focus lies on the highly instructive spectral collocation method (SCM).^48–50 The SCM was previously shown to be efficient and robust to compute dispersion of purely elastic waves.^51,52 Its main advantage is the straightforward implementation based on the underlying partial differential equations in their conventional strong form. The SCM code by Kiefer⁵³ is extended to piezoelectric plates, and we make our implementation publicly available under the name GEW piezo plate.⁵⁴ As an alternative, highly related method, we also introduce the spectral element method (SEM)^55,56 in  Appendix A. While it is numerically superior to the SCM, it is at the same time conceptually more intricate. Our SEM implementation will be provided in GEWtool,⁵⁷ a software for guided wave analysis. Overall, we hope to demonstrate with this contribution the power of spectral methods, the simplicity of SCM, and contribute to a wider adoption of semi-analytical methods for dispersion computations. Furthermore, we validate our model with LUS experiments and contribute a compilation of multi-modal guided wave spectra in $LiNbO 3$ wafers, including two different propagation directions and three different electrical BCs.

This paper is organized in the following way: in Sec. II, we derive the guided wave problem in piezoelectric media, with focus on the treatment of electrical BCs. The implementation of the SCM is explained in Sec. III, and computational results are presented in Sec. IV. Our LUS system is described in Sec. V. Last, measurement results are presented and discussed in Sec. VI before we conclude in Sec. VII.

The computational method proceeds in three steps:

1. Re-write the governing equations (4) and the Neumann BCs (5) and (8) in matrix form.

2. Discretize both differential systems.

3. Incorporate the discrete BCs into the governing system. Optionally, replace the Neumann BCs with Dirichlet BCs.

Second, a large number of numerical discretization techniques could be used to obtain algebraic approximations to the boundary-value problem (10) and (12). For the numerically superior SEM, we refer to  Appendix A and instead focus on the more instructive SCM in the following. The SCM is particularly simple to implement as it does not require a weak formulation nor previous meshing of the geometry. Furthermore, it exhibits exponential convergence when increasing the degrees of freedom, leading to small matrices. The Chebyshev SCM approximates an unknown function $ϕ(z)$ by a weighted superposition of the first $N$ Chebyshev polynomials, and the error is then minimized on $N$ so-called Chebyshev–Gauß–Lobatto collocation points $z i$⁠, $i∈{1…N}$⁠. For details on the SCM, refer to the standard literature.^48,49 For implementation purposes, it is important that the method allows for the explicit construction of a differentiation matrix $D z$⁠. To explain its meaning, assume that $ϕ d=[ϕ( z i)]$ is the mathematical vector that collects the values of $ϕ(z)$ at all collocation points $z i$⁠, and $ϕ d ′$ is the vector collecting the corresponding values of $∂ zϕ(z)$⁠. The differentiation matrix relates these two vectors through $ϕ d ′≈ D z ϕ d$⁠, providing a discrete counterpart to the differentiation operation. We use the MATLAB package DMSUITE⁵⁰ by Weideman and Reddy to generate the differentiation matrices $D z$ and $D z 2$ of size $N×N$⁠.

Equation (19) represents an algebraic eigenvalue problem. By fixing $k= k 0$⁠, we may compute the eigenvalues $ω n 2$ with standard techniques, e.g., eig in MATLAB. We remark that $M$ is a singular matrix, leading to some restrictions on the choice of the eigenvalue solver. Nonetheless, modern methods appropriately handle this case. It is also possible to do the reverse and fix the angular frequency $ω= ω 0$ in order to compute the eigenvalues $i k n$⁠. This quadratic eigenvalue problem can be solved by standard companion linearization techniques as implemented in MATLAB’s polyeig. Note that even for real frequencies, this approach yields a complex-valued wavenumber spectrum;¹ see  Appendix D. In either case, the eigenvectors $Ψ d$ provide access to the mechanical and electrical field distributions $u( z i)$ and $ϕ( z i)$⁠.

Our computations show very good agreement with those published by Kuznetsova et al.;²³ see  Appendix B. Note that the higher the modal order (i.e., the higher the frequency), the larger the number of discretization points $N$ needs to be. To ensure six-digit accuracy with the SCM, we use $N=20$ in all subsequent computations. An extended convergence analysis is presented in  Appendix B. All numerical experiments are performed on one core of an Apple M1 Pro processor. Our SCM implementation uses the QZ-algorithm from MATLAB’s eig-function to compute a set of dispersion curves at 120 $k$-values in 0.7 s (6 ms per $k$-value). The SEM implementation in GEWtool computes the same example in 0.6 s (5 ms per $k$-value), which is slightly faster thanks to the symmetry of the SEM matrices. However, per default GEWtool exploits subspace methods (eigs) and sparse matrix representation. 25 modes are located in this way in only 0.25 s (2 ms per $k$-value).

The presented computational method is rather easily extended to multilayered plates.⁵² This is particularly interesting for MEMS devices, which often consist of layered structures. The layers can include piezoelectric, dielectric, and metallic media. Electrical currents in the latter can lead to losses, and this can be accounted for by modeling a complex effective permittivity. These extensions are concisely discussed in  Appendix C for the interested reader.

We now investigate guided waves in a $510μ m$-thick monocrystalline $128 °$ Y-cut lithium niobate (⁠ $LiNbO 3$⁠) wafer. This material is often used for MEMS as it exhibits large electromechanical coupling.^61,63 The material parameters are given in  Appendix E. To account for the crystal orientation in the computation, we need to rotate the tensors $c$⁠, $e$⁠, and $ϵ$—which are initially given⁶⁴ in the orthonormal $e X e Y e Z$-material coordinate system (Fig. 3)—to the wave propagation frame $e x e y e z$ (Fig. 1). Here, $e X$ lies within the wafer surface and is perpendicular to the wafer flat. The material’s $e Z$ axis (⁠ $e Y$ axis) has been rotated out of the plate’s normal $e z$ by an angle of $μ= 38 °$ (⁠ $μ= 128 °$⁠) around $e X$⁠.

Two directions of propagation are considered: (i) along $X$ and (ii) along $X+ 90 °$⁠, i.e., along the wafer flat; see Fig. 3. Furthermore, three different combinations of electrical BCs are systematically investigated: shorted–shorted, open–shorted, and open–open. The different BCs are compared in Fig. 2(a) for propagation along $X$ and in Fig. 2(b) for propagation in $X+ 90 °$⁠. Several regions of the dispersion curves are observed to be quite strongly affected by the BCs, especially for propagation in $X+ 90 °$⁠. Strong electromechanical coupling can be expected in these regions. The SCM computation also provides the mechanical and electrical modal field distributions. Three examples picked at the marks in Fig. 2(b) are depicted in Fig. 2(c). Thereby, the field was extruded in the propagation direction according to the harmonic ansatz from (3).

The samples investigated were double-side polished SAW-grade $128 °$ Y-cut $LiNbO 3$ wafers with a nominal thickness of $500μ m$ and an average roughness specified below 1 nm.⁶⁵ The shorted BCs were produced by depositing 3 nm titanium and 400 nm gold by physical vapor deposition carried out in a Balzers PLS 570, metallizing either one or both sides of the wafer.

We use a LUS system to record the spatiotemporal response of the plate. The LUS setup is sketched in Fig. 3(a). To generate waves, we use a pulsed laser (Bright Solutions Wedge HB 1064) with a wavelength of 1064 nm, a pulse duration of about 1.5 ns, pulse energies of approximately 2 mJ, and a pulse repetition rate set to 1000 Hz. The laser is focused by a cylindrical lens with 130 mm focal length and deflected toward the sample by $90 °$ with a dichroic mirror, which is mounted on an automated translation stage. By scanning the mirror, the distance between generation and detection is varied. We deliberately moved the cylindrical lens out of the focal position by about 3 mm with the intention to couple efficiently into modes of lower $k$ values and to avoid ablation. The cylindrical lens produces a line shaped thermo-elastic source, predominantly emitting plane waves with a wave vector perpendicular to the line.

The out-of-plane component of the resulting surface displacement $u z(x,t)$ is recorded with a two-wave-mixing vibrometer (Sound & Bright Tempo-FS200) operating at 532 nm wavelength connected to an oscilloscope (Teledyne LeCroy WaveRunner HRO 66Zi). The detection unit was focused onto the sample through the dichroic mirror with a lens of 100 mm focal length. A line of 14 mm was sampled along $x$ by scanning the dichroic mirror with a pitch of $50μ m$⁠. To improve the signal-to-noise ratio, 5000 traces were averaged for each measured position, and the bandwidth of the oscilloscope was limited to 50 MHz.

Figure 4(a) shows the recorded $u z(x,t)$ for the open–shorted case with propagation in the X-direction. Applying a spatiotemporal Fourier transform results in a representation of the wavefield in the $k$- $f$-domain⁶⁶ and is shown in Fig. 4(b), where we plot the spectral velocity magnitude $| v z(k,f)|=ω| u z(k,f)|$⁠.

To enhance the experimental dispersion curves, we apply a window to the signals in the $x$- $t$-domain. The window is designed to select contributions with energy velocities between 500 and 7000 m/s by cutting out a wedge of the data. To reduce sidelobes, the window $w(x,t)$ is smoothed by convolution with a Hanning window in the $x$ and $t$ dimensions. Figure 4(c) shows the windowed signals $u z(x,t)w(x,t)$ and also indicates the area selected by the limiting energy velocities as dashed lines. The resulting dispersion maps with reduced noise and artifacts around $k=0$ are shown in Fig. 4(d). A portion of the dispersion spectrum around 6 MHz has a negative wavenumber but positive energy velocity and is, thus, detected. Such backward waves⁴⁵ appear in plates next to zero-group-velocity points. For simplicity, for Fig. 5, we symmetrize the measured spectra by $v z(k,f)+ v z(−k,f)$ and only show positive values for $k$⁠.

The open–open samples are almost transparent to the 1064 nm generation laser, resulting in insufficient coupling into guided waves. We followed two separate strategies to obtain measurements in this case: First, we coated the samples with black paint, which in some attempts provided good results. Here, we reduced the generation laser energy to avoid ablation (we otherwise observed thickness resonances as the predominant contributions in the spectrum). To enhance detection on the black painted sample, we added a small dot of silver paint, which we focused the vibrometer on. Second, we used a 355 nm laser (Teem photonics PowerChip UV 355 nm) in order to exploit the significantly higher absorption of $LiNbO 3$ in the ultraviolet (UV) spectral range. Its pulse energy was $25μ J$⁠, the pulse length was about 1 ns, and the repetition rate was set to 250 Hz. The UV laser was guided to the sample by a flexible periscope with seven mirrors. To scan the laser along the sample, the periscope’s end-piece together with a 100 mm focal length spherical lens was mounted at a linear scanning stage and pointed at the sample under an angle of $≈ 45 °$⁠. This way, the primary setup described above could remain, while providing a workaround for the transparent sample. The workaround is shown in Fig. 3(b).

With these two strategies, we were able to produce dispersion maps of sufficient quality for comparison with the SCM. The best results are shown in Fig. 5 (UV laser in (c), black paint in (f)). We note that Yang and Tsai⁴⁶ successfully measured this case using a 532 nm laser and a high pulse energy of 200 mJ.

The experimental and theoretical results are superposed in Fig. 5. The computations with the SCM are shown as red dots on top of the spectral magnitude obtained by the LUS experiments. The insets indicate the direction of propagation and whether the surfaces are metallized or not. Shear-horizontal waves, which are included in the SCM result (and shown faintly in Fig. 2), are neither generated nor detected with the experimental setup and, thus, do not show. For clarity, these modes have been identified in the SCM solutions and omitted in the plot. For the calculations shown, we used measured thickness values of the samples obtained with a micrometer screw, ranging between $500μ m$ and $510μ m$⁠. For the coated wafers, we used the overall thickness, including the coating.

Based on Fig. 5, excellent agreement is found between theory and experiment. Small deviations are found in the highest order modes for the coated cases. These are well explained by the finite thickness of the metallization, which is not represented in the SCM solutions here, but in an additional solution in  Appendix C. Moreover, small experimental deviations from the theoretical directions of propagation cannot be excluded. We stress that the effect of different electrical BCs clearly shows in the experiment. For an assessment of the differences, refer to Fig. 2. The results obtained on the uncoated plates [Figs. 5(c) and 5(f)] exhibit a lower signal-to-noise ratio due to the workarounds described in Sec. V, but the acquired modes match well with the theory.

The configuration with the UV-laser results in an oval spot on the sample surface, rather than a line source. In this case [Fig. 5(c)], waves with arbitrary wave vector orientations are generated. All of them are potentially detected on the scan line, even when the wave vectors are not collinear to the scanning direction.⁶⁷ While the main features of the measured dispersion show good agreement with the calculations at fixed wave vector orientation, the blurred out regions in the experiment may be due to these wave skewing effects.

We presented the SCM to calculate dispersion curves and mode shapes for guided electroelastic waves in piezoelectric plates. LUS experiments were conducted to measure dispersion curves on a $LiNbO 3$ wafer. Comparison of theory and experiment for different directions and electrical BCs showed excellent agreement. From these results, material characterization based on an inverse problem seems feasible. A need for new characterization techniques is especially present in the case of deposited thin films and at high frequencies. The properties of thin films often differ from the bulk material and depend on the growing process and layer thickness. We consider the possibility to measure a sample under different electrical BCs particularly promising to disentangle elastic, piezoelectric, and electric material properties.

In the context of multilayered structures, the SCM could be used in combination with LUS measurements on different production steps of a stacked MEMS device, e.g., on the pure substrate of a SAW device and then on the deposited stack, including electrodes. Previous work has proven that LUS can indeed be used to measure guided waves and resonances in the GHz range for material characterization purposes.^44,68 With the optical spot sizes that are achievable in practice, effective generation and detection of guided waves down to $1μ m$ to $2μ m$ wavelength is feasible.

As an alternative to SCM, we discussed the SEM in  Appendix A. It is based on the weak form of the waveguide problem and is numerically superior to SCM but also somewhat more difficult to understand. An implementation of the SEM is available in GEWtool⁵⁷ for multilayered piezoelectric plates.

This research was funded in whole or in part by the Austrian Science Fund FWF (No. P 33764). This project is co-financed by research subsidies granted by the government of Upper Austria (No. Wi-2021-303205/13-Au). D. A. Kiefer has received support under the program “Investissements d’Avenir” launched by the French Government under Reference No. ANR-10-LABX-24. We thank Istvan Veres (Qorvo), Thomas Berer (Qorvo), and Jannis Bulling (BAM) for valuable discussions on the presented matter. For the purpose of open access, the author has applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.

The authors have no conflicts to disclose.

D. A. Kiefer: Formal analysis (lead); Funding acquisition (supporting); Methodology (equal); Software (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). G. Watzl: Data curation (lead); Investigation (equal); Writing – review & editing (supporting). K. Burgholzer: Resources (lead). M. Ryzy: Writing – review & editing (equal). C. Grünsteidl: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Software (equal); Supervision (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are openly available on Zenodo at https://doi.org/10.5281/zenodo.13828765, Ref. 69.

The spectral element method (SEM) is numerically superior to SCM because it (i) conserves the hermitian symmetry of the problem, (ii) preserves the regularity of the operators, and (iii) converges faster. These properties also facilitate the postprocessing, e.g., computation of group velocities⁷⁰ and zero-group-velocity points⁵⁹ as well as mode-sorting.⁷¹ SEM is akin to the finite element method but uses high-order polynomials defined on the entire domain as a basis; i.e., no meshing is necessary. In this section, we sketch the SEM-discretization of the waveguide problem without providing theoretical details, for which we rather refer to the literature.^56,72

A comparison of phase velocity $c p=ω/k$ dispersion curves calculated with the SCM (and the SEM from  Appendix A) to those published by Kuznetsova et al.²³ is shown in Fig. 6 and shows very good agreement. The latter authors employ root-finding of the characteristic equation. Two different crystal cuts of $LiNbO 3$ and propagation directions are depicted. Material constants are taken from Ref. 64 and rotated as needed before solving. Note that the root-searching method that serves as reference missed solutions, which is an inherent drawback of this method that is resolved by the semi-analytical methods proposed here.

The convergence of solutions is studied in Fig. 7 for an arbitrarily selected mode. The wavenumber is fixed at 20 rad/mm, and we compute the angular frequency $ω$ of the 12th mode (ordered in increasing $ω$⁠). As reference solution $ω 0$⁠, we use a high number of discretization points of $N=50$⁠. The convergence is compared to that of the classical Finite Element Method (FEM) with first order Lagrange polynomial basis functions (implemented in this case in GEWtool using a multilayered plate). While the FEM is known to converge algebraically, spectral methods converge exponentially as long as the solution is sufficiently smooth.^48,49 Hence, the relative error $|ω− ω 0|/ ω 0$ is seen to decrease rapidly for the spectral methods when increasing the discretization $N$⁠. Particularly impressive is the rapid convergence of the SEM implemented in GEWtool ( Appendix A). For this rather high-order mode, 12 accurate digits are attained already with $N≥26$ for the SCM and $N≥20$ for the SEM.

It is relatively straight forward to extend both the SCM and the SEM to multilayered plates. For the SCM, the approach explained in Ref. 52 can be extended to include electrical interface conditions (continuity of $ϕ ~$ and the normal component of $D ~$⁠). For the SEM, each layer is represented by one element and a conventional element assembly is performed,⁷² which is implemented in GEWtool.

The latter approach is used in the following to assess the influence of the gold (Au) coatings. We compare solutions ignoring the gold layers (⁠ $500.8μ m$ $LiNbO 3$ shorted–shorted) to those accounting for them (⁠ $0.4μ m$ Au– $500μ m$ $LiNbO 3$– $0.4μ m$ Au open–open). For gold, we use isotropic material properties (⁠ $C 11=111 GPa$⁠, $C 44=25 GPa$⁠, $ρ=19300 kg/ m 3$⁠, $ϵ= 6.9 ϵ 0$ and $σ=45.2 MS/ m$⁠). Gold is an excellent conductor exhibiting $σ/( ω maxϵ)$ $≈ 4.7× 10 9≫1.$ Hence, the resulting wavenumbers are almost real-valued, and these solutions are shown in Fig. 8. The gold layers tend to decrease the frequency due to the effect of mass loading,¹ which is especially visible for the higher-order modes. This explains the small mismatch observed between LUS and SCM above 20 MHz in Fig. 5(d).

A major advantage of semi-analytical methods is the straight-forward computation of the complex wavenumber spectrum. To this end, we fix real-valued frequencies in (19) and solve the quadratic eigenvalue problem in $ik$⁠, which is naturally complex-valued. An example is implemented in GEW piezo plate,⁵⁴ and the result is reproduced in Fig. 9. The companion linearization used to solve the polynomial eigenvalue problem doubles the problem size. For this reason, the 400 frequency points in this example compute in 23 s (57 ms per $ω$-point). A more efficient implementation is provided in GEWtool’s⁵⁷  computeK-function, which achieves the result in 4.4 s (11 ms per $ω$-point). Instead of relying on MATLAB’s polyeig-function, the companion linearization is solved with the subspace methods from MATLAB’s eigs-function, in this case for 30 modes.