Spectral Element Method for the Elastic/Acoustic Waveguide Problem in Anisotropic Metamaterials

In order to simulate elastic wave propagation in a complex structure with inhomogeneous media, we often need to obtain the propagating eigenmodes of an elastic waveguide. As the waveguide is assumed uniform in one direction, the original 3-D problem can be converted into a so-called 2.5-D problem by using the Fourier transform in that direction. However, the introduction of elastic metamaterials (EMM) broadens the horizon of this subject, and new features are required in EMM waveguides that cannot be obtained by most traditional waveguide solvers. In this work, a spectral element method (SEM) is developed to simulate the elastic/acoustic waveguide problem in anisotropic media with anisotropic mass density and/or negative index parameters. To the best of our knowledge, the SEM has not been introduced previously for such a waveguide problem. For waveguides with anisotropic density that cannot be solved by the FEM in most of commercial software packages, we design an anisotropic density EMM waveguide with our SEM solver to demonstrate some intriguing phenomena. The spectral element results are verified by several numerical examples through comparison with the traditional finite element method (FEM) to show its significant advantages in term of accuracy and computation efficiency.

In order to simulate elastic wave propagation in a complex structure with inhomogeneous media, we often need to obtain the propagating eigenmodes of an elastic waveguide. As the waveguide is assumed uniform in one direction, the original 3-D problem can be converted into a so-called 2.5-D problem by using the Fourier transform in that direction. However, the introduction of elastic metamaterials (EMM) broadens the horizon of this subject, and new features are required in EMM waveguides that cannot be obtained by most traditional waveguide solvers. In this work, a spectral element method (SEM) is developed to simulate the elastic/acoustic waveguide problem in anisotropic media with anisotropic mass density and/or negative index parameters. To the best of our knowledge, the SEM has not been introduced previously for such a waveguide problem. For waveguides with anisotropic density that cannot be solved by the FEM in most of commercial software packages, we design an anisotropic density EMM waveguide with our SEM solver to demonstrate some intriguing phenomena. The spectral element results are verified by several numerical examples through comparison with the traditional finite element method (FEM) to show

Introduction
Recently, elastic waveguide problems have gained much attention due to various engineering applications, for example, ultrasound characterization, non-destructive testing, and structural health monitoring 1,2,3 . For various types of elastic waveguides 4 , the mode analysis is an important research topic, because wave propagation and scattering phenomena in a waveguide can be described as the superposition of all of the propagation modes and evanescent modes. In this class of problems, one is interested in solving the propagation constants and the corresponding field distributions of individual modes in a given waveguide structure.
As an extension of the electromagnetic and acoustic waveguides 5 Numerical examples show the significant advantages of the SEM in term of accuracy and computation efficiency compared with the conventional FEM as implemented by COMSOL. We also design the simulation of waveguides with anisotropic density that cannot be solved by commercial solvers.
The organization of this paper is as follows. In Section II, we will present the detailed weak formulation of the elastic waveguide. The discretization by the SEM is shown in Section III. Finally, the accuracy and efficiency of the SEM are demonstrated by several examples in Section IV.

Governing Equations
For a general anisotropic and inhomogeneous elastic metamaterial with a potentially anisotropic mass density, elastic wave equations in frequency domain read = (∇u + ∇u T )/2 (1c) where ρ = (ρ ij ) 3×3 is the anisotropic mass density; ω denotes the angular frequency, u is the particle displacement; , τ are the 2nd-order strain and stress tensors; c is the 4th-order elastic tensor. For Voigt notation, c ijkl can be converted to second-order tensor (C rs ) 6×6 . The subscripts of C and c satisfy the relations between (r, s) and (i, j, k, l): 1 ↔ 11, 2 ↔ 22, 3 ↔ 33, 4 ↔ 23, 5 ↔ 13 and 6 ↔ 12. Therefore, the constitutive equation (1b) can be transformed into a matrix form 32 . And, the divergence of τ computed by the left divergence operator can be expressed as where Einstein's convention has been adapted, with the repeated indices implying summation. The operators "⊗" and "·" represent the diadic product the dot product, respectively. Substituting (2) into (1a) yields The Latin subscripts i, j, k, . . . represent three-dimensional indices and the Greek subscripts α, β, . . . are two-dimensional indices. It is well known that the waveguide problem is actually a 2.5-dimensional problem, where the field is three-dimensional depending on (x 1 , x 2 , x 3 ) = (x, y, z) but the materials are two-dimensional depending on (x 1 , x 2 ). When the propagation is along the +z-axis and the cross section of the waveguide is uniform in the z-direction, the phasor expression for displacement field u and the operator ∇ can be written explicitly as for any given waveguide mode, whereê i and u i (x, y) denote the i-th unit vector and its corresponding component of the displacement field in Cartesian coordinates, respectively, and γ z = jk z = α z + jβ z is the complex propagation constant (the real variables α z and β z are called the attenuation constant and phase constant respectively). k z is the z-component of the wave vector. In the following formulations, the time convention e jωt is omitted. Therefore, substituting (4) into (1c), the strain tensor can be written as On the other hand, the stress tensor can be expressed as τ = e −γzz τ ijêi ⊗ê j , and the components τ ij are indicated as where the subscripts of C come from the elements of a symmetrical constant matrix r defined by Obviously, τ is also symmetrical. Inserting (4)-(6) into (3), we can obtain the governing equation of the elastic waveguide as follows • Tensor formulation: • Component formulation: aê i ⊗ê j ·ê k =ê i δ jk ,ê i ·ê j ⊗ê k = δ ijêk andê i ⊗êα :ê β ⊗ê j = δ ij δ αβ b For a given k, for example k = 3, 3ijmê3 ⊗ê i ⊗ê j ⊗êm is a 3rd-order tensor. Thus, we write it as 3ijmêi ⊗ê j ⊗êm for short.
Spectral Element Method for the Elastic/Acoustic Waveguide Problem in Anisotropic Metamaterials 5 • The components of the coefficient tensor = kijmêk ⊗ê i ⊗ê j ⊗ê m are denoted by It is easy to see that (8) is a quadratic eigenvalue problem, where γ z = jk z = α z + jβ z is the eigenvalue and u denotes the corresponding eigenvector. The goal of this work is to develop the SEM for solving the eigenpairs (γ z , u).

Weak Formulation
Based on the framework of FEM, multiplying equation (7) by the test function ϕ and integrating and using the integration by parts for the second and fourth integrations, we arrive at the weak form equation for the solid region The above bilinear functions can be written in detail as where Γ is the cross section of waveguide; ∂Γ denotes the boundary of Γ;n = n αêα represents the unit outward normal vector at the point on the edge ∂Γ; the superscript " †" means the complex conjugate. Because of the existence of boundary integral items I 1 and I 2 , we need one or the combination of the following boundary conditions.

Boundary Conditions
In order to solve the propagation constant γ z within a given waveguide, for the external boundary integral −γ z I 1 + I 2 , we need the suitable boundary conditions, such as the hard BC, the soft BC, the BPBC and the ABC.
(2) Soft BC readsn · τ = 0. According to (6), we have Substituting (12) into the external boundary integral, we have −γ z I 1 + I 2 = ∂Γn · τ dx in view ofn · τ = n j τ ji and n 3 = 0 in the waveguide problem. Therefore, the boundary integrals vanish and the region of the integration is unchanged.
(a) By the Bloch theorem 33 , we first obtain where k = k t +ẑk z is the Bloch wave vector, r and a are the position vectors on the boundary ∂Γ and the lattice translation vector, respectively. Define the Bloch periodic subspace where For any ϕ, u j belonging to H B p , from (11d) and (11e), note that the normal vectors defined on a pair of periodic boundary have opposite directions, we can obtain the external boundary integral I 1 = I 2 = 0. Consequently, the external boundary integrations vanish for all the opposite boundaries by using the BPBC on the external boundary. (b) Meanwhile, the BPBC waveguide problem can be transformed into the equivalent periodic boundary conditions (PBC) waveguide problem. When u j are written as the plane wave form u j (k t , r) = u j (k t , r)e −jkt·r , we can obtain the periodic boundary conditions u j (k t , r + a) = u j (k t , r) arising from (13). Therefore, the corresponding periodic subspace can be defined by For any ϕ, u j belonging to H p , it is easy to check that the external boundary integrals are still zero. By replacing the operator ∇ t with ∇ t − jk t in (10), we arrive at a new scheme where the bilinear functions can be written as In our work, the ABC is used to truncate the infinite external boundary when the waveguide is unbounded in the transverse directions. From 25 , the ABC for the scalar mass density is expressed as where t is the boundary traction, v is the velocity field on the surface, and C L , C T represent the material bulk speed of longitudinal and transverse waves in the background outside the computational region Γ, respectively. However, in this work we treat anisotropic density, so the boundary function is rewritten aŝ Accordingly, the scalar expressions with the Einstein notation are expressed as Moreover, noting that n 3 = 0, the external boundary integrations in (10) are replaced by The above bilinear function can be written in detail as The weak formulations of the elastic waveguide for the pure solid model with the BPBC and the ABC are expressed compactly as The above is for the case where all materials are solid in the waveguide. If part of the waveguide is made of fluid, we need to consider the special coupling between fluid and solid in the eigenvalue problem.

The Fluid-Solid Coupling System
When the waveguide is filled with an inhomogeneous medium including parts of fluid and solid, we need to consider the fluid-solid coupling. The symbol "f " is introduced to denote the local fluid region f , which shares the common interface ∂Γ f s with the local solid region s shown in Fig 1. For the fluid-solid coupling system, we not only derive the weak form in the fluid region, but also give the continuity condition of the fluid-solid at the interface defined by ∂Γ f s . First, the governing equation for the potential χ, defined as in the fluid region is introduced from 26 where ρ f is the density of fluid/acoustic materials and κ is the bulk modulus. The phasor expression for potential is shown as χ = χ(x, y)e −γzz . Multiplying (24) by the test function ψ and integrating, after using the integration by parts, we obtain the weak form for the fluid region The above bilinear functions can be expressed as follows Second, when the boundary integral −γ z I 1 + I 2 in (10) is restricted to the interface ∂Γ f s between the solid region and the fluid region, because of the continuity condition of the tractionn · τ = jωχn 26 and the fact n 3 = 0, it follows that Finally, for the fluid region, similarly, the external boundary integration I 4 satisfies On the other hand, by replacing the normal component of the velocityn·v fluid =n·(ρ −1 f ∇χ) in the fluid region with the normal component of the velocityn · (jωu) in the solid region 26 , I f s 4 can be derived as Spectral Element Method for the Elastic/Acoustic Waveguide Problem in Anisotropic Metamaterials 9 Compactly, the weak formulations of the fluid-solid coupling system are shown as Note that, when the cladding medium outside a core of the waveguide is filled with solid, −γ z I 1 + I 2 is replaced by I 3 shown in (23b) and I 4 = 0. Conversely, when the cladding is a fluid region, −γ z I 1 + I 2 = 0 and I 4 is shown in (27). The above completes the formulation of elastic waveguide and its weak formulations. In the next section, we will introduce the discretization scheme to calculate the propagation constants γ z of the waveguide and their corresponding modes (eigenvectors).

Basis Functions
In order to approximate the unknown field component u j , we apply the GLL polynomials as the basis functions. The N th-order 1D GLL polynomials are defined as where the interpolating points ξ r ∈ [−1, 1], and they are chosen as the GLL points which are the roots of equation (1−ξ 2 r )L N (ξ r ) = 0, and L N (ξ) is the derivative of the N th-order Legendre polynomial. Note that, the Legendre polynomials are orthogonal polynomials that allow to reduce the interpolation errors compared to the standard Lagrange polynomials used in the FEM. u j can be approximated by using the tensor-product ϕ s (η) of two 1D nodal basis functions, where the subscript p is the compound index of (r, s). Let the physical domain be subdivided into a number of non-overlapping quadrilateral elements, so that each element can be mapped into the reference element [−1, 1] × [−1, 1] by the mapping x(ξ, η), y(ξ, η) 5,34 . For example, the irregular element κ with curved edges can be mapped to the reference elementκ by using the curvilinear mapping shown in the corresponding invertible mappings are applied to the basis function ϕ(x, y) = ϕ(ξ, η) is the Jacobian matrix, as derived in 24,30,31 .

Discrete Forms
In general, three unknown components u j of the displacement field can be approximated by where n sj represents the number of nodal degrees of freedom (DOF) of the component u j for the solid region. Thus, the total number of DOF in the solid region is Inserting (31) into (18), we arrive at the quadratic eigenvalue problems 1 , · · · , u j,n sj ), the subscript "s" means the solid region andT is the boundary integral matrix, which is equal to zero when using the Hard BC, the Soft BC and the BPBC, and nonzero for the ABC. After the invertible mapping, the elemental matrices consist of the following parts where p, q = 1, 2 · · · , N + 1, the superscript "(κ)" means the reference element and J b arises from the mapping from any edges to reference domain [−1, 1]. Meanwhile,T s arising from the ABC is expressed in (33). Similarly, u j , χ are written as Spectral Element Method for the Elastic/Acoustic Waveguide Problem in Anisotropic Metamaterials 11 where n f denotes the number of the total nodal DOF for the fluid region. Substituting (34) into (29), we arrive at the fluid-solid coupling eigenvalue problem based on the BPBC ], w (w 1 , · · · , w n f ) . The elemental matrices are given as following For the remaining elemental matrices, they can be obtained by replacing the superscript (κ) of (33) with (s). After these matrices are assembled, the quadratic eigenvalue problems (35) is converted to a first order generalized eigenvalue problem in (37) for γ z referring to 35 , where I and 0 denote the identity matrix and zero matrix, and then it can be solved by using the eigenvalue solver "eigs" in MATLAB based on ARPACK library routines.

Numerical Results
In this section, several examples are presented to verify the high accuracy and efficiency of the SEM for simulating the elastic waveguide problems. The memory, the number of degrees of freedom (DOF) and the accuracy for our method are compared with the commercial FEM solver COMSOL. Finally, we conduct a numerical experiment on an elastic matematerial (EMM) core which cannot be solved by COMSOL, because of the presence of anisotropic density. Before the experiments, there are some preparations. First, in the simulation of the BPBC waveguide problem, the wave vector is defined by k = k(x sin θ cos φ +ŷ sin θ sin φ +ẑ cos θ) where k = ω/v i , v i is the velocity of the P wave or S wave in the background medium and (θ, φ) are the elevation and azimuthal angles of the propagation direction. Second, for convenience, we introduce the notations in our tables and figures: 1) v p is the velocity of the P wave (longitudinal wave). 2) v s is the velocity of the S wave (transversal waves).
3) ρ is the mass density. 4) λ and µ are Lamè constants. 5) k N i,z is the i-th eigenmode wavenumber k z obtained in the z direction by the N -th order SEM. 6) The reference valuē k 10 i,z is the solution of the 10th-order SEM with an extremely fine mesh. The relative error is calculated by |k N i,z −k 10 i,z |/|k 10 i,z |. 7) The computational time and memory are displayed with the "tic", "toc" function and "memory" function in Matlab, respectively. Last but not least, for the quadratic eigenvalue problems, the solver will provide two opposite eignvalues (γ z and −γ z ). To determine the correct sign of the propagation constant, we introduce two quantities: the time averaged Poynting vector 36 p = Re(−jωu·τ † )/2 and the corresponding power P z = Γẑ · pdxdy in the cross section. The positive P z is the criterion for choosing the correct sign of γ z in the following numerical examples. The SEM method is implemented by using Matlab on a MacBook Pro 2018 PC with 16 GB Memory and Intel Core i7 CPU. COMSOL was used for comparison verification on the same PC.

Bloch periodic unit cell
In order to verify the accuracy and convergence of the proposed SEM, we first consider a simple inhomogeneous anisotropic waveguide with the BPBC. The Bloch periodic unit cell has many applications in lithography and the design of elastic metasurfaces, which act as a plate-like waveguides connecting two elastic half-spaces 37 . The configuration of the unit cell is shown in Fig. 3, where nine circular lead cores are embedded in the zinc square lattice. These circles with different radius are spaced one millimeter apart. The material properties are λ Pb = 3.142 × 10 10 N/m 2 , µ Pb = 5.986 × 10 9 N/m 2 , ρ Pb = 11340 kg/m 3 . The cladding is a transversely isotropic material with c 11 = 16.5 GPa, c 12 = 3.1 GPa, c 13 = 5.0 GPa, c 33 = 6.2 GPa, c 55 = 3.96 GPa and ρ Zn = 2700kg/m 3 . The frequency f = 5 MHz and the unit cell is 2 cm containing multiple wavelengths, so that it is a large scale problem.
The numerical results of k N i,z obtained by the 5th-order SEM, the 5th-order FEM in COMSOL and the 10th-order SEM with an extremely fine mesh are shown in Table 1 (the negligible imaginary part is not shown). They are denoted by SEM-K5, COMSOL-K5 and SEM-K10, respectively, and in view of the maximum interpolation order of COMSOL is only 5, SEM-K10 is taken as the reference value. It is observed that SEM-K5 matches excellently with both COMSOL-K5 and SEM-K10. On the other hand, as illustrated in Table 2, to achieve similar accuracy, COMSOL requires more 2.02 times DOFs, 2.35 times CPU time and 1.54 times memory than the SEM. We can also see that when the numbers of element and DOF are taken to be similar, COMSOL is not as accurate as the SEM and requires a little more computational costs. Thus, the proposed SEM is more efficient than the FEM method, mainly because of the spectral accuracy (the exponential convergence) of the SEM shown in Fig. 5. The magnitude distributions of u(x, y) for the 1st, 7th, 18th mode are displayed in Fig. 4. All of them propagate in the lead core with different radii.

Resonant structure of an EMM
In order to verify that our SEM solver is accurate and efficient for the inhomogeneous solidfluid coupling BPBC waveguide, we first consider a resonant structure in the building block of a left-handed material (LHM) proposed in 18 . This kind of resonant structure will bring negative elastic parameters within a certain frequency range and the cross section of the unit cell, a rubber coated water cylinder embedded in a foam host, is shown in Fig. 6. The lattice constant is a and the radius of the rubber and water is 0.32a and 0.24a, respectively.  When we set a = 1 m, the corresponding frequency is chosen as 34.887 Hz referring to 18 . In addition, the material parameters are listed in Table 3 and the BPBC is used in the example, (θ, φ) = (0, 0) and k = ω/v s , v s is the velocity of S wave in the foam. Table 4 shows that the numerical solutions of the inhomogeneous isotropic BPBC waveguide obtained by the SEM and COMSOL agrees well. On the other hand, as illustrated in Table 5, the proposed SEM is more efficient than the FEM in terms of the DOF and memory. Moreover, from the subgraph (a) and (c) of Fig. 7, we can see that there is a quadrupolar resonance in the rubber region for the first mode produced by the P wave, due to the much smaller v p of the rubber than those in the background foam and the water core. Besides, as shown in the subgraph (b) and (d) of Fig. 7, a total reflection occurs at the boundary between the rubber region and the water core for the second mode, because of the much larger v p of the water in Table 3.

Optical fiber model
Next, to verify the ABC formulation for an open (unbounded) inhomogeneous isotropic waveguide, we consider the optical fiber. It is a common optical waveguide consisting of the cladding and the fiber core, and its elastic waveguide properties are of significant interest Spectral Element Method for the Elastic/Acoustic Waveguide Problem in Anisotropic Metamaterials 15 foam rubber water Fig. 6. The cross section of the resonant structure for a left-hand material with a rubber coated water cylinder embedded in a foam host, with their material properties listed in Table 3. Table 3. Parameters for the Resonant Structures in Figure 6.   38,39 . The cross section of the optical fiber is shown in Fig. 8, which consists of the core and the cladding. The radius of the core and the cladding is a = 4.1 µm and 3a = 12.3 µm, respectively; the cladding is pure SiO 2 and the core is filled with one of the three different materials as shown in Table 6. The SEM is employed to simulate the elastic waveguide properties of this optical fiber. Besides, to verify the accuracy and effectiveness of the SEM for solving the solid-fluid system, a fluid cladding is also considered. The material parameters are included in Table 6. In order to simulate the unbounded waveguide structure, the ABC is used to truncate the cladding so that the simulated structure mimics an infinite cladding region.
(1) Normal Elastic Materials First, to verify the accuracy of the SEM solver for the inhomogeneous open waveguide problems, we conduct a numerical experiment on the actual optical quartz fiber model consisting of cladding 1 and core 1, and the frequency is chosen as 3 GHz as in realistic application 40 . The agreement among the three results in Table 7 verifies the accuracy of our scheme. Besides, we observe that the real part of the higher-order mode gradually decreases while the imaginary part falling into different orders of magnitude gradually Table 5. The comparison of FEM and SEM for the elastic resonant structure waveguide in Figure 6. core cladding a 2a x y Fig. 8. The cross section of the optical fiber (a=4.1 µm), with an unbounded cladding truncated by an ABC at r = 3a. The core and cladding materials can take the combination of materials listed in Table 6.
increases. The phenomenon indicates that the energy loss of the higher-order mode gradually increases. The relative errors and computational costs of the SEM and the FEM are illustrated in Table 8. For a similar mesh, the proposed 3rd-order SEM and 4th-order FEM can achieve similar accuracy (3E-6). The corresponding memory used by FEM is more than SEM, illustrating the proposed SEM is more efficient than FEM. Moreover, it can be observed that the 6th-order SEM can achieve higher accuracy (3E-7) with less memory, due to the spectral accuracy of the SEM solver. Incidentally, no impurity is present in our contour maps in Fig. 9, which indirectly indicates that no spurious modes exist in our method as discussed in 40 . In addition, there are convincing Table 6. Parameters for the Cores and Claddings of the Optical Fiber in Figure 8. explanations for the spurious modes. In general, spurious modes are obtained in the following two cases. One is that basis functions cannot describe the physical properties of solutions. In this manuscript, the GLL polynomials are employed to construct the basis functions which are obviously continuous at the interpolation points. The other one is that a discrete space cannot compactly approximate the solution space H 1 (Γ). For our method, the discrete space Q N,h = span{φ 1 , φ 2 , ·, φ N sj } is used to approximate the solution space H 1 (Γ), so that it is compact. In conclusion, there are no spurious modes in our method. Moreover, waves are well absorbed at the outer absorbing boundary. Table 7. kz (Mrad/s) of the elastic fiber-optics waveguide in Figure 8 obtained by the SEM and COMSOL for Core 1 and Cladding 1 listed in Table 6.  Table 8. The comparison of FEM and SEM for the elastic fiberoptics waveguide in Figure 8 with Core 1 and Cladding 1 listed in Table 6. (2) Double Negative Index Elastic Metamatrial (EMM) Core Second, for the same size model, we now consider the effects of the EMM core with a negative index. We design an example on a simultaneously negative mass density and bulk modulus EMM core 2 constructed by reference 15 , embedded in Cladding 1 in Table 6. The frequency is chosen as 0.3 GHz. Through calculating the velocities of P-wave and S-wave respectively shown in Table 6, we find the velocity of S-wave is an imaginary number, thus the S-wave is forbidden in this material. Again, k N i,z (N =5,10) obtained by the two methods are shown in Table 9. It is observed that the SEM solution matches excellently with the reference results and the COMSOL's results, verifying that our scheme is suitable for the negative index materials. Meanwhile, as illustrated in Table 10, DOF and the memory used by COMSOL (N =5) is 2 and 1.5 times more than SEM (N =5) to achieve the similar accuracy (8E-7). Evidently, it shows the high computation efficiency of the SEM. Furthermore, the magnitude distributions of u(x, y) corresponding to k i,z (i = 1, · · · , 6) with EMM core 2 are plotted in Fig. 10. On the other hand, instead of core 2, we conduct another experiment on core 3, whose density and bulk modulus are positive. The agreement in Table 10 verifies the accuracy of the results. Fig. 11 plots the distribution of u corresponding to k i,z (i = 1, 2, 3) with core 3. In contrast to the previous configuration, we observe in Fig. 10 that these modes in the waveguide of EMM core 2 propagate only at the interface between the core and cladding because of the presence of the negative index material.
Furthermore, we notice that one propagation mode exists under low frequencies (f d/c ∈ [6.5 × 10 −6 , 9 × 10 −3 ], d = 2a) as shown in Fig. 12. Within this frequency range, different from normal elastic materials, the increase of the frequency does not alter the distribution interval of the real part of k z (the phase constant β z = −5.23 × 10 5 ) with the negligible imaginary part (the attenuation constant α z ) on the basis of the positive Table 9. kz (Mrad/s) of the elastic waveguide with either EMM Core 2 or Core 3 and Cladding 1 in Figure 8 obtained by the SEM and COMSOL. EMM core 2 core 3  Fig. 10. Magnitude distributions of u(x, y) for eignmodes corresponding to k i,z (i = 1, · · · , 6) obtained in the open fiber-optics waveguide problem with EMM core 2 in Cladding 1 in Figure 8. (a)-(f) correspond to the first to the sixth mode. All of them propagate only at the interface between the core and cladding because of the presence of the negative index material, different from the normal material Core 3 in Fig. 11. P z . Through the observation in Fig. 13, we can find the propagation mode is caused by the P wave and concentrated in the core. In addition, in order to explain the existence of this mode, the phase velocity v p = ω/β z is shown in Fig. 12 (b). It can be found that this mode exhibits backward wave propagation in the cross section, which is defined as the phase velocity direction (−ẑ) antiparallel to the Poynting vector (+ẑ), caused by the negative-index materials 41 . Hence, different from normal elastic materials, the application of EMMs will bring some special eigenmodes in the elastic waveguide. The phase velocity. The negative value means the direction of the phase velocity is −ẑ, antiparallel to the +ẑ. Thus, the backward wave propagation phenomenon is found in this mode.

(3) Solid-Fluid Coupling Model
In the previous case, the cladding was assumed unbounded, which may not be realistic. Actually, the external medium of the most practical open waveguide problems is fluid (for example, either air or water). Therefore, here we examine the same size model to verify the fluid-solid coupling system of the open waveguide problems at the frequency of 60 MHz. In this case, the cross section consists of the solid core 3 and the fluid cladding 2. Noting that COMSOL does not provide the ABC in the modal analysis of the acoustic module. So for comparison, we set the impedance value of the plane wave as an approximation in COMSOL when the outer boundary is far enough. The agreement is good as illustrated in Table 11, demonstrating that the proposed SEM is capable of treating the ABC solid-fluid problem. Besides, the relative errors obtained by different orders of SEM confirm the exponential convergence in Fig. 14. Next, we give a detail discussion about the third mode, whose attenuation constant is almost zero. First, as observed in Fig. 15, different from other modes, the propagation of this mode concentrates in the core. The reason for this phenomenon is that the third mode may be caused by the transversal wave, which cannot be transmitted into the fluid region. Moreover, same as the investigation in 30 , for an exact integration ( (N +1)th-order GLL quadrature in each element) of the second-order geometrical modeling, the errors of mode 3 are straight lines if one groups the even and odd orders separately, and the even and odd orders have different offsets. Therefore the relative error of this mode is reasonable.

EMM with Anisotropic Density
In addition to the metamaterials with negative index discussed above, the metamaterials with anisotropic mass density have attracted more and more attention recently. Because the equivalent model with effective anisotropic mass density can describe the dynamic behavior of the original lattice system in all directions. Hence, we conduct one numerical experiment on one anisotropic density core that cannot be simulated by some traditional numerical methods. Besides, in this section, the formulation of phase velocity obtained by elastic tensor C and isotropic density ρ in literature 42 is extended to one suitable for anisotropic densityρ = (ρ ij ) 3×3 . Note that the explanations of the symbols are referred to reference 42 . Starting with time domain governing equation Multiplying both sides by the inverse ofρ, the scalar expressions for a homogeneous medium are arrived at After denoting Γ ik = c ijkl n j n l , where n j is the component of the unit propagation vector n and multiplying both sides by p l , the components of the unit polarization vector colinear with the displacement. The final eigenvalue formulation is obtained where v 2 is the eigenvalue. The cross section centered at (0,0) m is shown in Fig. 16, the width of the square cladding and the square core is 0.5 m and 0.11 m, respectively. The cladding is zinc with isotropic material parameters {v p , v s } = {4820.7, 2361.6} m/s and the mass density is 7100 kg/m 3 . The core is an anisotropic EMM with the effective elastic coefficients C 11 = 36.63 GPa, C 12 = 5.57 GPa, C 13 = 13.53 GPa, C 22 = 18.83 GPa, C 23 = 7.84 GPa, C 33 = 48.38 GPa, C 44 = 12.41 GPa, C 55 = 6.69 GPa, C 66 = 2.272 GPa. The frequency we choose is f = 16 kHz and corresponding effective anisotropic mass density represent ρ EMM = diag{6277, 3168, 2700} kg/m 3 according to 17 . The anisotropic mass density is frequency-dependent and caused by the different locally resonant frequencies along different directions in the microstructure design, depending on the inverse proportional function ρ eff,i = a + b/(ω 2 i − ω 2 ) 17,43 , where a, b are the positive constants given by the detailed model parameters and ω i is the locally resonance frequency along the i direction (i = x, y, z). Note that ω 1 is the smallest, the ρ EMM is certainly produced by the frequency below and close to the ω 1 , leading to the resonance phenomena dominated by the u x . Besides, the velocity of EMM alongn = (0, 0, 1) calculated through equation (38) is {v p , v s1 , v s2 } = {4233, 1979.2, 1032.4} m/s, smaller than the cladding. Therefore, the ABC is used to truncate the cladding. First, the good agreement between the 5th-order SEM numerical results and the 10th-order results of the extremely fine mesh is demonstrated in Table 12. Then, the relative errors of the first three modes plotted in Fig. 17 indicate the exponential convergence. Moreover, the magnitude distributions of u are plotted in Fig. 18 and all of them in the xy plane are along the x direction and dominated by the x-component of the u. On the other hand, we conduct another experiment on an normal anisotropic elastic core with the isotropic mass density ρ = 3772.5 kg/m 3 (the geometric average 3 3 i=1 ρ ii ) and the same elastic coefficients C. The magnitude distributions of u in this waveguide with normal materials are plotted in Fig. 19 and they are dominated by the components along three principal axis respectively. In contrast to the configuration with an EMM core, it can be found that the propagation mode dominated by u y , u z shown in Fig.  19 (c),(d) cannot be obtained in the example with the anisotropic mass density core. The phenomena are due to the difference between the EMM core and the normal core in view of the locally resonance frequencies in each principal axis, which is caused by the different mass density tensors. In conclusion, the above explains the phenomena caused by the use of the EMM core with the anisotropic mass density and demonstrates the rationality of our results.

Conclusions
This paper presents a SEM solver for the general EMM waveguide problems with negative index and anisotropic mass density as well as normal materials. The solver can treat inhomogeneous and anisotropic solids, but also include the fluid-solid coupling. Meanwhile, Zinc EMM y 0 x Fig. 16. The cross section of the anisotropic density waveguide consist of a square EMM core and an unbounded Zinc cladding truncated by a square outer ABC boundary. Table 12. k n i,z of the anisotropic mass density EMM waveguide problem in Figure 16.   17. Relative errors of the first three modes for the anisotropic mass density EMM core waveguide problem in Figure 16.
the discussions about four boundary conditions (the hard BC, the soft BC, the BPBC, the ABC) are provided. Both excellent agreement between results and those from the commercial FEM solver COMSOL and less computational costs are demonstrated in the numerical validations. Moreover, some interesting phenomena brought by the application of the EMM can be observed in the numerical experiments, for instance, unusual modes generated by the negative refractive index or common modes eliminated by the anisotropic mass density.  Figure 16 replaced by an isotropic mass density ρ = 3772.5 kg/m 3 . (a)-(d) correspond to the first to the fourth modes. Note that the 2nd, 4th modes dominated by uy and uz respectively are absent in the EMM waveguide in Fig. 18.