## 2 Materials and Methods
### Analytical Approach
#### 2.2.2 Band Structure Solution  
To calculate the band structure of the piezoelectric phononic crystal with arbitrary cavities, it is necessary to discretize by using the finite element method. For the displacement of any point in the element [35]:  
\[u= \begin{bmatrix}u_{x}\\ u_{y}\\ u_{z}\end{bmatrix} =\mathbf{N}_{u}\mathbf{u}^{e}, \tag{2}\]  
where \(\mathbf{N}_{u}\) is the shape function matrix of the displacement. According to general knowledge of elasticity, the strain at any point in the unit can be defined by:  
\[S=\begin{bmatrix}\frac{\partial u_{x}}{\  
where \(\mathbf{N}_{\varphi}\) is the shape function of the potential. It is understood from general knowledge of electrostatics that the expression of the electric field at any point in the unit can be written as:  
\[E= \begin{bmatrix}E_{x}\\ E_{y}\\ E_{z}\end{bmatrix} = \begin{bmatrix}-\frac{\partial\varphi}{\partial x}\\ -\frac{\partial\varphi}{\partial y}\\ -\frac{\partial\varphi}{\partial x}\end{bmatrix} =- \begin{bmatrix}\frac{ \partial}{\partial x}\\ \frac{\partial}{\partial y}\\ \frac{\partial}{\partial z}\end{bmatrix} \mathbf{N}_{\varphi}\mathbf{\varphi} ^{\varepsilon}=\mathbf{B}_{\varphi}\mathbf{\varphi}^{\varepsilon}. \tag{7}\]  
In Equation (7), there is:  
\[\mathbf{B}_{\varphi}=- \begin{bmatrix}\frac{\partial}{\partial x}\\ \frac{\partial}{\partial y}\\ \frac{\partial}{\partial z}\end{bmatrix} \mathbf{N}_{\varphi}, \tag{8}\]  
therefore, the kinetic energy of the unit can be represented by:  
\[T^{\varepsilon}=\int\limits_{V_{1}}\frac{1}{2}\rho_{1}\dot{u}^{\mathrm{T}} \dot{u}dV+\int\limits_{V_{2}}\frac{1}{2}\rho_{2}\dot{u}^{\mathrm{T}}\dot{u}dV, \tag{9}\]  
where \(\rho_{1}\) and \(\rho_{2}\) represent the densities of non-piezoelectric and piezoelectric materials, while \(V_{1}\) and \(V_{2}\) represent the volumes of non-piezoelectric and piezoelectric materials. The potential energy of the unit can be written as:  
\[U^{\varepsilon}=\int\limits_{V_{1}}\frac{1}{2}S^{\mathrm{T}}TdV+\int\limits_{ V_{2}}\frac{1}{2}S^{\mathrm{T}}TdV-\int\limits_{V_{2}}\frac{1}{2}EDdV. \tag{10}\]  
According to [PERSON]'s principle, it can be seen that:  
\[\delta\int\left(T^{\varepsilon}-U^{\varepsilon}\right)dt=0, \tag{11}\]  
then Equations (9) and (10) are brought into Equation (11) for variational calculations. Thus, the finite element equation of motion of the unit can be acquired:  
\[\mathbf{M}_{uu}^{\varepsilon}\ddot{\mathbf{u}}^{\varepsilon}+\mathbf{K}_{uu} ^{\varepsilon}\mathbf{u}^{\varepsilon}+\mathbf{K}_{u\varphi}^{\varepsilon} \mathbf{\varphi}^{\varepsilon}=0, \tag{12}\]  
\[\mathbf{K}_{qu}^{\varepsilon}\mathbf{u}^{\varepsilon}+\mathbf{K}_{\varphi \varphi}^{\varepsilon}\mathbf{\varphi}^{\varepsilon}=0. \tag{13}\]  
In the above formula, \(\mathbf{M}_{uu}^{\varepsilon}\) and \(\mathbf{K}_{uu}^{\varepsilon}\) represent the mass matrix and the elastic stiffness matrix of the unit, respectively. \(\mathbf{K}_{u\varphi}^{\varepsilon}\) and \(\mathbf{K}_{\varphi u}^{\varepsilon}\) are the piezoelectric-elastic coupling stiffness matrices, and they are transposed to each other. The expression of each unit matrix is:  
\[\mathbf{M}_{uu}^{\varepsilon}=\int\limits_{V}\mathbf{N}_{u}^{\mathrm{T}} \rho\mathbf{N}_{u}dV, \tag{14}\]  
\[\mathbf{K}_{uu}^{\varepsilon}=\int\limits_{V}\mathbf{B}_{u}^{\mathrm{T}} \mathbf{c}\mathbf{B}_{u}dV, \tag{15}\]  
\[\mathbf{K}_{u\varphi}^{\varepsilon}=\int\limits_{V}\mathbf{B}_{u}^{\mathrm{T}} \mathbf{c}\mathbf{B}_{\varphi}dV, \tag{16}\]  
\[\mathbf{K}_{\varphi\varphi}^{\varepsilon}=\int\limits_{V}\mathbf{B}_{\varphi }^{\mathrm{T}}\mathbf{c}\mathbf{B}_{\varphi}dV, \tag{17}\]where \(\mathbf{c}\) is the elasticity coefficient matrix of the unit, \(\mathbf{e}\) is the piezoelectric strain constant matrix, and \(\mathbf{\varepsilon}\) is the dielectric constant matrix.