Recursive Asymptotic Hybrid Matrix Method for Acoustic Waves in Multilayered Piezoelectric Media

This paper presents the recursive asymptotic hybrid matrix method for acoustic waves in multilayered piezoelectric media. The hybrid matrix method preserves the numerical stability and accuracy across large and small thicknesses. For discussion and comparison, the scattering matrix method is also presented in physics-based form and coherent form. The latter form resembles closely that of hybrid matrix method and helps to highlight their relationship and distinction. For both scattering and hybrid matrix methods, their formulations in terms of eigenwaves solution are provided concisely. Making use of the hybrid matrix, the recursive asymptotic method without eigenwaves solution is described and discussed. The method bypasses the intricacies of eigenvalue-eigenvector approach and requires only elementary matrix operations along with thinlayer asymptotic approximation. It can be used to determine Green’s function matrix readily and facilitates the trade-off between computation efficiency and accuracy.


Introduction
For many years there has been considerable interest in the study of acoustic wave propagation in multilayered piezoelectric media.Many techniques have been developed for analysis of such media, including transfer matrix method [1], impedance/stiffness matrix method [2][3][4], scattering/reflection matrix method [5][6][7] and hybrid matrix method [8].A comprehensive review of these methods has been provided in [9] along with their variants, numerical stability, computational efficiency, usefulness and deficiency.Since the transfer matrix method becomes unstable toward large thicknesses, while the impedance matrix method is inaccurate toward small thicknesses, they are not to be discussed further below.On the other hand, owing to their numerical stability and accuracy, both scattering and hybrid matrix methods deserve to be exploited further as mentioned or demonstrated in some recent works [10][11][12].In particular, the scattering matrix methods so far have been presented more in physics-based form (in terms of reflections and transmissions), which has motivated the unified matrix formalism in [10].However, with the unified formalism therein, it is still not clear about any relationship or distinction with hybrid matrix method.Moreover, most matrix methods thus far rely on the eigenwaves solution in their basic building blocks.Since eigensolver often takes substantial computations, it is useful to consider other methods without the need for eigenwaves solution [9].
In this paper, we present the recursive asymptotic hybrid matrix method for acoustic waves in multilayered piezoelectric media.The method is extended from the non-piezoelectric case [8] and exploits the hybrid matrix which preserves the numerical stability and accuracy across large and small thicknesses (instead of stiffness matrix [13] that may become inaccurate).For discussion and comparison, we also present the scattering matrix method in physics-based form and coherent form.The latter form resembles closely that of hybrid matrix method and helps to highlight their relationship and distinction.For both scattering and hybrid matrix methods, their formulations in terms of eigenwaves solution are provided concisely.Making use of the hybrid matrix, the recursive asymptotic method without eigenwaves solution is described and discussed.The method bypasses the intricacies of eigenvalue-eigenvector approach and requires only elementary matrix operations along with thin-layer asymptotic approximation.It can be used to determine Green's function matrix readily and facilitates the trade-off between computation efficiency and accuracy.

Acoustic Waves in Multilayered
Piezoelectric Media

Problem Formulation
Figure 1 shows a planar multilayered structure comprising N piezoelectric layers stratified along ẑ direction (within optional external layers 0 and N+1).For each layer f of thickness f h , its upper and lower interfaces/ boundaries are denoted by f Z  and f Z  , respectively.Let the fields in each layer f be described by field vector Assuming plane harmonic wave with exp( ) j t  time dependence and transverse wavenumber t t k s   , the field vector f f satisfies a first-order differential equation as where the layer system matrix is f A consists of the material parameters of layer f specified in terms of mass density f  and various stiffness constants, piezoelectric stress constants and permittivity via f  's (see [1]).

Solution with Eigenwaves
Equation (3) can be written as an eigenvalue problem ( ) whose solutions represent the eigenwaves within each layer f.For convenience, the normal wavenumbers zf k and their associated eigenvectors f  can be grouped into the following matrices: Here, ( ) . The superscripts '>' and '<' stand for "upward-bounded" and "downward-bounded" partitions, which correspond to upwardbounded and downward-bounded eigenwaves respectively (cf.boundedness/radiation condition).In line with field vector (1), the eigenwave matrix σ υ partitions.Each of these may be further partitioned in accordance with their compositions in (2) as , (Note that our convention in ( 7)-( 8) is that the notation without superscript '>' or '<' represents fields while the same notation with such superscript represents waves of upward-bounded or downward-bounded type.)Using the matrices above, the field vector solution can be expressed as f c is the coefficient vector (to be determined), while ( ) ( ) is the wave amplitude vector that lumps the exponential terms together.Following the upward-bounded and downward-bounded associations above, these vectors can be partitioned into ( ) , ( ) ( ) (and their decompositions) are functions of z, while f ψ , f c (and their decompositions) are not.Furthermore, the field vector f f is continuous across the interface of two different layers, so we have . However, the wave amplitude vector . Thus, it is important to specify exactly the z location of the interface to be within which of the two adjacent layers.

Scattering Matrix Method
Using the eigenwaves in each layer f, one can proceed to determine the solution for a stack of multilayered media.
To that end, we first define the local interface scattering matrix that better describes the physics of wave scattering (reflection/transmission) at the interface of layers f and f + 1: , 1 denote the local reflection and transmission matrices for waves incident from layer f to f + 1, while 1, denote those for incidence from layer f +1 to f.These matrices can be derived directly in terms of the eigenwaves of both layers as Based on the local interface scattering matrix, one can determine the scattering matrix for additional layers (one at a time) of a stack using certain recursive algorithm.In particular, consider the downward-bounded waves incident from layer f +1 toward layer 0. The stack reflection and transmission matrices, , can be obtained from the local interface scattering matrix and the preceding ,0 f r and ,0 f t using the recursive algorithm (cf.( 23) and (25) of [9]): Likewise, the stack reflection and transmission matrices 0, 1 f  r and 0, 1 f  t for incidence of upward-bounded waves from layer 0 toward layer f +1 can also determined via recursive algorithm: 0, 1 0, ,0 The form of ( 13)-( 16) facilitates the physics-based description of wave multiple reflections in the stack of multilayered media.
As an alternative, it is instructive to define the matrix relating the wave amplitude vectors in the form Such matrix has been denoted as layer-interface scatterer [9], since it combines the layer scatterers ( ) with interface scattering matrix as To be in coherent form, the stack scattering matrix [1: ] f S is also defined in place of , r t , which embeds (within layers 0 and f+1) the stack from layer 1 to f (denoted by the superscript [1:f]), i.e.
In essence, they represent a full matrix variant of algorithm A3 in Table 1 of [9].

Hybrid Matrix Method
The scattering matrix method in the previous section involves relations among wave amplitude vectors f  w and f  w .As mentioned earlier, since these vectors are not continuous across interfaces, it should be more convenient to work directly with field variables instead.In this aspect, a variety of definitions and algorithms are possible including the transfer and impedance matrix methods.These methods are not unconditionally stable since they may cause numerical instability or inaccuracy problem for very large or very small layer thickness.Such problem can be overcome altogether by resorting to f H is called layer hybrid matrix since it has a mixture of impedance, admittance and transfer elements.Using the eigenwaves in each layer, the layer hybrid matrix can be determined as 11 12 21 22 It can be analytically shown that f H is still numerically stable even when the layer thickness tends to infinity or zero.Indeed, assuming at least slight loss as in practice, when the layer thickness tends to infinity, i.e.
On the other hand, when the layer thickness tends to zero, it is evident that ( ) ( ) , thus the hybrid matrix in (26) becomes 11 12 Therefore, the hybrid matrix preserves the numerical stability and accuracy across large and small thicknesses.
For multilayered media, the stack hybrid matrix [1: ] f H for a stack from layer 1 to f (denoted by the superscript This stack matrix can be obtained by incorporating the layer matrix f H into the recursive algorithm: ( ) ( ) ( ) ( ) Notice that the form of ( 30)-( 33) resembles closely that of ( 21)-( 24), which helps to highlight their relationship and distinction.In particular, both scattering matrix and hybrid matrix do not differ much in their recursive algorithms for a stack of multilayered media.However, besides relating different entities (waves f w vs. fields f f ), their basic matrices are distinct, i.e., f S involves eigenwaves of two layers in ( 12) and (18); while f H involves eigenwaves of individual layer only in (26).

Solution without Eigenwaves-Recursive Asymptotic Method
Thus far both scattering and hybrid matrix methods rely on the eigenwaves (of two or one layers) as the input (for f S and f H ) in each recursion to arrive at [1: ] f S and [1: ] f H .For such eigenwaves solution, there exist various intricacies of solving the eigenvalues and eigenvectors including complex root searching, degeneracy treatment and upward/downward eigenvector sorting or selection.
To obviate the need for eigenwaves, we resort to the recursive asymptotic hybrid matrix method.The method bypasses the intricacies of eigenvalue-eigenvector approach and requires only elementary matrix operations along with thin-layer asymptotic approximation as described below.
For each individual layer f, we geometrically subdivide the layer into n+1 sublayers having thicknesses as shown in Figure 2.For the thinnest sublayer n+1, its hybrid matrix is obtained directly by thin-layer asymptotic approximation:  Starting with this matrix, we implement self-recursions as This recursive algorithm proceeds until i = 1 and the layer hybrid matrix is found as f  H H . Throughout the procedure, there is no need to solve any eigenproblem and the hybrid matrix can be computed stably and accurately even for very thick or very thin layer.

Discussion and Numerical Results
The previous sections have discussed some algorithms for scattering and hybrid matrix methods.For concise comparison, Table 1 lists each of the algorithms and its pertaining equations involved for each major step represented by an arrow.In each major step, there is at least one (dense) matrix inversion to be dealt with, which often constitutes the most time-consuming operation.For the scattering matrix method, we list the algorithms in physics-based form as well as coherent form.The latter form helps to bring out the close resemblance with the algorithm of hybrid matrix method.Also listed in Table 1 is the input required for each algorithm.Since the scattering matrix relates wave amplitude vectors across interfaces, the input ought to be eigenwaves of two layers.As for the hybrid matrix that relates field variables, the input may need only the eigenwaves of individual layer.Through the recursive asymptotic method, the input does not invoke any eigenwaves at all.
To highlight the distinctions between the hybrid ma- With eigenwaves (5) ( 7) Hybrid matrix method: (26) (30 Recursive asymptotic hybrid matrix method: (34) (38) (30 trix method with eigenwaves and the recursive asymptotic method without eigenwaves, we further list down below the key steps in their respective procedure.In particular, the procedure with eigenwaves is i) Solve the eigenvalue problem (5) for wavenumbers and eigenvectors ii) Perform upward/downward-bounded eigenvectors sorting or selection in ( 6)-( 7), noting the boundedness/ radiation condition and degeneracy treatment if needed iii) Derive the layer hybrid matrix using (26).
Step i) is often time-consuming, while step ii) deserves much careful attention and could be rather bothersome in practice.On the other hand, the procedure without eigenwaves via the recursive asymptotic method is i') Initialize the thin-layer asymptotic approximation (34) directly from f A ii') Perform self-recursions (35)-(38) until i = 1 iii') The layer hybrid matrix is found as (1) f  H H .All steps here are straightforward and involve elementary matrix operations only.
To assess the accuracy of recursive asymptotic hybrid matrix method, we investigate the relative error changes with the number of geometric subdivisions n+1 or equivalently, the recursion number n.We arbitrarily take a ZnO layer of 1 μm thick at 1 GHz as an example.H represent the layer hybrid matrix obtained from the recursive asymptotic method and eigenwaves solution, respectively.The error is calculated by taking the average over a range of transverse wavenumbers.Notice that the error decreases initially due to smaller truncation error for smaller initial sublayer thickness n d .After certain minimum point, the error increases slightly and reaches a plateau without increasing further.To illustrate the usefulness of recursive asymptotic hybrid matrix method, let us consider a ZnO/ diamond/Si structure at 2 GHz.The thicknesses of ZnO and diamond layers are 1.2 and 10 μm respectively, while Si substrate and vacuum are assumed semi-infinite.For analysis of surface acoustic wave (SAW) on such structure, one can derive the generalized Green's function matrix G defined by where s  is the charge density on the surface.G can be formulated using the scattering matrix with eigenwaves in a robust manner, see [5].Alternatively, one can also determine G using the stack hybrid matrix [1: ] N H (for a stack from layer 1 to N) as where 0  is the permittivity for vacuum (layer N+1), ( ) and sub Z is the characteristic surface impedance for Si substrate (layer 0) The stack hybrid matrix [1: ] N H can be obtained with or without eigenwaves solution as mentioned earlier.Figure 4 shows the Green's function element computed with and without eigenwaves.In the latter case, we apply the recursive asymptotic hybrid matrix method with n=6.Although this recursion number is rather small, the results agree quite well and the plots are barely distinguishable.Referring to Figure 3, one can select higher recursion number for better accu- racy, although this may not be needed in many cases (e.g. when material data is not that accurate).In general, the computation efficiency is improved for lower accuracy required and also for thinner layer with fewer geometric subdivisions.Therefore the method provides a very convenient way that facilitates the trade-off between computation efficiency and accuracy.Note that the efficiency improvement here is meant for every layer and one will gain substantial savings in the total computation time when there are many layers in the stack for modeling inhomogeneous media.Moreover, the method is very useful for being simple enough since it does not require any eigenwaves for all layers (even semi-infinite substrate).Thus, it may be applicable even when the eigensolver package is not readily accessible, such as on light-weight multi-thread processors (e.g.GPUs).

Conclusions
This paper has presented the recursive asymptotic hybrid matrix method for acoustic waves in multilayered piezoelectric media.The hybrid matrix method preserves the numerical stability and accuracy across large and small thicknesses.For discussion and comparison, the scattering matrix method has also been presented in physics-based form and coherent form.The latter form resembles closely that of hybrid matrix method and helps to highlight their relationship and distinction.For both scattering and hybrid matrix methods, their formulations in terms of eigenwaves solution have been provided concisely.Making use of the hybrid matrix, the recursive asymptotic method without eigenwaves solution has been described and discussed.The method bypasses the intricacies of eigenvalue-eigenvector approach and requires only elementary matrix operations along with thin-layer asymptotic approximation.It can be used to determine Green's function matrix readily and facilitates the trade-off between computation efficiency and accuracy.

ff
formed by generalized stress vector f σ (comprising normal stress f τ and normal electric displacement zf D ) and generalized velocity vector f υ (comprising velocity f v and the rate of change of electric potential

Fig- ure 3
shows the average relative error versus recursion number n.The relative error is measured by -

Figure 3 .
Figure 3. Average relative error vs. recursion number n.

Figure 4 .
Figure 4. Green function element computed with and without eigenwaves (via recursive asymptotic hybrid matrix method with n = 6).