Efﬁcient and Accurate Modeling of the Surface Deﬂection of Thin Layers on Composite Substrates with Applications to Piezoelectric Parameter Measurements

: The electrical and mechanical response of multilayered structures involving a piezoelectric layer and bull’s eye shaped electrodes is investigated. A boundary element method is employed based on spectral domain Green’s functions. With this method, the electric ﬁeld distribution is determined ﬁrst, and the local mechanical displacement in a second step. As will be shown, this allows us to exploit cylindrical symmetry for the electric surface charge distribution, but not for the vertical surface displacements. The effect of substrate bending due to in plane-stress, introduced by the piezoelectric constant d 31 , and the beneﬁts of using bull’s eye electrode geometries with thick metallic backplates intended to reduce this effect are studied. A rigorous analysis and a largely simpliﬁed, but accurate approximation are compared. The application of this technique is demonstrated on a practical example for highly efﬁcient and accurate determination of selected piezoelectric coefﬁcients from surface topography measurements on such structures.


Introduction
Laser optical measurements of the vertical displacement profile on grown piezoelectric films on silicon substrates is a common technique to determine the piezoelectric parameters of the deposited piezoelectric layer [1,2]. However, the relation between piezoelectric parameters and the measurement is not straightforward, so that a common approach is to gradually adjust the parameters of a finite element model to achieve agreement between simulation and measurement [3]. The cylindrical symmetry of circular-or bull's eye-shaped electrode geometries on multilayered piezoelectric composites (see Figure 1) is commonly exploited to achieve numerically efficient finite element method (FEM) models of these structures. However, the justification for this simplification is not guaranteed by the geometry alone, but instead, depends also on whether the influence of anisotropy of the materials involved is negligible or not. As will be shown, the lateral anisotropy has a significant impact on the considered problem, and therefore a reduction to a symmetric representation in order to reduce complexity is not justified in general.
In this work, an efficient and accurate boundary integral technique based on spectral Green's functions (i.e., wavenumber domain representations) is employed to model the transverse deflection of abovementioned structures. A similar approach has been extensively used in the past for modeling wave propagation in surface acoustic wave sensors [4][5][6][7].
Unfortunately, the implementation of the proposed method in its rigorous form requires a relatively high level of modeling effort, which the practitioner, who is mainly interested in the results, usually wants to avoid. Fortunately, the characteristic properties of such structures allow for a significant reduction in complexity such that an accurate and fast calculation can be obtained, which may also appeal to an end user. The simplified calculation has been reported in the authors' earlier work [8] showing that the problem can be reduced to the determination of a spectral Green's function and numerical integration. In this work, the rigorous calculation and the subsequent simplifications leading to the simplified form are shown. Researchers who are mainly interested in the application can therefore may skip the rigorous treatment in Section 2 and jump directly to the reduced form in Section 5.

Basic Modeling Approach
The modeling approach in this work uses a boundary element method (BEM) implemented in the wavenumber domain based on Green's functions (GF) that is similar to the method discussed in [9]. It can be summarized as follows: we use two Green's functions G ϕ and G u , where the indices ϕ and u refer the electric potential and the mechanical displacement, respectively. G ϕ relates electrode charge distribution and potential and is used to determine the electric charge for a prescribed electrode potential. The mechanical displacement is calculated based on the charge distribution using G u . Interesting properties of the problem can readily be investigated by closer inspection of the Green's functions (see Section 2.1 below). For instance, whether a reduction to circular symmetry is permissible or not is most easily determined. Furthermore, the decay lengths of the fields, and thus the influence of the vertical extension of the modeling domain, e.g., the air half-space above the electrodes, and the material thicknesses can be estimated and may provide useful information and benchmarks for the developers of FEM simulation models as well.
The Green's functions used are specified in the spectral domain, to take into account the time-harmonic measurement mode of the laser Doppler vibrometer used to record the deflection. The method is therefore inherently dynamic, and it can be directly observed in the GF at which oscillation frequencies wave propagation phenomena will occur. Dispersion relations are shown in Section 2.1 revealing that Lamb-type waves exist that have no cut-off frequency, and must be suppressed by limiting the radial extension of the structure clearly below the critical wavelength, in order to avoid interference. A closer examination of the GF in the next section also discloses if the oscillation frequency in an experiment considered quasistatic is sufficiently low, or if a static model must be replaced by a dynamic one.

Symmetry Considerations
Some of the basic features of Green's functions will be discussed before they are treated in detail in Section 3.1. The symmetry properties of the GF indicated in Figure 2 will be used to simplify the calculations in advance. For laterally (x 1 , x 2 ) infinitely extended and homogeneous layers, the GF are shift invariant in the x 1 /x 2 plane, i.e., G(r, r ) = G(∆r), with ∆r = r − r and r = x 1 e x1 + x 2 e x2 . Here, r, r , and e i are the vector to the field point, the vector to the source point, and the unity vectors in i-direction, respectively. For cylindrical symmetry, G(∆r) = G[|∆r| cos (φ)e x1 + |∆r| sin (φ)e x2 ] is required, or equivalentlyG(k) = G[|k| cos(ψ)e k1 + |k| sin(ψ)e k2 ] in the spectral domain, where φ and ψ are the azimuthal angle in the spatial and wave domain, respectively. In the remainder of this paper, spectral domain quantities are marked with a tilde symbol˜. The calculation of the GF is performed in the wavenumber (k 1 , k 2 ) domain, with details given in [5,9,10], for instance, where k is the wave vector and k i its components in e ki directions, respectively. Results are shown for electrostaticG ϕ (k) and the mechanical GFG u (k). The key results are that anisotropy, boundary conditions, and layer thicknesses influenceG ϕ andG u differently, allowing certain simplifications to be applied. These results may also be beneficial for developers of finite element models in order to check for reduction in complexity without impairing validity or accuracy. The spectral electrostatic GFG ϕ (k) is shown in Figure 2 (left) , where only negligible dependence on the angle ψ is observed. On the contrary, a noticeable change over ψ is visible for the electromechanical GFG u (k). Therefore, a simplification to radial symmetry is only valid for the calculation of the electrical fields, but not for the mechanical displacement.

Dispersion Relations
The vertical displacement of the in the x 1 /x 2 -plane infinitely extended layer composite due to a surface charge distributionρ at its top is given bỹ in the spectral domain. At certain critical wavevectors k * (ω),G u [k * (ω)] shows poles indicating that a wave propagating in k direction is possible, even for vanishing electrical excitation. The associated critical wavelengths λ * and phase velocities v * are given by The dispersions of critical wavelengths for different aluminum backplane thicknesses are shown in Figure 3. It can be observed that at frequencies as low as 1 Hz, the critical wavelengths are very long and in the order of meters for the anti-symmetric Lamb wave. In case of high frequency operation and larger extended substrates, interference due to standing waves has to be considered. Although the dispersion is dominated by the aluminum plate, it varies considerably from the well known dispersion diagram of single aluminum plates (see, e.g., [11]), which can be attributed to the additional layers. The symmetric mode with its extraordinarily long wavelength (λ ≈ 5.8 m at f = 1 kHz) is not considered, as the limited extension of the structure does not support these waves. Therefore, fromG u , it can already be determined if standing wave interference will be of concern in the measurement setup.

Boundary Element Method
Defined piecewise constant electric potentials of ϕ(r) of 1 V and −1 V are applied to the central and outer ring electrodes, respectively, while the aluminum backplane is grounded (0 V). The resulting electric charge distribution induces mechanical displacements of the piezoelectric AlN film and the non-piezoelectric underlying layers. In the first step, the resulting surface charge distribution on the electrodes ρ(r) is determined by the method of mean weighted residuals. In a second step, the displacement field in out-of-plane direction u z (r) is calculated.

Method of Mean Weighted Residuals (MWR)
The method of MWR can be summarized as following [9]: the potentialφ(k) and the surface charge densityρ(k) are directly related by the spectral electrostatic Green's functioñ with k representing the radial wavenumber k 2 = k 2 1 + k 2 2 for the radially symmetric fields. The transformation to the spatial domain (with spatial radial coordinate r) is obtained by inverse Hankel transform (A scaling factor of 1/(2π) has been introduced for agreement with the two dimensional inverse Fourier transform of radially symmetric functions).
where J n (·) denote in the following Bessel functions of the first kind and nth order. The charge distribution is synthesized by an infinite series of trial functions B i (r) scaled by factors a i in the spatial and spectral domain. The residuals are weighted with functions T j (r) · r according to For convenience, such integrals can be shortened as an inner product with a distinction between spatial an spectral integrals according to From (6)-(9), a linear system or in matrix notation for the unknown coefficients a i of the series expansion (5) arises.

Displacement Field Calculation
Once the coefficients a i are calculated, the charge distribution follows from (5) in both domains. The vertical displacement u z can subsequently be derived from [7] whereG u (k, ψ) denotes a spectral Green's function relating mechanical displacement and surface charge density, which is in contrast toG ϕ (k) not symmetric as it varies with the angle ψ as previously illustrated in Figure 2.

Application
The basic approach outlined in Section 2 is now applied to the particular problem of the layered structure depicted in Figure 1. The respective material parameters of the layers used for calculation are summarized in Appendix A.1. The explicit derivation of the Green's function is, e.g., collected in [12,13]. First, the construction of the Green's function and an approximation forG ϕ based on a non-piezoelectric material are shown. The integrals in (8) and (9) are computed using a finite number of rectangular basis and weighting functions for which closed-form solutions are available. Therefore, very thin elements, for which the integrals converge only very slowly, can be calculated in closed form, such that the elements can be refined at the electrode edges where gradients are largest. The purpose of this rather cumbersome approach is to have an exact solution at hand to validate the applicability of simplifications. For example, it will be shown that the calculation of the charge distribution using the MWR approach can be avoided at all for the practical thicknesses of piezoelectric AlN layers typically in the range of a few hundreds on nanometers. A very efficient calculation requiring only the evaluation of one Green's functions and direct numerical integration can be obtained.

Electrostatic Green's Function
The crucial part of this method is the accurate knowledge of the Green's functions. Unfortunately, a closed form expression for the electrostatic Green's function considering all material layers yields a too complex function for direct analytical integration. Fortunately, the high conductivity of the comparably thick doped silicon layer leads to a strong concentration of the electrical field within the piezoelectric layer and thus only negligible electrical influence of the layer composition below the piezoelectric layer results. However, the dominant electric field in thickness direction and the non-zero piezoelectric stress constant e 31 , e 32 , or equivalently non-zero piezoelectric constants d 31 , d 32 produce shear stresses on the silicon surface and thus influenceG ϕ (k) noticeably. Our approach is therefore to consider the AlN layer alone and to neglect mechanical loading effects for the moment and account for it later by adjusting the parameters of the approximate model.
According to the IEEE standard on piezoelectricity [14], the two equivalent sets of constitutive piezoelectric equations in Voigt notation [15] are the stress-charge form and the strain-charge form They can be converted to each other using the relations Here, T, S, c, s, ε, E and D, are the tensors of stress, strain, elastic stiffness, compliance, permittivity the vectors of the electric field strength and the electric displacement, respectively. In addition, d and e refer to as the piezoelectric tensor in strain-charge and stress-charge form, respectively. Conversion between both systems is regularly required because the calculation of Green's function, e.g., in [12] is performed using (13), whereas piezoelectric materials are mostly characterized by the parameters of (14). The material tensors of AlN (wurtzite, class 6 mm, see, e.g., pp. 300 in [16]) are given by (Note that the structures of the material tensors in system (13) are equal to (14) for this class).
The representation (14) is more suitable in the stress free case T = 0, andG ϕ can be calculated from D = −ε T · ∇ϕ and ∇ · D = 0 with respective boundary conditions, yielding This result is very similar to the GF for isotropic materials as can be found, e.g., in [17]. The function f (h I , k) accounts for the extension of the air domain (I) with different boundary conditions at h I Infinite air halfspaces are easily implemented by the BEM. However, for FEM models, where the air domain has to be discretized in a vertical direction, the required extension is of interest and can be approximated with the following consideration. For a disk of radius R, the lowest dominant wavenumber is k ≈ 2π/(4R) and the hyperbolic functions are approximately 1 at h I k = π. Therefore, it is advised to use h I > 2R.
The stress free solution in (21) is exact for d = 0 pm/V (or e = 0 C/m 2 ). Deviations due to piezoelectricity can be well approximated by adjusting the ε T 11 and ε T 33 parameters to a certain extent as shown in It can therefore be concluded that approximations in closed form are possible. For strong piezoelectric coupling, the agreement suffers at low wavenumbers, but the asymptotic behavior as k → ∞ is always accurately reproduced, which is a crucial property as it allows an efficient implementation of the BEM in (6), which is discussed in the following section.

Calculation of the Electric Field Using MWR
The field calculation in cylindrical symmetry always requires the calculation of integrals with strongly oscillating kernels. As shown in Figure 4,G ϕ decays slowly (dominated by k −1 according to (21)) such that a numerical computation of the integral is exhaustive. In order to achieve a numerically efficient implementation, our approach is to use simple trialB i (k) and weighting functionsT j (k), which allow us to use closed form expressions for the infinite integrals. An alternative approach is to use specialized orthogonal functions that already approximate the expected charge distribution reasonably, such that a lower number of trial functions yield good results, but involves complicated integrals. For instance, this approach is used by [18] to formulate the problem as dual integral equations and solved using techniques, e.g., found in [19]. An interesting feature of our presented approach is that when the trial functions of [18] are used and the weighting functionsT j (k) are set equal to the trial functionsB i (k), then (6) yields the same integral and solution system as presented in [18]. The approach presented in this paper is therefore equivalent to the dual integral approach.  In the following, we use the first method with simple non-overlapping hollow cylindrical trial and weighting functions centered at the radii r i,j with wall thickness w i,j , shown in Figure 5 for which the integrals can be computed in closed form. The spectral domain representations are obtained by Hankel transform bỹ with G ϕ can also be expanded into an infinite series of exponentials with The right-hand-side of (6) can be represented by a sum of four so-called Lipschitz-Hankel integrals I(µ, ν; λ) defined by (A1) of the type type I(1, According to page 389 in [21] or page 315 in [22], the improper integral on the left side of (30) with highly oscillating kernel can be replaced by the finite integral over a smooth function stated in the right side of (30) that can be efficiently evaluated. Moreover, this integral can be expressed in closed form by means of elliptic integrals. To the best of the authors' knowledge, it is not straight-forward to find this representation in the common literature on integrals. It is therefore derived in the Appendix. Efficient implementations for elliptic integrals are reported in [23] and are made available for MATLAB, e.g., by T. Hoffend [24]. The sum in (29) converges for b > 2, which is always given. The amount of summands required to achieve a relative accuracy of f can be estimated by N = log( f )/ log(a 2 − 1) and is larger for piezoelectric materials with high permittivity (see (27)). By using the convergence acceleration method for alternating series of Longman [25], N can be drastically reduced while retaining high numerical accuracy. For the calculation of the expansion coefficients a i , the matrix M in (10) has to be evaluated for all j and i first, using the right-hand-side integral in (6). The charge distribution is calculated using a truncated version of (5) The symmetric matrix M summarizes the various interactions between potentials at segments r j and surface charge at r i . The mutual influence of the elements decreases with distance |r j − r i |, i.e., the self action contributions at the diagonal are dominant. For piezoelectric layers being thin compared to the element width w i,j , this decrease is stronger, and it is sufficient to calculate only a limited amount of upper diagonals of M and use the symmetry property. In Figures 6 and 7, the resulting surface charge densities and the respective electrical potentials for different thicknesses h II of 0.5 µm, 5 µm, and 50 µm of the piezoelectric layers are shown for the inner and outer electrode. It is observed that the charge distribution is nearly uniformly distributed and the charge accumulations at the electrode edges are not dominant for the 500 nm thick AlN layer used in the later shown experiment. This picture changes for thicker layers, especially when the thickness is comparable to the radii of the electrodes. Figure 6 also shows the simple approximation for constant charge densities and the simulated (C i , C o ) and approximated (Ĉ o ,Ĉ i ) capacitances of the inner and outer electrodes calculated by This further substantiates the validity of approximating the charge densities as constant over each electrode, when applying thin AlN layers.

Calculation of the Displacement Field
The displacement field can be calculated using (12), but the computation of the integral is inconvenient. Expanding the ψ-dependence into a Fourier series (G u Using the identity 2π 0 e −jkr cos(φ−ψ) e jnψ dψ = 2πe jn(φ−π/2) J n (kr) valid for integer n, the integral in (12) can be replaced by In Figure 8,G u (k, ψ) over ψ is shown for some wavenumbers k. It can be observed that the variation over ψ can be approximated by a single cosine function with period 4, and therefore a good approximation is found with (Note thatG u,4 =G u,−4 and J −n (x) = (−1) n J n (x) for integer n).G u (k, ψ) =G u,0 (k) + 2G u,4 (k) cos(4ψ) . (37) The approximate mechanical displacement field follows from i (k)k(G u,0 (k)J 0 (kr) The evaluation of this integral is still difficult because of the oscillatory behavior of the integrand. Again, an approximation for the asymptotic behavior ofG u that can be integrated is crucial for efficient numerical evaluation. In Figure 9,G u is shown at low (left) and high (right) wavenumbers. At low wavenumbers,G u depends strongly on the thickness of the aluminum layer (IV). On the contrary, at high wavenumbers, this dependence is not noticeable and an approximation using damped exponential functions is possible, where δ i denotes the damping factor. The assocatated parameters can be determined using Prony's method [26]. For example, the function in Figure 9 (right) can be well approximated by using only the three terms (i ∈ − 1, 0, 1) of (39), i.e., which is indicated in Figure 9 by the dashed line. By using this approximation and the trial functions in (23) for the calculation of the vertical displacement, integrals of type I(0, 1, −1) of the Lipschitz-Hankel integrals (see (A1) in the Appendix) are encountered. Their conversion to well behaved functions for efficient numerical evaluation is stated in the Appendix. The functionG u,4 is shown in Figure 10. The spectral decay rate is in the order of k −12 , i.e., for this function, integrals of I(1, 4, −12) must be evaluated. A conversion to smooth integrals according to [22] is only possible for I(1, 4, λ) with λ > −4, however, the original integral itself is already smooth enough for direct numerical integration using either the Romberg rule or by considering the method proposed in [27].

Results
The strong influence of d 13 on surface bending is best demonstrated by setting d 13 = 0 first. The vertical surface deflection is calculated using (38) and is shown in Figure 11 for different thicknesses of the AlN layer. The total deflection u z is dominated by d 33 and can be approximated by for thin layers. Noticeable deviations occur only at thicker layers. No substrate bending occurs, and therefore the lateral extension of the structure is not critical. However, the picture changes if d 13 is active and significant surface bending due to in-plane stress is introduced (see Figure 12). The substrate bending of the structure is dominated by the behavior ofG u at low wavenumbers and depends strongly on the thickness of the aluminum plate if the plate is below a certain thickness, as can be observed in Figure 9. Above a thickness of h IV = 2 mm, no significant change results for thicker plates. The propagation of the Lamb waves does only interfere noticeably for structures, which are laterally wider than the critical wavelengths shown in Figure 3. For the calculation of u z (r), the integration ofG u,0 is split into two parts. For the asymptotic behavior at high wavenumbers, the approximation obtained by Prony's method with subsequent least squares optimization given in (40) is used. The residual part ofG u,0 decays rapidly and can be evaluated accurately by straightforward numerical integration.

The Influence of the Radial Boundary Condition
In order to study the influence of the lateral extension of the structure when d 13 is active, a zero normal displacement u z = 0 at certain radii R c (R c larger than the outer electrode, but lower than the critical wavelength λ * ) is enforced by using a dyadic Green function with additional surface forcef z influence functionsH u andH ϕ : The surface force densities are discretized using trial functions D n (r) with spectral representationD n (k) according tof The anisotropy of the structure is neglected in this case, and therefore the results are approximations which show the qualitative influence of fixed displacement at an outer ring. Instead of the system (10) given by b = M · a, the new system must be solved with additional unknown expansion coefficients c. The components of the additional matrix elements are given by The system (44) can be rearranged yielding the unknown coefficients a and c. An additional displacement contribution to the integral (38) as a result of the forces required to maintain zero displacement at the substrate edges results in the form of Figure 13 shows the effect of such boundary conditions. The dashed lines represent the required force distribution (100 rectangular elements each). For better visibility of the effects, the displacement at the disk center is subtracted for each simulation. It can be observed that the effect of the lateral boundary condition on the actual shape in the electrode region is hardly noticeable, but an offset depending on R c occurs. It is therefore concluded that the resulting displacement fields are not sensitive to the lateral boundary conditions. This simplifies the calculation considerably, especially when anisotropy is implemented. The effect of anisotropy for a pinned boundary condition at R c = 1.5 mm is shown in Figure 12 at an azimuth φ = 0 • .

Efficient Calculation for Thin Layers
The charge distributions in Figure 6 are nearly constant and their magnitude can be estimated using ε eff . At the electrode edges, charges accumulate, but their impact is small for thin AlN layers. It can therefore be assumed that the charge distribution can be specified as known as long as the edge charges have no significant impact on surface bending. The spatial influence of a circular charge singularity is shown in Figure 14 and shows that the charge singularities act in a very limited range on the surface bending and can therefore be neglected. The charge distribution can therefore be estimated in advance, and the calculation of the MWR can be reduced to the evaluation of the integral in (36) with one trial functions for inner and outer electrode each. A further simplification is achieved when the infinite integral over the wavenumber is replaced by a Fourier-Bessel series [28] where η m,n represents the mth zero of J n (x). Considering only the 0th and 4th order ofG u and assuming piecewise constant charge distributions, the calculation can be approximated by k n = k 0,n and k m = k 4,m resulting in In Figure 15 simulation results using the described method are compared to measurements reported in [3], where the same structure was analyzed. From the measured raw data, an offset of 1.4 pm was subtracted. This offset has shown to be random for repeated scans. In this simulation, the results are very sensitive to the d 31 component, and a small modification from 2.1 to 2.128 pm/V results in good agreement between measurement and simulation. For the rigorous computation, the electric charge distribution was determined using the mean weighted residual approach. However, as was demonstrated, the numerically costly charge calculation can be avoided by approximating it using the effective dielectric constant of the AlN layer. This is permissibly due to the comparably high electrical conductivity of the n-doped silicon wafer (≈2 mS/m), which confines the electrical fields to the thin AlN layer at low oscillation frequencies. It is furthermore shown, that the radial mechanical boundary condition (e.g., a clamping at a certain radius) has, apart from causing an offset in vertical direction, negligible influence on the vertical displacements in the electrode regions. 128 pm/V. With anisotropy considered, the measurements (u z,meas ) agree well with the rigorous (u z,rigor ) and the approximate (u z,appr ) computation.

Conclusions
The effect of the piezoelectric constant d 31 on the surface deflection and the insensitivity to the lateral boundary condition could be shown. In absence of d 31 and for thin layers up to some microns, the electric and the mechanical field distributions follow the one-dimensional results very closely, i. e., ρ ≈ ε eff /h II V and u z ≈ −d 33 V. The influence of d 31 introduces a lateral stress, which deflects the surface in an anisotropic way, although the electric fields maintain their rotational symmetry. Simulations with neglected anisotropy may show noticeable deviations, especially at the outer electrode of the bull's eye structure. The implementation of the anisotropy is relatively simple, due to the strong spectral decrease of the Fourier coefficients in azimuthal direction. The use of the aluminum backplane to reduce the effect of surface bending has shown to be clearly beneficial. It could also be demonstrated by analyzing the Green's functions, that when a certain thickness is reached (in our case 2 mm), surface bending could not be reduced by a further increase of the backplane thickness h IV . From the Green's functions, the symmetric and anti-symmetric Lamb-type waves can also be determined. It is concluded that the antisymmetric mode may interfere with the measurement at frequencies above hundreds of Hertz and for thin backplanes.  αβ sin 2 (θ) For the calculation of the displacement fields, integrals like I(1, 0; −1) are required that can be represented by For I(1, 0; −1), a representation involving elliptic integrals exists and can also be found in [20].