Complex Rayleigh Waves in Nonhomogeneous Magneto-Electro-Elastic Half-Spaces

The complex Rayleigh waves play an important role in the energy conversion efficiency of magneto-electro-elastic devices, so it is necessary to explore the wave propagation characteristics for the better applications in engineering. This paper modifies the Laguerre orthogonal polynomial to investigate the complex Rayleigh waves propagating in nonhomogeneous magneto-electro-elastic half-spaces. The improved method simplifies the calculation process by incorporating boundary conditions into the constitutive relations, shortens the solving time by transforming the solution of wave equation to an eigenvalue problem, and obtains all wave modes, including real and imaginary and complex wavenumbers. The three-dimensional curves of full frequency spectrum and phase velocities are presented for the better description of the conversion from complex Rayleigh wave modes to real wave ones; besides, the displacement distributions, electric and magnetic potential curves are obtained in thickness and propagation directions, respectively. Numerical results are analyzed and discussed elaborately in three cases: variation of nonhomogeneous coefficients, absence of magnetism, and absence of electricity. The results can be used to optimize and fabricate the acoustic surface wave devices of the nonhomogeneous magneto-electro-elastic materials.


Introduction
The past few decades have seen considerable time and effort invested in studying the mechanics of magneto-electro-elastic composites. With their particular properties, such as the mutual conversion of magneto-electro-elastic energy [1], these materials are widely applied in sensors, ultrasonic nondestructive evaluation, and vibration control devices [2][3][4][5]. They are so promising for many applications that they have become the next generation multi-functional materials [6]. To improve their thermal and corrosive resistance, the research has gradually expanded to the field of nonhomogeneous material; thereafter, the study of the waves propagating in nonhomogeneous magneto-electro-elastic composites has received more and more attention [7,8], and some researchers have successfully fabricated an actuator with small mechanical damping and heat generation by using functionally graded piezoelectric material [9].
The study on the propagation behaviors in these materials is so important and urgent that many researchers have been working in this field. Pang et al. [10] investigated Rayleightype surface waves in a piezoelectric-piezomagnetic layered half-space. Wei et al. [11] explored the propagation characteristics of shear horizontal surface waves propagating in a layered piezoelectric-piezomagnetic half-space and discussed the sources that effect on the dispersion curves and phase velocities. For investigating the Love waves propagation in layered magneto-electro-elastic half-space, Bou Matar et al. [12] combined the Legendre and Laguerre polynomial technique to obtain their solutions. Du [13] not only investigated the Love waves propagating in layered magneto-electro-elastic half-space with initial stress conditions using the Heaviside step function; by solving the characteristic matrix, the eigenvalues give the solution of propagating, evanescent, and complex wavenumber waves in nonhomogeneous magneto-electro-elastic half-spaces; moreover, the eigenvectors are obtained to present the mechanical displacement, electric, and magnetic potentials. For examining the effect of different material constants on the complex Rayleigh waves, three cases are taken account of: nonhomogeneous magneto-electro-elastic half-space, nonhomogeneous piezoelectric half-space, and nonhomogeneous piezomagnetic half-space, and their results are compared and discussed in detail. This paper also investigated the effect of the nonhomogeneous variation on the complex Rayleigh wave propagation characteristics. Furthermore, the 3D dispersion curves are obtained for discussing the cut-off frequencies, and the phase velocity curves are given for exploring the wave conversion.

Basic Equations
Consider an orthotropic nonhomogeneous magneto-electro-elastic half-space, where the material coefficients vary continuously along the z direction. The Cartesian coordinate system (as shown in Figure 1) is designated as follows: the upper surface of the half-space is set as the horizontal coordinate plane x − y; the z axis is downward and perpendicular to the x − y plane; the half-space occupies the region [0, ∞) in the z direction. Rayleigh surface waves propagate in the x-axis direction. Let u i represent the components of particle displacement in the ith direction, and ρ be the mass density. When the electric and magnetic sources are absent, the governing equations in the half-space are in which .j denotes differentiation with x j , theü i represents the second derivative of u i over time, B is the magnetic induction, D is the electric displacement, and σ denotes the stress. The Einstein summation notation is applied here. For a linear, anisotropic and coupled mechanical-electric-magnetic material medium, the constitutive equations [38,39] are where E k , H k , and ε kl are the electric field, magnetic field, and strain, respectively; ik , d ik , and µ ik are, respectively, the dielectric, magnetoelectric, and permeability coefficients; and e kij , q kij , and c ijkl are the piezoelectric, piezomagnetic, and elastic coefficients, respectively. Based on the elasticity theory and the quasi-static Maxwell equation, the strain, electric, and magnetic field are related to the mechanical displacement u i , the electric potential φ, and magnetic potential ϕ, as follows: We just consider the waves propagating in the nonhomogeneous magneto-electroelastic half-space polarized along the z-axis direction. The upper surface are magnetically short and electrically open [24] in the present paper. Therefore, the mechanical and electricmagnetic boundary conditions at the upper surface are So, combining Equations (1) and (4) will obtained the solution of the Rayleigh waves propagating in the magneto-electro-elastic nonhomogeneous half-space.

Solution of the Problem
In this study, the variables of field are independent of y because the Rayleigh waves propagate along x-axis; therefore, Equation (1) can be written: By introducing the unit step function h(z): The above boundary conditions described in Equation (4) can be automatically incorporated into Equation (2), which leads to Assuming that the material properties are functions of depth, then the mechanical and electric-magnetic coefficients become the functions related to f (z) (it is assumed that they all are controlled by a same function for the convenience of expression). In practice, every material coefficient can be applied distinct continuous function, such as c(z) = c · f 1 (z), e(z) = e · f 2 (z) · · · ), and they can be defined as where e, q, , d, µ, and c are, respectively, the piezoelectric, piezomagnetic, dielectric, magnetoelectric permeability, and elastic constants in the homogenous material.
Assuming a Rayleigh surface wave propagating along x-axis, then the fields of mechanical displacement, electric, and magnetic potential can be written in the following form: Firstly, substituting Equations (7)-(9a-d) into Equation (5), and then reorganizing all these equations according to the orders of wavenumbers and categorizing them by the terms of k 2 , k 1 , and k 0 , therefore, Equation (5) can be expressed as: where the prime mark ( ) represents a derivative, and using it twice denotes the second derivative. When the Laguerre orthogonal polynomial series are introduced, U, W, X, and Y can be written as: in which, p i m and r i m (i = 1, 2) represent expansion coefficients and where L m (z) is the mth order Laguerre polynomial, and functions Q m (z) are orthonormal and complete from 0 to ∞, so appropriate for calculating in semi-infinite interval. Even though m can be any non-negative integer (0, 1, 2, 3, . . . ), practically, the summation of Equation (11) has halted at some value m = M because the solutions have converged well within a finite number of terms.
To begin with, substituting Equation (11) into Equation (10), and then using the complex conjugate Q * j (z) to multiply the obtained equations (j are integers from 0 to infinity), will lead to 4 × (M + 1) equations because any equation in Equation (10) will give (M + 1) equations. Thirdly, integrate the product over z in the range [0, ∞), and rearrange them according to k 2 , k 1 , and k 0 . Last, their matrix representation is achieved, as follows: (13) in which the vector and the square matrices can be expressed as the following Equation (14), and their elements can be obtained from Equation (10).
To find the solution of wave number k in Equation (13), here, we introduce a new column vector N: After inserting Equation (15) into Equation (13), the resulting equations can be rewritten as: Separating k and P in Equation (16), the result becomes: Considering N = kP in Equation (15) and introducing an identity matrix I, thus, Equation (17) becomes: where I has the same dimension as A. Letting finally, Equation (18) can be further simplified to: Equation (19) indicates that the wave number k is an eigenvalue, and R is an eigenvector corresponding to the k in this equation. So, giving a series of angular frequency values ω, a series of corresponding wave numbers k will be generated in the form of eigenvalues, which also can give the phase velocity with ease. The field distributions can be calculated by the eigenvectors R. As can be seen, this equation considers the angular frequency ω as independent variables, and k as the eigenvalues. The wavenumber of the equation exists in three forms, namely k 2 , k 1 , and k 0 . As a result, the equation can present all kinds of solutions, including the propagating and non-propagating wave numbers.

Numerical Results and Discussion
In this section, the propagation behaviors of complex Rayleigh waves in the nonhomogeneous magneto-electro-elastic half-space are studied numerically based on the formulations derived from the previous section. The material properties (at z = 0) can be seen in Table 1, which vary gradually in the thickness direction, and the nonhomogeneous function f (z) = (1 + z) n in Equation (8), where n is the nonhomogeneous coefficient. The mass density ρ is 7500 kg/m 3 at z = 0. The material constants in Table 1 are used for all examples, except for the section of approach validation, which used the material lead zirconate titanate (PZT-4).

Approach Validation
The verification for the present approach is carried out in two parts: comparing the mechanical displacement with earlier works and validating the boundary conditions. To the authors' knowledge, so far, the complex waves in nonhomogeneous magneto-electro-elastic half space have not been studied, so, here, we only present the displacement amplitude distribution in homogeneous PZT-4 half-space to make comparison with the available numerical results. The used material properties are taken from the literature [41]. The mechanical displacement amplitude distributions for the purely real solutions are given in Figure 2 for the order truncation M = 10. The curves in this figure represent normalized mechanical displacement (u/w 0 and w/w 0 ). The solid curves are obtained from the present paper when the wavenumber k = 2, and the dotted lines are taken from the literature [42]. Obviously, the excellent agreement between them solidly validates the accuracy of the present method.
Similarly, the normalized electrical displacement distributions are illustrated in Figure 3, in which its curves are normalized using the electric displacement on the upper surface. The figure gives the electric displacements of M at different values. As can be seen, although different M will generate different curves, their electric displacements all are zeros at z = 0. This is consistent with the description of Equation (4). All these show that the boundary conditions are met well. On the other hand, the curves converge with M increases, and the result when M is 9 is exactly the same as when M is 10, which shows that the solutions have converged well.

Full Frequency Spectra
By examining Equation (19), one can see that, if a series of angular frequency values are given, a series of corresponding wave numbers will be generated in the form of eigenvalues. As is well known, ω/2π will always give the real frequency f . The wave numbers obtained in the present paper can be categorized into three kinds of different roots: the purely real, the pure imaginary, and the complex. Taking the real wavenumbers as the x-axis, imaginary wavenumbers as the y-axis, and f as the z-axis, the full frequency spectra for a nonhomogeneous magneto-electro-elastic half-space will be constructed just like Figure 4. ξ = k × 10 2 in this figure is available everywhere in the present paper and will not be explained anymore. Figure 4a presents the full frequency spectra with the nonhomogeneous coefficient n = 1 in the function f (z) = (1 + z) n , the material constants used to calculate in the figure is listed in Table 1. Which shows that the wavenumber roots occur in pairs of opposite signs, i.e., if a + bi (a, b ∈ R) is a solution of Equation (19), then ±a ± bi all will be the solutions of it, and they are all distributed symmetrically in the top four quadrants as illustrated in the figure. For b/a = 0, the branches represent real waves (with blue dotted curves); For a/b = 0, the branches represent non-propagating waves (with red dotted curves in the figure), in which amplitudes decay exponentially; for ab = 0, the branches represent complex waves (with green dotted curves in the figure) traveling in the exponentially damped trigonometric wave forms with propagating distance. In 3D complex space, the purely real branches generally show the trajectories of first decline and then rise, while the purely imaginary and complex branches are very complicated. They do not show obvious regularity, and some of them are converted into propagating modes.  Figure 4. Full dispersion spectra of complex Rayleigh waves for a nonhomogeneous magneto-electroelastic half-space: (a) 3D full frequency spectra; (b) one quarter of the 3D full frequency spectra. Nonhomogeneous coefficient n = 1, ξ = k × 10 2 ; blue branches for real solutions, red for imaginary, green for complex.
Considering the symmetry of the wave's solution, here, we only present one quarter of its curves as shown in Figure 4b, then project all the curves onto the imaginary and real planes, respectively, as shown in Figure 5. It can be found that in this figure: For the purely imaginary branches, the first three branches I1, I2, and I3 are joined to purely real branches smoothly where the wavenumbers become zero; the purely imaginary branches become more complicated with the increase of frequency and the decrease of wavenumber. For the complex branches, some branches convert to purely real modes and continue propagating just as C1 branch; some branches convert to purely imaginary modes with the increase of frequency just as C2 branch. The real branches corresponding to propagating waves have been studied very deeply, so they do not necessitate the time to be discussed here. In a word, Figure 5 not only presents the values of wavenumbers and frequencies of the wave's solution in the nonhomogeneous magneto-electro-elastic half-space, but it also clearly illustrates the relationship between the three kinds of solutions.  6-8 present the 2D frequency spectra for the nonhomogeneous magnetoelectro-elastic half-space with different cases. Their material constants all come from Table 1, but there exist differences in the nonhomogeneous, electric and magnetic coefficients. The material constants in Figure 6 are the same as that in Figure 5, except for the nonhomogeneous coefficient n = 2; Figure 7 shows the frequency spectra without the magnetic field, and the material used in Figure 8 without the electric field. All waves propagate in the x-axis and their thickness directions are all along z-axis in these figures. The same as in the 3D figures, the point curves in blue correspond to real solutions, those in red correspond to imaginary solutions, and green to complex.   Comparing Figures 5-8 demonstrates that: nonhomogeneous coefficients have great influence not only on the cut-off frequencies of propagating waves, but, also on those of the complex and imaginary waves, the cut-off frequencies increase as the nonhomogeneous coefficients increase, e.g., the cut-off frequency of R1 branch is 24 when nonhomogeneous coefficient is 1, but it becomes 29 when the nonhomogeneous coefficient increases to 2; the cut-off frequency difference between the two branches R2 and R3 becomes larger with the increasing of nonhomogeneous coefficients, so that R3 crosses with the fourth branch. However, the electric or magnetic coefficients have little effect on the cut-off frequencies, which is verified in Figures 7 and 8. All these show that the influence of nonhomogeneous coefficients on cut-off frequencies is more significant than that of magnetic or electric coefficient.
Another finding is: there are more branches of green curves appeared in the right half field of Figure 7, while the quantity of this color curves has barely changed in the left half field compared with other three figures. This means that the ratio Im(ξ)/Re(ξ) of the complex solution becomes larger in this case. This indirectly shows that the magnetic coefficients have notable effect on the ratio. However, the ratio almost does not change in Figures 6 and 8. Therefore, the ratio Im(ξ)/Re(ξ) of the complex wave number can be changed by changing the volume fraction of the magnetic part in a composite material. The numbers of the imaginary branches of the complex solution (green curves in the left half field of figure) are basically the same in Figures 5-8, which shows that the magnetic, electric and nonhomogeneous coefficients have little effect on their imaginary part of the complex solutions in the magneto-electro-elastic nonhomogeneous half-space.

The Analysis of Displacement and Field Distributions
To explore the complex wave structures in the nonhomogeneous magneto-electroelastic half-spaces, this section will present their displacement, magnetic, and electric potential field distributions, at the same time compare them with those of the propagating waves. Here, we examine the points in the close vicinity f = 24 Hz, where is the transition area from complex waves to propagating waves, which can be seen in Figures 4 and 5. The exact solution can be obtained by the following steps: firstly, deriving the eigenvector R and wavenumber k from Equation (19), and immediately getting the numerical representation of vector P in Equation (14); and further getting the displacement, magnetic, and electric potential field distributions by Equations (9) and (11).
Next, two specific cases will be given to investigate the complex wave structures, their material coefficients are given in Table 1 and their nonhomogeneous coefficient n = 1. One is ξ = 2.2152, when f = 24.10 Hz, and another is ξ = 2.1365 + 0.1694i, when f = 24.07 Hz, which can be obtained through the method mentioned in the previous paragraph. Figures 9 and 10 show the displacement distributions along z and x in Figure 1, Their subfigure (a) are obtained at x is zero and (b) at z is zero, when f = 24.10 Hz, ξ = 2.2152 and f = 24.07 Hz, ξ = 2.1365 + 0.1694i, repectively. Figures 11 and 12 present the electric potential distributions, Figures 13 and 14 the magnetic potential distributions. The result, by comparing Figures 9a and 10a, this illustrates that there is barely change of the mechanical displacement field distributions in the thickness direction when the complex wave approaches to the propagating wave, and the similar phenomenon is also found in electric and magnetic potential distributions in the thickness direction. Figures 9a-14a demonstrate that the displacement, electric, and magnetic potential all decay to their final zero points at about z = 75, even though the attenuation trajectories of mechanical displacement are obviously different from that of the electric and magnetic potential. Figures 9b-14b show that: for complex waves, the electric and magnetic potential, same as in displacement, decay in the propagation direction x following a damped sinusoidal distribution rather than only exponential decay; the complex waves are able to travel farther and farther with their approaching purely real solution.

Effects of Nonhomogeneous and Magnetoelectric Properties on Phase Velocity
This section presents the phase velocities of complex waves for different cases, further discusses the effects of nonhomogeneous, magnetic, and electric coefficients on phase velocity. The phase velocity is defined as c p = ω/Re(ξ) in the present paper. Figure 15 presents the 3D phase velocity curves for a nonhomogeneous magneto-electro-elastic halfspace. In this figure, the real ξ is set to x-axis, the imaginary ξ is set to y-axis, and the phase velocity c p is set to z-axis; then, the purely real modes and the complex modes can be plotted in a three-dimensional coordinate system with ease, but the purely imaginary modes are not drawn here because their corresponding waves cannot propagate. The 3D phase velocity curves not only clearly present the purely and complex wave modes but also describe the process that the complex waves degenerate to purely real wave modes, e.g., the marked C point in Figure 15 illustrates where the complex wave mode converts to purely real mode. To express the phase velocity of converting point more clearly and facilitate the subsequent comparisons with different cases, here, we project all the curves in Figure 15 to Re(ξ) − c p plane, as shown in Figure 16. As a result, the phase velocity C x and the converting point C can be easily marked and are suitable to compare.  Figure 17 shows the 2D phase velocity curves for the nonhomogeneous magnetoelectro-elastic half-space, the result obtained at the nonhomogeneous coefficient n = 1 with the applied material constants taken from Table 1, which illustrates that the phase velocity of the converting point is 7.1 km/s. Figures 18-20 present the phase velocity curves with different cases. Figure 18 describes the phase velocity curves same as Figure 17, but, for nonhomogeneous coefficient n = 2, Figure 19 gives the phase velocity curves same to as Figure 17, but excluding the magnetic constants, and Figure 20 excluding the electric constants. It can be seen that the phase velocity of converting point soars to 12.7 km/s when the nonhomogeneous coefficient n = 2, which illustrates that nonhomogeneous coefficient has significant effect on converting point. Figure 19 presents the frequency spectra in the nonhomogeneous piezoelectric half-space, in which the phase velocity of converting point is 6.75 km/s near to 7.1 km/s, so magnetic coefficients have little effect on the converting point. The case in Figure 20 has similar characteristics to the previous case with respect to the converting point, but the value of converting point goes down slightly instead of increasing. Therefore, changing nonhomogeneous coefficient is the most straightforward method of changing the phase velocity of converting point, compared with changing electric or magnetic coefficients.

Conclusions
The Laguerre orthogonal polynomial expansion approach was improved to investigate the complex Rayleigh waves propagating in the nonhomogeneous magneto-electro-elastic half-spaces. The proposed method can present not only the solution of propagation waves but also that of complex or evanescent waves. The following conclusions can be drawn according to the numerical examples: 1.
The developed method takes angular frequency as independent variable and constructs characteristic equation, furthermore transforms the problem of finding the solution of acoustic waves into an eigenvalue problem, and successfully obtains the real, imaginary and complex wave modes.

2.
In this configuration and material, although electric and magnetic properties have less effect on the cut-off frequencies of complex solutions than nonhomogeneous coefficient, magnetic properties have more effect on the ratio of the real part to the imaginary part of complex solutions than other material properties. 3.
The energies of mechanical, electric and magnetic fields will decay to zero after traveling to the same point in the thickness direction, but the attenuation trajectories of mechanical displacement are obviously different from those of electric and magnetic potential. In the propagation direction, displacement, electric, and magnetic potential have similar attenuation tendencies and converting process.

4.
Nonhomogeneous coefficient has significant effect on the phase velocity of complex converting point, but the electric and magnetic properties do not, which, on the contrary, make the velocity slightly smaller.
The proposed method can be extended to other more complicated nonhomogeneous conditions with ease, even though here employed only one nonhomogeneous function for the convenience of expression. In practice, the function can be any derivative continuous function, so giving distinct nonhomogeneous function, respectively, to different material coefficient at the same time will make the present method more universe. Moreover, the propagation characteristics and behaviors of Rayleigh waves in layered half-space play an important role in nondestructive evaluation, ocean acoustics, seismic prospecting, etc., so further attention will be paid on the complex surface waves propagating in a layered half-space.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.