Size effects on the mixed modes and defect modes for a nano-scale phononic crystal slab

The size-dependent band structure of an Si phononic crystal (PnC) slab with an air hole is studied by utilizing the non-classic wave equations of the nonlocal strain gradient theory (NSGT). The three-dimensional (3D) non-classic wave equations for the anisotropic material are derived according to the differential form of the NSGT. Based on the the general form of partial differential equation modules in COMSOL, a method is proposed to solve the non-classic wave equations. The bands of the in-plane modes and mixed modes are identified. The in-plane size effect and thickness effect on the band structure of the PnC slab are compared. It is found that the thickness effect only acts on the mixed modes. The relative width of the band gap is widened by the thickness effect. The effects of the geometric parameters on the thickness effect of the mixed modes are further studied, and a defect is introduced to the PnC supercell to reveal the influence of the size effects with stiffness-softening and stiffness-hardening on the defect modes. This study paves the way for studying and designing PnC slabs at nano-scale.


Introduction
Phononic crystal (PnC) holds tremendous promise for novel devices and applications due to the characteristics of controlling the propagation of acoustic waves, especially the band gap property. It has been reported that PnC can be made as wave detectors [1] , waveguides [2] , filters [3][4] , transducers [5] , acoustic lenses [6] , acoustic focusing and imaging [7][8] , vibration isolators [9][10] , acoustic tweezers [11] , Klein tunneling [12] , and negative refraction [13][14] . With the development of modern communication techniques, acoustic waves at gigahertz (GHz) and even terahertz (THz) frequencies have become an essential aspect of research, e.g., acoustooptic coupling [15][16] , heat phonon propagation in thermoelectric devices [17] , acoustic radiation shield [18][19] , and phonon laser [20][21][22] . Consequently, the dimensions of these acoustic structures have to approach nano-scale. Moreover, the present nanotechnology makes it possible to fabricate and measure micro-or nano-scale PnC [23][24][25] . Due to the high surface-to-volume ratio of these nano-structures, the in-plane size effect or surface/interface effect becomes nonnegligible [26][27] . As we all know, the size effect cannot be evaluated by the classical elastic (CE) continuum theory. In order to overcome the weakness of the CE continuum theory, several non-classical continuum theories related to size-dependence were developed, e.g., the nonlocal elastic theory [28] , the strain gradient theory [29] , the surface elastic theory [30] , and the couple stress theory [31] . So far, both the in-plane size effect and the surface/interface effect have been studied extensively. Chen and Wang [32] , Chen et al. [33] , and Yan et al. [34] studied the band structures of nano-scale one-dimensional (1D) multilayer PnC with the transfer matrix method. The results showed that the nonlocal elastic continuum solution deviated from the CE continuum one finally approached the first-principles result as the thickness of each individual layer decreased [32] , and a cut-off frequency was found, beyond which the elastic wave could not propagate when the size effect was considered [33][34] . The recently developed meshless methods [35][36][37][38] provide high computational efficiency for calculating the band structure of PnC. Based on the nonlocal elasticity (NLE) theory, Zheng et al. [39] studied the size effect of nano-scale PnC with the meshless collocation method, and verified the results of Chen and Wang [32] . Hereafter, the surface elastic theory was used to calculate the band structure for two-dimensional (2D) PnC with air holes [40][41][42][43][44] . Zhang et al. [45] studied the surface effect of magneto-elastic PnC according to the Kirchhoff plate theory and surface elastic theory. Zhang et al. [46] and Zhang and Gao [47] studied the band gaps of 2D and three-dimensional (3D) PnC with the plane wave expansion method. In Ref. [48], the in-plane size effect on 2D PnC was investigated based on the nonlocal strain gradient theory (NSGT) with the finite element method (FEM), and both stiffness enhancement and soften behaviors were witnessed.
For PnC devices with slab structures, the thickness is a critical size parameter, especially for hole slabs. The corresponding band gap is opened when the thickness is less than the lattice constant. However, this phenomenon cannot be predicted by classical theories. Moreover, up to now, no research has been reported to study the thickness effect of the PnC slab. Therefore, in order to comprehensively understand the mechanism of size effects on the properties of the PnC slab, it is urgently necessary to study its thickness effect. Nano-scale PnC defect modes play an important role in cavity optomechanics [49] , and have received much attention for tremendous promising phenomena such as laser cooling nano-resonator to its ground state [50] and electromagnetically induced transparency [51] . However, there is no attention paid to the size effects on the defect modes. Thus, the research of the size effects on the defect modes is essential for promoting the development in cavity optomechanics.
In this study, the band structure of a nano-scale PnC slab is investigated based on the NSGT with the FEM. Combined with the constitutive equation for anisotropic material, the 3D nonclassic wave equations of the PnC slab are obtained. The non-classic wave equations are solved by using the general form partial differential equation module in COMSOL. The band structures of a circle hole PnC slab are obtained. First, the size-dependent behavior of the band structure is investigated, and the dependence of the first band gap upon the nonlocal parameters is focused. Then, the geometric parameters affecting the thickness effect are determined. Finally, the size effects on the defect modes are studied. at nano-scale. As a higher-order theory, the NSGT meets this requirement [55] . Moreover, the NSGT can be reduced to the former two theories and even the CE theory by fixing one or two nonlocal parameters as zero. Recently, the equation of motion was derived based on the NSGT, which successfully described the acoustic propagation behavior in 2D nano-PnC [48] . The equation of motion has the form as follows: where i, j = 1, 2, 3. f i and u i represent the body force and the displacement vector, respectively. ρ denotes the mass density. The total stress T ij is defined by where t ij and t ijm denote the nonlocal stress tensor and the nonlocal higher-order strain gradient stress tensor, respectively. The boundary conditions are where the traction t is prescribed on the boundary Ω T . The displacement u i or displacement gradient projection u (1) i in the outward normal direction of the surface is prescribed on the boundary Ω u . For the acoustic wave propagation in the anisotropic material, without external work, the harmonic wave equation can be expressed as follows: where C ijkl denotes the stiffness tensor, ε kl denotes the strain tensor, and ω denotes the angular frequency. l and ξ represent the nonlocal parameters, i.e., the length scale parameters. ∇ denotes the Hamiltonian operator. The detailed expansion of the wave equation is given in Appendix A.
Due to the introduction of the strain gradient, the wave equation becomes a system of 4th-order partial differential equations. However, the partial derivatives, which can be directly expressed in COMSOL, are at most second order. For order reduction, 27 intermediate variables are defined to represent the 2nd-order partial derivatives as follows: All 4th-order partial derivatives of displacements can be represented by the 2nd-order partial derivatives of these intermediate variables. In Ref. [48], the method of solving the non-classic wave equation based on COMSOL was verified, the 2D PnC was further calculated without considering the size effect, and the band structure agreed well with that obtained by the Dirichletto-Neumann map [34] .

Numerical results and discussion
Si is chosen as the material for the PnC slab. It is also a common optical material. Hence, the structures have the potential for fabricating phoxonic crystal (PxC) exhibiting both phononic and photonic band gaps simultaneously [16,22,56] . The relevant material parameters for Si are ρ = 2 331 kg · m −3 , C 11 = 165.7 GPa, C 12 = 63.9 GPa, C 44 = 79.62 GPa, = ωa/(2πc t ), and c t = 5 844 m · s, where ρ is the mass density, C 11 , C 12 , and C 44 are the elastic coefficients, is the normalized frequency, and c t is the transverse wave velocity. It is difficult to obtain the analytical solution for nonlocal strain gradient wave problems. However, with the general form of partial differential equation modules, the established non-classic wave equation of the NSGT in COMSOL can be easily solved [48] . First, the thickness effect on the PnC slab is studied without considering the strain gradient. Then, investigation is focused on the thickness effect on the band structure, as well as the band gap of the PnC slab with variable geometric and nonlocal parameters. In addition, for the PnC supercell slab with a point defect, the size effects on the defect modes are revealed.

Modal classification according to the energy ratio and symmetry
Consider an Si PnC slab of the square lattice with a circular hole. The lattice constant a = 10 nm, the radius of the circular hole r = 4.8 nm, corresponding to the filling ratio f = πr 2 /a 2 = 72%. The slab thickness is fixed as t = 3 nm. For this structure, the band gaps are shown in Fig. 1. Calculated by the CE theory, the normalized frequency ωa/(2πc t ) of the band gap ranges from 0.396 5 to 0.465 0. The in-plane mode and mixed mode can be classified by their bands with respect to the polarization [57] . The kinetic energy ratio in the If e = 0, the band is an in-plane mode; otherwise, it is a mixed mode containing both inplane and out-of-plane movements. In Fig. 1(a), the kinetic energy ratio e in the x 3 -direction is represented by the color of the band. The blue curves denote the in-plane modes, while the red curves represent the mixed modes. Next, attention is paid on the band gap and its edges. Figure 1(b) shows the mode shape of the edges (M 1 and M 2 ) and the possible edge (M 3 ) of the band gap. The modes M 1 , M 2 , and M 3 are located on the points Γ, X, and M of the irreducible Brillouin zone, respectively. The deformation of M 1 mainly concentrates on the x 3 -direction as the shearing motion. M 2 corresponds to the torsional motion in three directions. For M 3 , the deformation is mainly in-plane as the breathing motion. Notably, the band gap is determined by the mixed bands with vibration, mainly along the x 3 -direction. As shown in Fig. 1(c), the modes of the bands can be classified to odd modes and even modes according to the symmetry with respect to the middle plane, i.e., the plane of x 3 = 0 . The even and odd symmetry modes are calculated by applying symmetric and antisymmetric boundary conditions to the middle plane of the PnC slab, respectively. The even and odd band gaps are denoted by blue and red shadow zones, respectively. Since the red shadow zone is completely contained in the blue one, the odd band gap is just the complete band gap. Moreover, in the normalized frequency range of [0, 0.8], the in-plane modes in Fig. 1(a) correspond to the even modes in Fig. 1(c), while the mixed modes in Fig. 1(a) correspond to the odd modes in Fig. 1(c). The root cause is that the symmetric boundary condition leads the displacement u 3 to zero at the boundary, while the antisymmetric boundary condition causes the in-plane displacements u 1 and u 2 both to be zero at the boundary.

Comparison of the thickness effect and in-plane size effect
In the previous research on the size effect of the PnC, little attention has been paid to the thickness effect. Here, the size effect on the PnC slab is investigated by the proposed non-classic model. Note that 30 variables are required to solve the NSGT wave equation in COMSOL. For convenience, only the strain-driven nonlocal effect is considered. Let l = 0 in Eq. (4). Then, the NSGT is reduced to the NLE theory. It should be pointed out that just with one nonlocal parameter, Eringen's nonlocal theory fails to simultaneously fit the longitudinal and transverse acoustic dispersion curves of some materials, e.g., Si, Au, and Pt [58] . Therefore, the nonlocal parameters can be different for longitudinal and transverse acoustic waves. To distinguish the thickness effect from the in-plane size effect, the nonlocal parameter ξ is defined as ξ 1 and ξ 2 , which are perpendicular or parallel to the x 3 -direction. Thus, the governing equations are where ξ 1 and ξ 2 correspond to the in-plane size effect and the thickness effect, respectively. The wave equation can be reduced to the CE wave equation when ξ 1 = ξ 2 = 0. Besides, the determination of the nonlocal parameters is beyond the scope of the present study. For convenience, these parameters are chosen within the reasonable range [59] . Based on the NLE theory, the band structure is obtained by considering the size effect (see the dotted curves in Fig. 2). For comparison, the band structure obtained by the CE theory is plotted by the solid curves. In Fig. 2, the upper row is the bands for the even symmetry, and the lower row is the bands for the odd symmetry. The complete band gap is the intersection of two kinds of band gaps with even symmetry and odd symmetry. The normalized frequency ωa/(2πc t ) of the complete band gap ranges [0.392 5, 0.420 0], [0.360 9, 0.425 0], and [0.358 3, 0.419 7] by overlapping Figs. 2(a), 2(d); 2(b), 2(e); and 2(c), 2(f), respectively. The thickness effect and the in-plane size effect are discussed separately, and the combined effect based on the two effects is then discussed.
(I) As a total, both effects are considered with ξ 1 = 1 nm and ξ 2 = 1 nm, as shown in Figs. 2(a) and 2(d). The frequency of all bands drops, and the band gap narrows.
(II) Only the thickness effect is considered with ξ 1 = 0 nm and ξ 2 = 1 nm, as shown in Figs. 2(b) and 2(e). All bands of even symmetry nearly stay still (see Fig. 2(b)), while the bands of odd symmetry drop (see Fig. 2(e)). (III) Only the in-plane size effect is considered with ξ 1 = 1 nm and ξ 2 = 0 nm, as shown in Figs. 2(c) and 2(f). In Fig. 2(c), the in-plane mode is an even symmetric mode, and the frequency of its bands above the band gap decreases greatly and becomes the cutoff frequency marked by a red dot finally. As a result, the band gaps become much narrower.
In brief, as even modes, the in-plane modes are mainly affected by the in-plane size effect. Similarly, as odd modes, the mixed modes are dominated by the thickness size effect.

Size effect on the complete band gap
The complete band gaps are of our great interest, since in the complete band gap, the acoustic wave with any polarization and any direction cannot propagate through the PnC. In order to reveal the impact of the thickness effect on the complete band gap, the first complete band gaps are calculated for varying nonlocal parameters. As shown in Fig. 3(a), the starting frequency f s and cutoff frequency f c of the band gap both reduce as the nonlocal parameters ξ 1 and ξ 2 both increase. This is because both the in-plane size effect and the thickness effect cause stiffness-softening.
The gap-to-midgap ratio is defined as R BG = 200%(f c − f s )/(f c + f s ) to measure the relative width of the band gap. It is noteworthy that f c and R BG both decrease drastically when (ξ 1 , ξ 2 ) approaches (1, 0). This is because even with an equal nonlocal parameter, the thickness effect and in-plane effect have different impacts on the modes, corresponding to f s and f c , respectively. By fixing ξ 2 = 0 nm, when ξ 1 increases from 0 to 1, the frequency of the in-plane mode M 3 decreases drastically, and becomes the cutoff frequency of the band gap. Therefore, the complete band gap is affected more by the in-plane effect than the thickness effect.

Size-dependent behavior
It has been demonstrated that the thickness effect on the band structure of the PnC slab cannot be ignored. The thickness effect is further focused in terms of the geometric parameters. Only the modes with odd symmetry are considered since the thickness effect only acts on the mixed modes. In these parameters, the lattice constant a and filling ratio f are in-plane parameters, and f increases with the decease in the radius r of the hole. Figure 4 depicts the band structure comparison between the CE theory and NLE theory with the thickness effect, where the holes of the PnC slabs are with different r. No matter whether the thickness effect is considered or not, the band gap only appears when r = 0.48a, while the band gap disappears with the decrease in f . These bands of the mixed modes all shift down when only the thickness effect is considered. However, as r increases, the frequency shift due to the thickness effect does not change significantly. Thus, f has little influence on the thickness effect. Nevertheless, f has a great size effect on the in-plane modes of the circular hole PnC [48] . Fig. 4 Band structures of the mixed modes obtained by the CE theory and NLE theory under various r, where ξ1 = 0 nm, and ξ2 = 1 nm (color online) Figure 5 displays the band structure comparison between the CE theory and NLE theory for the lattice constants a = 10 nm, 20 nm, 50 nm, where f = 72% and t = 3 nm. The band structures obtained by the CE theory of these three PnCs are totally different. When a increases from 10 nm to 50 nm, the first band gap between the 3rd and 4th bands is gradually closed. However, as shown in Fig. 5(c), two new band gaps arise between the 8th and 9th bands and the 10th and 11th bands, respectively. Furthermore, comparing the band structure calculated by the CE theory and NLE theory, almost all band frequencies of the NLE theory drop from those of the CE theory for different a. Accordingly, the frequencies Fig. 5 Band structures of the mixed modes for different a, where the red shadow zones denote the band gaps obtained by the CE theory, the green shadow zones denote the band gaps obtained by the NLE theory, ξ1 = 0 nm, and ξ2 = 1 nm (color online) of the first band gaps decrease, i.e., the red area descends to the green area. This is also a softening phenomenon due to the thickness effect. Moreover, for the PnC slab, the smaller the lattice constant is, the larger its frequency drops, and thus the stronger the thickness effect is. Finally, the effect of the thickness t on the band structure is further studied. To ensure that the air hole PnC slab has a band gap, t is designed to be smaller than a. t is chosen as 0.3a, 0.6a, and 0.9a, respectively, as shown in Fig. 6. The first band gaps obtained by the CE theory appear for the PnC slabs with all t, while the second band gaps obtained by the CE theory only arise for the PnC slabs with two larger t. The thickness effect is then taken into consideration, for t = 0.6a in Fig. 6(b), and the first band gap is widened since the starting frequency shifts down while the cutoff frequency remains unchanged. Besides, the second band gap is closed. For t = 0.9a in Fig. 6(c), the thickness effect can be ignored when only the first two band gaps are concerned. Therefore, for PnC slabs, the thinner the thickness is, the stronger the thickness effect is. Fig. 6 Band structures of the mixed modes for different t, where the red shadow zones denote the band gaps obtained by the CE theory, the green shadow zones denote the band gaps obtained by the NLE theory, ξ1 = 0 nm, and ξ2 = 1 nm (color online)

Size effect on the defect modes
PxC can own phononic and photonic band gaps simultaneously. The defect PxC confines sound and light waves to defects, in which the acousto-optic coupling or optomechanical interaction is enhanced [56,[60][61][62] . However, so far, the research on the size effect of defect PnC has not been reported. Consequently, a point defect is introduced into a 7 × 7 PnC supercell to study the size effect on the defect modes (see Fig. 7(b)). As shown in Fig. 7(a), the bands of even and odd symmetries are plotted by blue circles and red dotted curves, respectively. Six flat bands of the defect modes, i.e., α, β, γ, σ, , and ζ, are found in the band gaps. They all belong to in-plane modes due to the even symmetry. These modes have an obvious feature, i.e., the displacement field is highly concentrated in the defect cavity (see Fig. 7(b)). Among them, the modes ε and ζ are a pair of degenerate modes sharing the same frequency, and their mode shapes can be obtained by rotating each other 90 degrees clockwise or counterclockwise. Based on the NSGT, the frequency shifts of these six defect modes vary with different nonlocal parameters (see Fig. 8). It should be pointed out that each nonlocal parameter is considered by fixing the others to zero. The frequency shifts of these modes are all proportional to the square of each nonlocal parameter. It can be noted that the frequency shift caused by ξ 1 is two orders of magnitude larger than that caused by ξ 2 . This is because the displacements of the defect modes are all concentrated in the x 1 x 2 -plane. Moreover, the degenerate modes and ζ have the same frequency shift caused by all three nonlocal parameters. In Fig. 8(a), although the frequency is the lowest among the six modes, α has the third highest frequency shift value. In Fig. 8(b), the frequency and its shift of these defect modes have different orders of magnitude. The order of the frequency shift is δ, β, α, γ, and (ζ) from small to large, while the order of the frequency is α, β, γ, δ, and (ζ) from small to large. The similar phenomenon can also be witnessed in Fig. 8 Figure 9 is further plotted to discuss the softening and hardening effects on the defect modes. The normalized frequency shifts ∆ωa/(2πc t ) of the six defect modes versus ξ and l are displayed, respectively. The normalized frequency shifts of the defect modes are represented by the color and contours. These contours are all parallel, indicating that the frequency shifts are all proportional to the square of ξ, as well as l. Besides, a contour with zero frequency shift is named as the neutral line. Above the neutral line, the modes all exert the stiffness-hardening effect. On the contrary, below the line, the modes all exert the stiffness-softening effect. At the neutral line, no size effect produces on the modes. Each frequency shift contour indicates that the frequency shifts are equal, which is the result of the joint action of stiffness hardening and softening. As shown in Fig. 9, for the six defect modes, the slopes of frequency shift contours are 0.01 times of 0.514, 0.235, 0.242, 0.386, 0.496, and 0.496, respectively. Consequently, even with the same nonlocal parameters, these modes exert different frequency shifts due to the size effect. It is totally different from the 1D beam NSGT wave problem, whose slopes of the frequency shift contours are 1 for all modes [63][64] .

Conclusions
A theoretical model for predicting the band structure of elastic wave propagation in a PnC slab is proposed based on the NSGT. The 3D wave equations are established to evaluate the size effect on the band gap. The modes are divided into in-plane modes and mixed modes according to the ratio of the kinetic energy contents in the x 3 -direction. Some conclusions can be drawn.
(i) The in-plane modes and mixed modes are obtained by applying symmetric and antisymmetric boundary conditions, respectively.
(ii) The frequency shifts of the in-plane modes and mixed modes are mainly dominated by the in-plane effect and the thickness size effect, respectively.
(iii) The strength of the thickness effect is stronger for smaller lattice constant or thickness.
(iv) For the complete band gap, the in-plane size effect is greater than the thickness size effect.
(v) The frequency shift of the defect modes caused by the in-plane nonlocal parameter is two orders of magnitude larger than that caused by the thickness nonlocal parameter since the defect modes are in-plane modes.
(vi) The frequency shifts of the defect modes are all proportional to the square of each stiffness-softening nonlocal parameter and stiffness-hardening nonlocal parameter. The contours of the normalized frequency shift versus two kinds of nonlocal parameters are all parallel. But different defect modes have different slopes of frequency shift contours.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.