Localization of elastic waves in one-dimensional detuned phononic crystals with flexoelectric effect

ABSTRACT

Although forbidden band effect in perfectly periodic phononic crystals (PC) is very attractive, random disordered (i.e. detuning) phenomenon is inevitable in engineering processing, thus exploring the effect of detuning on the wave characteristics of PC becomes a necessity. In this study, fundamental governing equations and boundary conditions are derived from the principle of virtual work. Wave characteristics and localization factor of one-dimensional (1D) detuned nano-PC are investigated based on the transfer matrix method, with flexoelectric effect duly accounted for. Subsequently, with $\mathrm{B}\mathrm{a}\mathrm{T}\mathrm{i}{\mathrm{O}}_3-\mathrm{S}\mathrm{r}\mathrm{T}\mathrm{i}{\mathrm{O}}_3$ nano-PC taken for illustration, forbidden band properties and localization factor of 1D elastic waves in harmonic and detuned states are systematically characterized. It is demonstrated that localization factor can characterize the energy band structure of 1D PC perfectly. Flexoelectric effect tends to widen the width of forbidden band with increasing detuning, and detuning is linearly related to the bandwidth. The forbidden band is more sensitive to flexoelectric coefficient detuning than thickness detuning. The research results provide useful theoretical guidance for designing high-frequency nanoscaled devices with the function of filtering based on band gap effect of nano-PC.

Graphical abstract

KEYWORDS:

• Flexoelectric effect
• phononic crystal
• detuning
• localization factor
• forbidden band

1. Introduction

The flexoelectric effect belongs to electromechanical coupling, but it is not a simple higher-order extension of the piezoelectric effect. The piezoelectric effect acts only in non-centrosymmetric materials, while the flexoelectric effect exists in all dielectric materials. Flexoelectricity is a symmetry independent electromechanical coupling phenomenon that outperforms piezoelectricity at micro and nanoscales, due mainly to its size-dependent behavior arising from gradient terms in its constitutive relations [1].

The flexoelectric effect describes the coupling between strain gradients and electrical polarization. A multitude of investigations has been carried out to understand the phenomenon and explore its underlying physical mechanisms. Ever since a phenomenological framework [2] was proposed to characterize the flexoelectric effect, theoretical research on flexoelectric effect has developed rapidly. For example, the microscopic mechanism of flexoelectric effect was proposed to explain the polarization signals resulting from shock loading a nonpiezoelectric material [3]. A polarization gradient was added to the conventional elastic dielectric energy function to establish a theory to characterize the flexoelectric effect [4]; nonetheless, as the effect of strain gradient was not taken into account, the theory was not perfect and could only partially explain the experiment phenomenon. With the rapid development of flexoelectric theories, the flexoelectric coefficient of cubic ion crystal has been continuously calculated [5]. Indenbom et al. [6] put forward the Landau theory of flexoelectric effect, dividing formally the flexoelectric effect and the nonlocal piezoelectric effect [7]. More recently, Maranganti et al. [8] developed a mathematical framework containing strain gradient and polarization coupling, proposed Green’s function solution of the governing equation, and solved the inclusion problem, while Shen et al. [9,10] established the dielectric variational principle by considering flexoelectric effect, surface effect, and electrostatic force and derived the corresponding governing equations and boundary conditions: the proposed theory can not only solve complex static electromechanical coupling problems in nanodielectrics but also analyze dynamic problems such as elastic wave propagation when kinetic energy is introduced into the dielectric variational principle. In general, a large number of studies have been carried out on the structures related to piezoelectric materials, including nonlinearity and nanosize [11]. However, the research on flexoelectricity is just beginning.

Wave propagation in phononic crystals (PC) has become another hot topic, for the emergence of band gap characteristics provides new ideas for energy control, leading to many important discoveries. The essence of studying PC is to investigate the propagation of elastic waves in periodic materials or structures. In fact, systematic studies on the propagation characteristics of waves in periodic materials [12] have been carried out long ago, but it was not until 1995 that the existence of elastic wave band gaps was firstly determined from the perspective of experiments [13]. Before 2000, the study of PC was limited to Bragg scattering principle, which requires that the lattice constant and band gap frequency should be in the same order of magnitude [14,15]. Afterward, a new PC structure with local resonant property was creatively proposed [16]. The lattice constant of this crystal is 2 orders of magnitude smaller than the wavelength corresponding to the band gap, which marks another major breakthrough in the study of PC. Currently, PC with increasingly smaller sizes is used to design complete acoustic band gap at high frequencies, thus providing a vibration-free environment for high-precision mechanical systems [15]. In addition, the propagation and resonance of horizontal shear waves in layered PC with strip cracks were studied [17], and the results reveal wave resonance and localization in band gap. Experiment measurements [18] were also made on the complete surface acoustic band gap of a two-dimensional piezoelectric PC. It was found that the frequency range of complete band gap was exactly the same as that predicted by theory. Nanomaterials with various dimensionalities (e.g. nanowires, nanofilms, two-dimensional materials, and three-dimensional nanostructures) have shown great potential in recent development of flexible electronics [19] and PC. However, most of existing studies were carried out by assuming perfect periodic PC at macroscopic scale, without no consideration of nanoscale as well as the detuning effect.

In general, the study of PC at nanoscale with flexoelectric effect considered is relatively lacking. Liu et al. [20] studied the dispersion relation of one-dimensional (1D) PC considering flexoelectric effect by using the transfer matrix method, while Yang et al. [21–23] developed theoretical models to analyze elastic waves in nanolayer PC, love waves in piezoelectric nanocomposites, and Lamb waves in infinite nanoplates. Further, Eliseev et al. [24] proved that shear horizontal waves can propagate near the smooth surface of nonpiezoelectric crystals by considering the flexoelectric effect, whereas Hu et al. [25] explored the propagation characteristics of bulk wave in elastic media when flexoelectric effect, micro-inertia effect, and strain gradient elasticity were considered simultaneously. In addition, the study of wave propagation in complex conditions and structures also has a lot of reference value. The wave propagation of functionally graded ceramic-metal plates [26], sandwich plate [27], and sandwich structure with a soft core and multi-hybrid nano-composite face sheets [28] have been studied in depth. However, as these theoretical studies are all based on the premise that the PC has perfect periodic nanostructure, the effect of detuning on wave propagation in nano-PC with flexoelectric effect considered remains elusive. In practice, the detuning phenomenon is inevitable during the processing of periodic structures. Moreover, the band structure of PC can be adjusted by deliberately introducing the detuning, which in turn controls the propagation behavior of elastic waves. It is therefore particularly important to study the wave propagation properties of detuned flexoelectric nano-PC so that theoretical guidance for practical processing can be obtained.

This paper focused on studying theoretically the propagation properties of elastic waves in harmonic and detuned nano-PC with idealized 1D layered structures. The flexoelectric theory is employed to establish the governing equation of wave propagation, with the effects of strain gradient and polarization gradient accounted for. The transfer matrix method and the Bloch theory are subsequently used to calculate the localization factor in the detuning state, which is then used to characterize the band structure and localization phenomenon of the detuned 1D PC. The influence of layer thickness detuning and flexoelectric coefficient detuning on energy band structure is quantified, and the underlying physical mechanisms are explored.

2. Problem statement

To study the wave characteristics of flexoelectric PC in detuning state, an idealized layered PC is selected for theoretical analysis, as shown in Figure 1. The 1D PC consisted of alternating sublayers A and B that are made of different materials, with layer thicknesses $d_1$ and $d_2$, respectively. Its unit cell has thence a thickness of $d=d_1+d_2$, i.e. the periodic thickness of the PC. In the present study, the thickness $d_1$ of sublayer A is not fixed but varied for random detuning.

Figure 1. Schematic of 1D layered PC.

To introduce the effect of flexoelectricity, a constructive and effective approach is extending the linear piezoelectric theory. For non-centrosymmetric dielectrics, the bulk internal energy density function relating to polarization $P_i$, polarization gradient $P_{i,j}$, small strain ${\epsilon }_{ij}$, and strain gradient $u_{j,kl}$ can be written as [10]:

(1) $U_{\mathrm{b}}=\frac{1}{2}{a}_{kl}{P}_k{P}_l+\frac{1}{2}{b}_{ijkl}{P}_{i,j}{P}_{k,l}+\frac{1}{2}{c}_{ijkl}{\epsilon }_{ij}{\epsilon }_{kl}+v_{ijk}{\epsilon }_{ij}{P}_k+e_{ijkl}{\epsilon }_{ij}{P}_{k,l}+f_{ijkl}{P}_i{u}_{j,kl}$(1)

Here, strain ${\epsilon }_{ij}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right)$, ($a\mathit{,}b\mathit{,}c\mathit{,}v\mathit{,}e\mathit{,}f$) are the material property tensors, $a$ represents the second-order dielectric tensor, $b$ is a high order dielectric constant, characterizing the coupling between the electric polarization gradients, $c$ is the fourth-order elastic tensor, $v$ is the third-order piezoelectric tensor, and $e$ was introduced by Mindlin [29] in his theory of polarization gradient to link the gradients of polarization to strains, while the tensor $f$ is related to strain gradient polarization coupling, its components commonly referred to as the flexoelectric coefficients [20,30].

To highlight the flexoelectric effect and neglect the influence of piezoelectric effect, it is assumed herein that all the subjects of the current study belong to central symmetric dielectrics, i.e. $v\mathit{=}\mathbf{0}$ (i.e. the third-order piezoelectric tensor $v$ is null in dielectrics with inversion symmetry). The bulk internal energy density function can thence be written as [10]:

(2) $U_{\mathrm{b}}=\frac{1}{2}{a}_{kl}{P}_k{P}_l+\frac{1}{2}{b}_{ijkl}{P}_{i,j}{P}_{k,l}+\frac{1}{2}{c}_{ijkl}{\epsilon }_{ij}{\epsilon }_{kl}+e_{ijkl}{\epsilon }_{ij}{P}_{k,l}+f_{ijkl}{P}_i{u}_{j,kl}$(2)

Under the assumption of infinitesimal deformation, the constitutive equation can be expressed as:

(3) ${\sigma }_{ij}=\frac{\mathrm{\partial }{U}_b}{\mathrm{\partial }{\epsilon }_{ij}}=c_{ijkl}{\epsilon }_{kl}+e_{ijkl}{P}_{k,l}$(3)

(4) ${\tau }_{ijm}=\frac{\mathrm{\partial }{U}_b}{\mathrm{\partial }{u}_{i,jm}}=f_{kijm}{P}_k$(4)

(5) $E_i=\frac{\mathrm{\partial }{U}_b}{\mathrm{\partial }{P}_i}={a_i}_j{P}_j+f_{ijkl}{u}_{j,kl}$(5)

(6) $V_{ij}=\frac{\mathrm{\partial }{U}_b}{\mathrm{\partial }{P}_{i,j}}=b_{ijkl}{P}_{k,l}+e_{klij}{u}_{k,l}$(6)

where ${\sigma }_{ij}$ is the stress tensor in classical elasticity, ${\tau }_{ijm}$ is the higher-order stress (couple stress) tensor, $E_i$ is the effective local electric field, and $V_{ij}$ is the higher-order local electric field. Note that ${\sigma }_{ij}={\sigma }_{ji}$ and ${\tau }_{ijm}={\tau }_{jim}$.

To incorporate the micro-inertia effect into mathematical model, the following expression for the kinetic energy density is considered as [31,32]:

(7) $K=\frac{1}{2}\rho {\dot{u}}_i{\dot{u}}_i+\frac{1}{2}\rho l^2{\ddot{u}}_{i,j}{\ddot{u}}_{i,j}$(7)

where $\rho$ is the mass density, $l$ is defined as the micro-inertia characteristic length (the scaling parameter for dynamics), and the dot over the vector component $u_i$ refers to time derivative.

Governing equations and boundary conditions are derived from the virtual work principle, which has the following form:

(8) $-\underset{\mathrm{\Omega }}{\int }\delta U_{b}d\mathrm{\Omega }+\underset{\mathrm{\Omega }}{\int }Kd\mathrm{\Omega }=0$(8)

where $\delta$ is the variational notation, $\mathrm{\Omega }$ is the volume occupied by the dielectric material.

Applying the virtual work principle leads to the following governing equations [20](20) $c\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}+h\frac{{\mathrm{\partial }}^{2}P}{\mathrm{\partial }{x}^2}=\rho \frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{t}^2}-\rho l^2\frac{{\mathrm{\partial }}^{4}u}{\mathrm{\partial }{x}^2\mathrm{\partial }{t}^2}$(20) :

(9) $({\sigma }_{ij}-{\tau }_{ijm,m}+\rho l^2{\ddot{u}}_{i,j}{)}_{,j}+F_i=\rho {\ddot{u}}_i$(9)

(10) $E_i-V_{ij,j}+{\phi }_{,i}=0$(10)

(11) $-{\epsilon }_0{\phi }_{,ii}+P_{i,i}=0$(11)

where ${\epsilon }_0$ represents the dielectric constant in vacuum, $F$ represents the external force, and $\phi$ is the potential of Maxwell self-field defined as:

(12) $E_i=-{\phi }_{,i}$(12)

Since ${\sigma }_{ij}-{\tau }_{ijm,m}$ appears in the force balance relation, it can be understood as a “physical stress” [8,33], namely:

(13) ${\sigma }_{ij}^{phys}={\sigma }_{ij}-{\tau }_{ijm,m}$(13)

Upon substituting the constitutive equations (3)–(3) ${\sigma }_{ij}=\frac{\mathrm{\partial }{U}_b}{\mathrm{\partial }{\epsilon }_{ij}}=c_{ijkl}{\epsilon }_{kl}+e_{ijkl}{P}_{k,l}$(3) (6(6) $V_{ij}=\frac{\mathrm{\partial }{U}_b}{\mathrm{\partial }{P}_{i,j}}=b_{ijkl}{P}_{k,l}+e_{klij}{u}_{k,l}$(6) ) into the governing equations (9)–(9) $({\sigma }_{ij}-{\tau }_{ijm,m}+\rho l^2{\ddot{u}}_{i,j}{)}_{,j}+F_i=\rho {\ddot{u}}_i$(9) (11(11) $-{\epsilon }_0{\phi }_{,ii}+P_{i,i}=0$(11) ). For 1D PC, the elastic field and electric fields vary only in the $x$-direction when a longitudinal wave is incident vertically, so that the inhomogeneous equations considering flexoelectric effect can be obtained as:

(14) $c{\mathrm{\nabla }}^{2}u+(e-f){\mathrm{\nabla }}^{2}P+\rho l^2{\mathrm{\nabla }}^{2}u+F=\rho u$(14)

(15) $(e-f){\mathrm{\nabla }}^{2}u+b{\mathrm{\nabla }}^{2}P-aP-\mathrm{\nabla }\phi =0$(15)

(16) $-{\epsilon }_0{\mathrm{\nabla }}^2\phi +\mathrm{\nabla }P=0$(16)

Since both tensor $e$ and tensor $f$ always appear in pairs, we define $h$ (i.e. $h=e-f$) instead of ($e-f$) for conciseness [8]. Neglect the influence of external force, so that Eqs. (14)–(14) $c{\mathrm{\nabla }}^{2}u+(e-f){\mathrm{\nabla }}^{2}P+\rho l^2{\mathrm{\nabla }}^{2}u+F=\rho u$(14) (16(16) $-{\epsilon }_0{\mathrm{\nabla }}^2\phi +\mathrm{\nabla }P=0$(16) ) can be reduced to inhomogeneous equations, as:

$c\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}+h\frac{{\mathrm{\partial }}^{2}P}{\mathrm{\partial }{x}^2}=\rho \frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{t}^2}-\rho l^2\frac{{\mathrm{\partial }}^{4}u}{\mathrm{\partial }{x}^2\mathrm{\partial }{t}^2}$ (17)

(18) $h\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}+b\frac{{\mathrm{\partial }}^{2}P}{\mathrm{\partial }{x}^2}-aP-\frac{\mathrm{\partial }\phi }{\mathrm{\partial }x}=0$(18)

(19) $\frac{\mathrm{\partial }P}{\mathrm{\partial }x}-{\epsilon }_0\frac{{\mathrm{\partial }}^2\phi }{\mathrm{\partial }{x}^2}=0$(19)

Under the condition of open circuit that $P-{\epsilon }_0(\mathrm{\partial }\phi /\mathrm{\partial }x)=0$, Eq. (19)(19) $\frac{\mathrm{\partial }P}{\mathrm{\partial }x}-{\epsilon }_0\frac{{\mathrm{\partial }}^2\phi }{\mathrm{\partial }{x}^2}=0$(19) is satisfied. Then, the governing equations are

(20) $c\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}+h\frac{{\mathrm{\partial }}^{2}P}{\mathrm{\partial }{x}^2}=\rho \frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{t}^2}-\rho l^2\frac{{\mathrm{\partial }}^{4}u}{\mathrm{\partial }{x}^2\mathrm{\partial }{t}^2}$(20)

(21) $h\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}+b\frac{{\mathrm{\partial }}^{2}P}{\mathrm{\partial }{x}^2}-\eta P=0$(21)

Where $\eta =a+{\epsilon }_0^{-1}$.

Assume that the vertical incident longitudinal wave is steady-state harmonic, i.e. its displacement and potential vary with time harmonically. The solution to the above system of equations (20)(20) $c\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}+h\frac{{\mathrm{\partial }}^{2}P}{\mathrm{\partial }{x}^2}=\rho \frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{t}^2}-\rho l^2\frac{{\mathrm{\partial }}^{4}u}{\mathrm{\partial }{x}^2\mathrm{\partial }{t}^2}$(20) and (21)(21) $h\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}+b\frac{{\mathrm{\partial }}^{2}P}{\mathrm{\partial }{x}^2}-\eta P=0$(21) can thence be expressed as:

(22) $u(x,t)=\bar{U}(x)exp(-i\omega t)$(22)

(23) $P(x,t)=\bar{P}(x)exp(-i\omega t)$(23)

where $i=\sqrt{-1}$. To facilitate the calculation, the dimensionless coordinate, ${\zeta }_i=\frac{x_j}{{\bar{d}}_1}$, is introduced, where ${\bar{d}}_1$ represents the average thickness of the first sublayer. The wave equations (20)(20) $c\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}+h\frac{{\mathrm{\partial }}^{2}P}{\mathrm{\partial }{x}^2}=\rho \frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{t}^2}-\rho l^2\frac{{\mathrm{\partial }}^{4}u}{\mathrm{\partial }{x}^2\mathrm{\partial }{t}^2}$(20) and (21)(21) $h\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}+b\frac{{\mathrm{\partial }}^{2}P}{\mathrm{\partial }{x}^2}-\eta P=0$(21) can thence be rewritten as:

(24) $\left(c+\rho l^2{\omega }^2\right)\frac{{\mathrm{\partial }}^2\bar{U}(\zeta )}{\mathrm{\partial }{\zeta }^2}+h\frac{{\mathrm{\partial }}^2\bar{P}(\zeta )}{\mathrm{\partial }{\zeta }^2}+{{\bar{d}}_1}^2\rho {\omega }^2\bar{U}(\zeta )=0$(24)

(25) $h\frac{{\mathrm{\partial }}^2\bar{U}(\zeta )}{\mathrm{\partial }{\zeta }^2}+b\frac{{\mathrm{\partial }}^2\bar{P}(\zeta )}{\mathrm{\partial }{\zeta }^2}-{{\bar{d}}_1}^2\eta \bar{P}(\zeta )=0$(25)

Detailed derivation process of Eqs. (24)(24) $\left(c+\rho l^2{\omega }^2\right)\frac{{\mathrm{\partial }}^2\bar{U}(\zeta )}{\mathrm{\partial }{\zeta }^2}+h\frac{{\mathrm{\partial }}^2\bar{P}(\zeta )}{\mathrm{\partial }{\zeta }^2}+{{\bar{d}}_1}^2\rho {\omega }^2\bar{U}(\zeta )=0$(24) and (25(25) $h\frac{{\mathrm{\partial }}^2\bar{U}(\zeta )}{\mathrm{\partial }{\zeta }^2}+b\frac{{\mathrm{\partial }}^2\bar{P}(\zeta )}{\mathrm{\partial }{\zeta }^2}-{{\bar{d}}_1}^2\eta \bar{P}(\zeta )=0$(25) ) is presented in the Appendix.

Eventually, the solution can be deduced by substituting Eqs. (22)(22) $u(x,t)=\bar{U}(x)exp(-i\omega t)$(22) and (23)(23) $P(x,t)=\bar{P}(x)exp(-i\omega t)$(23) into Eqs. (24)(24) $\left(c+\rho l^2{\omega }^2\right)\frac{{\mathrm{\partial }}^2\bar{U}(\zeta )}{\mathrm{\partial }{\zeta }^2}+h\frac{{\mathrm{\partial }}^2\bar{P}(\zeta )}{\mathrm{\partial }{\zeta }^2}+{{\bar{d}}_1}^2\rho {\omega }^2\bar{U}(\zeta )=0$(24) and (25(25) $h\frac{{\mathrm{\partial }}^2\bar{U}(\zeta )}{\mathrm{\partial }{\zeta }^2}+b\frac{{\mathrm{\partial }}^2\bar{P}(\zeta )}{\mathrm{\partial }{\zeta }^2}-{{\bar{d}}_1}^2\eta \bar{P}(\zeta )=0$(25) ), as:

(26) ${\bar{U}}_j(\zeta )=C_{1j}{e}^{l_{1j}\zeta }+C_{2j}{e}^{-l_{1j}\zeta }+C_{3j}{e}^{i{l}_{2j}\zeta }+C_{4j}{e}^{-i{l}_{2j}\zeta }$(26)

(27) $\begin{array}{c}{\bar{P}}_j(\zeta )=-\frac{1}{h_j}[C_{1j}(c_j+{\rho }_j{l}_j^2{\omega }^2+\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{1j}^2})e^{l_{1j}\zeta }+C_{2j}(c_j+{\rho }_j{l}_j^2{\omega }^2+\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{1j}^2})e^{-l_{1j}\zeta }\\ +C_{3j}(c_j+{\rho }_j{l}_j^2{\omega }^2-\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{2j}^2})e^{i{l}_{2j}\zeta }+C_{4j}(c_j+{\rho }_j{l}_j^2{\omega }^2-\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{2j}^2})e^{-i{l}_{2j}\zeta }]\end{array}$(27)

Here, $C_{1j}$, $C_{2j}$, $C_{3j}$, and $C_{4j}$ are undetermined coefficients, $j=1$ and $2$ represent the first and second sublayer, respectively, and $l_{1j}$ and $l_{2j}$ are given by:

(28) $l_{1j}=\sqrt{\frac{\left(c_j+{\rho }_j{l}_j^2{\omega }^2\right){\eta }_j{{\bar{d}}_1}^2-{{\bar{d}}_1}^2{\rho }_j{\omega }^2{b}_j+\sqrt{{\left(\left(c_j+{\rho }_j{l}_j^2{\omega }^2\right){\eta }_j{{\bar{d}}_1}^2+{{\bar{d}}_1}^2{\rho }_j{\omega }^2{b}_j\right)}^2-4{{\bar{d}}_1}^4{h}_j^2{\rho }_j{\omega }^2{\eta }_j}}{2\left(b_j\left(c_j+{\rho }_j{l}_j^2{\omega }^2\right)-h_j^2\right)}}$(28)

(29) $l_{2j}=\sqrt{\frac{2{\rho }_j{\omega }^2{\eta }_j{{\bar{d}}_1}^4}{\left(c_j+{\rho }_j{l}_j^2{\omega }^2\right){\eta }_j{{\bar{d}}_1}^2-{{\bar{d}}_1}^2{\rho }_j{\omega }^2{b}_j+\sqrt{{\left(\left(c_j+{\rho }_j{l}_j^2{\omega }^2\right){\eta }_j{{\bar{d}}_1}^2+{{\bar{d}}_1}^2{\rho }_j{\omega }^2{b}_j\right)}^2-4{{\bar{d}}_1}^4{b}_j^2{\rho }_j{\omega }^2{\eta }_j}}}$(29)

3. Transfer matrix method

The transfer matrix method is efficient in dealing with dynamic problems of 1D structure. In the present study, the transfer matrix between adjacent cells in the 1D PC can be derived from the continuity of mechanical and electrical boundary conditions. Specifically, in accordance with the definition in Ref [8,20], the physical stress ${\sigma }^{phys}$ and electric tensors $V^{phys}$ can be expressed, respectively, as:

(30) ${\sigma }^{phys}=c\frac{\mathrm{\partial }u}{\mathrm{\partial }\zeta }+h\frac{\mathrm{\partial }P}{\mathrm{\partial }\zeta }$(30)

(31) $V^{phys}=h\frac{\mathrm{\partial }u}{\mathrm{\partial }\zeta }+b\frac{\mathrm{\partial }P}{\mathrm{\partial }\zeta }$(31)

For the 1D case of Figure 1, according to the derivation of the transfer matrix in Ref [34], the state vector is given by:

(32) $v=M_{\mathrm{j}}(\zeta )\{\bar{U},\bar{P},{\sigma }^{phys},V^{phys}{\}}^T$(32)

It follows that the matrix equations at the left and right sides of the first sublayer in the $n\mathrm{t}\mathrm{h}$ unit cell are:

(33) $v_{1L}^{(n)}=M_1(0)\{C_{11},C_{21},C_{31},C_{41}{\}}^T$(33)

(34) $v_{1R}^{(n)}=M_1(d_1)\{C_{11},C_{21},C_{31},C_{41}{\}}^T$(34)

where the subscripts $L$ and $R$ represent the left and right sides of the sublayer, respectively. Similarly, the matrix equations at the left and right sides of the second sublayer in the $n\mathrm{t}\mathrm{h}$ unit cell are:

(35) $v_{2L}^{(n)}=M_2(d_1)\{C_{12},C_{22},C_{32},C_{42}{\}}^T$(35)

(36) $v_{2R}^{(n)}=M_2(d)\{C_{12},C_{22},C_{32},C_{42}{\}}^T$(36)

Continuity at the interface between adjacent sublayers dictates that:

$v_{1R}^{(n)}=v_{2L}^{(n)}$ (37)

From Eqs. (32)(32) $v=M_{\mathrm{j}}(\zeta )\{\bar{U},\bar{P},{\sigma }^{phys},V^{phys}{\}}^T$(32) –(37), the relationship between the left and right sides of the $n\mathrm{t}\mathrm{h}$ unit cell is obtained as:

(38) $v_{2R}^{(n)}=M_2(d)[M_2(d_1){]}^{-1}{M}_1(d_1)[M_1(0){]}^{-1}{v}_{1R}^{(n)}$(38)

It follows thence that the continuity condition between the right side of the $(n-1)\mathrm{t}\mathrm{h}$ unit cell and the left side of the $n\mathrm{t}\mathrm{h}$ unit cell is:

$v_{1L}^{(n)}=v_{2R}^{(n-1)}$ (39)

The relationship between the $(n-1)\mathrm{t}\mathrm{h}$ and $n\mathrm{t}\mathrm{h}$ unit cells is further derived as:

(40) $v_{2R}^{(n)}=M_2(d)[M_2(d_1){]}^{-1}{M}_1(d_1)[M_1(0){]}^{-1}{v}_{2R}^{(n-1)}$(40)

Therefore, $T_n=M_2(d)[M_2(d_1){]}^{-1}{M}_1(d_1)[M_1(0){]}^{-1}$ is the transfer matrix between two connected unit cells in the PC. Due to the transferability of the structure, ${\mathit{T}}_{\mathit{n}}$ remains unchanged for any $n$ when detuning is not considered and hence can be expressed as T. In other words, T can be viewed as the transfer matrix of the periodic PC, which is a function of sublayer thickness and frequency.

The problem of elastic wave propagation in a periodic system such as the one depicted in Figure 1 can be simplified to an eigenvalue problem by using the Bloch theory, as [35]:

(41) $\left|T-e^{ikd}I\right|=0$(41)

where $I$ is a $4\times 4$ unit matrix, $k$ is the Bloch wave number, and $d=d_1+d_2$ is the thickness of the unit cell. By calculating the frequency $\omega$ corresponding to all wavenumbers $k$ in the irreducible Brillouin zone, the dispersion curve of the 1D flexoelectric PC can be obtained.

4. Prediction of localization factor

Localization factor is an important concept introduced to describe vibration localization in near-periodic structures [36]. It describes the average exponential decay of fluctuation amplitude when an elastic wave propagates in a detuning periodic structure and can be used to characterize the influence of localization phenomenon on structural vibration and strength.

To calculate the localization factor, two approaches are commonly adopted: one is the Lyapunov method [37] and the other is the matrix eigenvalue method [38]. By contrast, the calculation with the method of matrix eigenvalue is fairly intuitive, while the calculation process of the Lyapunov method is more complex and cumbersome. For simplicity, the matrix eigenvalue method is employed in the present study.

With 1D detuned PC is assumed, the random variables can be chosen in a range of options, such as sublayer thickness detuning, elastic modulus detuning, density detuning, flexoelectric coefficient detuning, and so on. For the convenience of analysis, the random variable of the $i\mathrm{t}\mathrm{h}$ unit cell is noted as:

(42) ${\epsilon }_i=({\epsilon }_{i1},{\epsilon }_{i2}\cdots {\epsilon }_{iq})$(42)

where q represents the number of random variables. According to the derivation of the transfer matrix above, the relationship between the $i\mathrm{t}\mathrm{h}$ and the $(i+1)\mathrm{t}\mathrm{h}$ state vectors can be written as:

(43) $v_{i+1}=T({\epsilon }_i)v_i$(43)

where $T({\epsilon }_i)$ is the transfer matrix, which is related to the random variable ${\epsilon }_i$. For a randomly detuned PC, the localization factor is an important characterization parameter that can be calculated from the eigenvalues of the transfer matrix, as:

(44) ${\lambda }_l=\underset{N\to \mathrm{\infty }}{lim}\frac{1}{N}\sum _{i=1}^{N}ln\left|{\lambda }_{i,l}({\epsilon }_i)\right|$(44)

where ${\lambda }_{i,l}({\epsilon }_i)$ is the $l\mathrm{t}\mathrm{h}$ eigenvalue of the transfer matrix $T({\epsilon }_i)$ of the $i\mathrm{t}\mathrm{h}$ unit cell. There exist multiple eigenvalues of the transfer matrix. By definition, the minimum positive ${\lambda }_l$ is selected as the localization factor [39]. Since ${\epsilon }_i$ is the random variable in a unit cell, it varies from one unit cell to another. If finite numbers of unit cells $N$ are selected, the calculated localization factor is also random and determined by the value of $N$. When $N\to \mathrm{\infty }$, the localization factor converges to a certain value. Therefore, a large number of unit cells are often needed to ensure the convergence of the localization factor. The present calculation shows that when $N$ is greater than 5000, the localization factor converges, and the localization factor corresponding to each frequency remains a fixed value. With random parameter given for each cell sample, an approximate localization factor can be obtained, namely:

(45) ${{\lambda }_l}^{(N)}=\frac{1}{N}\sum _{i=1}^{N}ln\left|{\lambda }_{i,l}({\epsilon }_i)\right|$(45)

5. Results and discussion

Although the flexoelectric effect can exist in crystals with arbitrary symmetry, it is more attractive in crystals that exhibit centrosymmetry in the absence of piezoelectric properties. In other words, the force-electric coupling phenomenon can be solely attributed to the flexoelectric effect in centrosymmetric crystals. Therefore, in this section, two typical centrosymmetric dielectric materials, $\mathrm{B}\mathrm{a}\mathrm{T}\mathrm{i}{\mathrm{O}}_3$ (cubic phase) and $\mathrm{S}\mathrm{r}\mathrm{T}\mathrm{i}{\mathrm{O}}_3$, are selected for subsequent calculations in order to illustrate the influence of flexoelectricity on the structure of energy band. Specific material and geometrical features [20,40] are summarized in Table 1.

Table 1. Material and geometry properties of 1D PC.

	$\mathrm{B}\mathrm{a}\mathrm{T}\mathrm{i}{\mathrm{O}}_3$	$\mathrm{S}\mathrm{r}\mathrm{T}\mathrm{i}{\mathrm{O}}_3$
$Relativepermittivity(\mathrm{a})$	$4\times {10}^3$	$3\times {10}^2$
$\rho (\mathrm{k}\mathrm{g}/{\mathrm{m}}^3)$	$6.02\times {10}^3$	$5.12\times {10}^3$
$b(\mathrm{N}{\mathrm{m}}^4/{\mathrm{C}}^2)$	$6.77\times {10}^{-6}$	$4.14\times {10}^{-6}$
$h(\mathrm{N}\mathrm{m}/\mathrm{C})$	$-1.03\times {10}^3$	$-1.20\times {10}^3$
$c(\mathrm{N}/{\mathrm{m}}^2)$	$1.62\times {10}^{11}$	$3.5\times {10}^{11}$
$\eta (\mathrm{N}{\mathrm{m}}^2/{\mathrm{C}}^2)$	$1.26\times {10}^{11}$	$1.13\times {10}^{11}$

It is worth noting that the value of micro-inertia characteristic length $l$ is set to the lattice constant according to the literature [41]. The lattice constant of $\mathrm{B}\mathrm{a}\mathrm{T}\mathrm{i}{\mathrm{O}}_3$ in the cubic phase and $\mathrm{S}\mathrm{r}\mathrm{T}\mathrm{i}{\mathrm{O}}_3$ is 0.39 nm.

When considering the first sublayer thickness detuning, it is assumed that $d_1$ obeys a uniform distribution with mean ${\bar{d}}_1$ and detuning degree (coefficient of variation) $\delta$. The range for the values of $d_1$ can be shown as:

(46) $d_1\in \left[{\bar{d}}_1(1-\sqrt{3}\delta ),{\bar{d}}_1(1+\sqrt{3}\delta )\right]$(46)

Further, let $r$ represent a random variable that is subjected to a standard uniform distribution, i.e. $r\in (0,1)$. Then, $d_1$ can be compactly expressed as:

(47) $d_1={\bar{d}}_1\left[1+\sqrt{3}\delta (2r-1)\right]$(47)

To facilitate the calculation, the thickness of the two sublayers in the unit cell can be nondimensionalized by ${\bar{d}}_1$. With random detuning considered, let ${\zeta }_1$ and ${\zeta }_2$ represent essentially the dimensionless widths of sublayers A and B. ${\zeta }_2=1$ means that only the detuning of sublayer A is discussed, and the width of sublayer B remains fixed. It is obvious that the detuned PC will turn into a classical harmonic PC if $\delta =0$. In the following, the influence of flexoelectric effect and parameters detuning (such as thickness, flexoelectric coefficient) on band gap and localization factors are analyzed and discussed.

5.1 Harmonic PC with and without flexoelectric effect

First, to explore the influence of flexoelectric effect, Figure 2 shows the energy band structure of $\mathrm{B}\mathrm{a}\mathrm{T}\mathrm{i}{\mathrm{O}}_3-\mathrm{S}\mathrm{r}\mathrm{T}\mathrm{i}{\mathrm{O}}_3$ bilayer PC and the corresponding localization factor, with the thickness of both sublayers taken as 0.8 nm. Specifically, as shown in Figure 2, the red and black scatter points represent the energy band results with flexoelectricity and without flexoelectricity, respectively. The red line represents the energy band result with flexoelectricity and micro-inertia effect. The detailed formulas of the transfer matrix are summarized in the Appendix. It is seen from Figure 2 that the band gap starts and ends at 0 or π, which belongs to a typical Bragg-type band gap. In the presence of flexoelectricity, the band gap curve is no longer straight, representing a significant dispersion phenomenon of wave. Besides, it is interesting to note that, with flexoelectricity effect considered, the band gap becomes higher and broader. In other words, flexoelectricity increases the center frequency and width of the first band gap compared to the classical elastic results. By the way, note that the slope of the dispersion curve corresponds to the phonon group velocity, which affects the heat transport characteristics of nanoscale materials. Therefore, the present results show that flexoelectricity leads to larger phonon group velocity and hence larger thermal conductivity [42]. Besides, the micro-inertial effect has little effect on the energy band structure of our PC model that already considers flexoelectricity, so we follow the approach in the literature [20,21] and ignore the micro-inertial effect in the following discussions.

Figure 2. Energy band structure and variation of localization factor with frequency ($d_{BaTiO3}=d_{SrTiO3}=0.8$ nm).

As shown in Figure 2, it is clear that the frequency range of the band gap coincides with that of the nonzero localization factor. In fact, the localization factor is another effective parameter that can well characterize the energy band structure in 1D harmonic PC [39]. According to the definition of localization factor, when it is nonzero, the propagation of an elastic wave is attenuated and the attenuation coefficient is $e^{-\lambda }$. Instead, when the localization factor is zero, the elastic wave can pass smoothly through the multilayer medium, which corresponds to the passband in the energy band structure. Besides, in the presence of flexoelectric effect, the localization factor is remarkably larger than that in the purely elastic state, as shown in Figure 2, which means that flexoelectric effect can increase the attenuation rate of elastic waves.

In conclusion, flexoelectric effect influences significantly the band gaps of nanoscale PC, and it is necessary to consider the influence of flexoelectricity. Localization factor is an effective parameter to characterize the band structure of 1D harmonic PC.

5.2 Influence of thickness detuning with fixed thickness ratios

The primary aim of the present study is to study wave localization in a detuning PC using the localization factor when flexoelectric effect is considered, since flexoelectric effect cannot be ignored as discussed above. For this purpose, the material thickness ratio of 1D PC is reduced from 1 to 0.5 to study the influence of different detuning degrees $\delta$ (i.e. 0.1 and 0.2), as shown in Figure 3. For reference, the corresponding variation curve of localization factor of a harmonic PC (i.e. $\delta =0$) is also presented. At first glance, it is seen that the peak of localization factor decreases significantly as the detuning degree is increased. Besides, it is interesting to note that, as the detuning degree $\delta$ is increased, the band gap of the PC exhibits a widening trend, i.e. the detuning degree can widen the band gap. As for different thickness ratios, the influence law of detuning degree is similar. In the next section, the influence of thickness ratio will be discussed in detail.

Figure 3. Variation curve of localization factor with thickness detuning: (a) thickness ratio $d_1/d_2=1$; (b) thickness ratio $d_1/d_2=0.5$.

In addition to localization factor and energy band, we further analyze the relationship between bandwidth and detuning degree. As shown in Figure 4, it is found that the first and second bandwidths show a clear increasing trend with the increase of detuning. Specifically, the bandwidth of the first band gap is increased from 51.5 GHz to 136.5 GHz, and the bandwidth of the second band gap is increased from 119.5 GHz to 306.5 GHz. It is interesting to note that the trend exhibits a positive linear correlation, i.e. there is a clear regularity in the increase of forbidden bandwidth with the change of detuning degree, which is important for designing instruments with specific bandwidths. According to the above rules, one can change the degree of detuning to modulate a PC with forbidden band to meet specific needs. By the way, the results of Figure 4 also show that the bandwidth of the second band gap is significantly larger than that of the first band gap.

Figure 4. Variation of bandwidth with detuning degree for $\mathrm{B}\mathrm{a}\mathrm{T}\mathrm{i}{\mathrm{O}}_3-\mathrm{S}\mathrm{r}\mathrm{T}\mathrm{i}{\mathrm{O}}_3$ PC with $d_1/d_2=1$.

Physically, detuning implies a disorder in the structural parameters (dimensions or materials, etc.) of the original perfectly periodic PC. The disorder leads to partial reflections of the wave at each bay, which can occur even at the passband frequencies of the corresponding periodic system [43]. Such multiple reflections lead, on the one hand, to a failure of the wave to travel smoothly downstream, thus showing an increase in forbidden band bandwidth with increasing detuning in the bandwidth diagram. On the other hand, the multiple reflections lead to a phenomenon where the energy is actually concentrated in a certain part of the structure, which is also commonly referred to as localization [44]. In other words, the increase of the forbidden band width and the localization of energy are two macroscopic phenomena caused by detuning, and the mechanism behind them is actually the same.

5.3 Influence of thickness ratio with fixed thickness detuning

In this section, the thickness ratio of 1D PC is varied from 0.5 via 1 to 2 when the thickness detuning is either fixed at 0 or 0.2. As shown in Figure 5, when the thickness detuning is fixed, decreasing the thickness ratio enables the beginning and end of the forbidden band to shift toward lower frequencies. This means that, with the thickness detuning fixed, 1D PC with relatively low-frequency forbidden band can be obtained by decreasing its thickness ratio. Besides, it is interesting to note that as the thickness ratio is decreased, the localization factor in the forbidden band zone tends to increase, and the localization phenomenon becomes more pronounced.

Figure 5. Variation curve of localization factor with thickness ratio: (a) thickness detuning degree $\delta =0$; (b) thickness detuning degree $\delta =0.2$.

5.4 Influence of flexoelectric coefficient detuning with fixed thickness ratio

The influence of flexoelectric efficient detuning is discussed when fixing the thickness ratio at 0.5, 1, and 2. Specifically, three detuning degrees (i.e. 0, 0.002, and 0.004) are selected, with the $\delta =0$ case of harmonic PC taken as reference. As shown in Figure 6, the localization factor decreases with increasing detuning of the flexoelectric coefficient. Upon comparing Figure 6 with Figure 3, it is seen that the effect of flexoelectric coefficient detuning on localization factor is similar to that of thickness-induced detuning within the entire frequency range considered. From a specific value of detuning, a very slight (about 0.001) fluctuation of the flexoelectric coefficient can lead to the localization phenomenon, implying that the sensitivity of flexoelectric effect to localization is much higher than that of the thickness parameter. This once again demonstrates that flexoelectric effect is important for dynamic analysis and cannot be ignored in nanoscale structures.

Figure 6. Influence of flexoelectric coefficient detuning on localization factor versus frequency curve with the thickness ratio fixed at (a) $d_1/d_2=0.5$ and (b) $d_1/d_2=1$.

5.5 Influence of thickness ratio with fixed flexoelectric coefficient detuning

Finally, the influence of thickness ratio on localization is studied by fixing the flexoelectric coefficient detuning at either 0 or 0.004. Three thickness ratios (i.e. 0.5, 1, and 2) are selected, and the results are displayed in Figure 7. Generally, when the detuning degree is fixed, the localization phenomenon due to flexoelectric efficient detuning is similar to that due to thickness detuning shown in Figure 5. In particular, decreasing the thickness ratio shifts the beginning and end of the forbidden band to lower frequencies. In other words, for 1D PC with given flexoelectric detuning, one can obtain PC with relatively low-frequency forbidden band by decreasing the thickness ratio. Besides, it is noted that as the thickness ratio is decreased, the localization factor in the forbidden band is also increased and the localization phenomenon will be more obvious.

Figure 7. Effect of thickness ratio on localization factor versus frequency curve with flexoelectric efficient detuning degree fixed at (a) $\delta =0$ and (b) $\delta =0.004$.

6. Conclusions

The variation of energy band structure of 1D harmonic and detuned PC has been analytically investigated, with flexoelectric effect accounted for. Generally, the calculated results show that detuning under flexoelectric effect affects significantly the band structure of the PC. First, the localization factor indeed reflects the energy attenuation rate and can effectively represent the band structure, including both the pass band and the forbidden band. Second, increasingly either the thickness detuning or flexoelectric coefficient detuning leads to linearly widened band gap, and the second band gap width increases faster than the first band gap. In addition, the sensitivity of flexoelectric coefficient detuning is much higher than that of thickness detuning. The present results provide a theoretical basis for designing high-frequency nanoscaled devices with the function of filtering based on band gap effect in PC.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (11502110, 11972185, and 12002156), the Open Fund of the State Key Laboratory for Strength and Vibration of Mechanical Structures, Xi’an Jiaotong University (No. SV2018-KF-01, No. SV2018-KF-22), and the Open Project of Space Structure and Mechanism Technology Laboratory of China Aerospace Science and Technology Group Co., Ltd.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the National Natural Science Foundation of China [11972185,11502110].

Appendix

1. Derivation of Equations (23) and (24)

The simplified inhomogeneous equations are:

$c\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}+h\frac{{\mathrm{\partial }}^{2}P}{\mathrm{\partial }{x}^2}=\rho \frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{t}^2}-\rho l^2\frac{{\mathrm{\partial }}^{4}u}{\mathrm{\partial }{x}^2\mathrm{\partial }{t}^2}$

$h\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{x}^2}+b\frac{{\mathrm{\partial }}^{2}P}{\mathrm{\partial }{x}^2}-\eta P=0$

The solution of the above equations can be expressed as:

$u(x,t)=\bar{U}(x)exp(-i\omega t)$

$P(x,t)=\bar{P}(x)exp(-i\omega t)$

where $i=\sqrt{-1}$. Substituting the assumed solution into the inhomogeneous equations leads to:

$\left(c+\rho l^2{\omega }^2\right)\frac{{\mathrm{\partial }}^2\bar{U}(x)}{\mathrm{\partial }{x}^2}+h\frac{{\mathrm{\partial }}^2\bar{P}(x)}{\mathrm{\partial }{x}^2}+\rho {\omega }^2\bar{U}(x)=0$

$h\frac{{\mathrm{\partial }}^2\bar{U}(x)}{\mathrm{\partial }{x}^2}+b\frac{{\mathrm{\partial }}^2\bar{P}(x)}{\mathrm{\partial }{x}^2}-\eta \bar{P}(x)=0$

To facilitate the calculation, the dimensionless coordinate, ${\zeta }_i=\frac{x_j}{{\bar{d}}_1}$, is introduced. where ${\bar{d}}_1$ represents the average thickness of the first sublayer. Therefore, the two equations above can be further simplified as:

$\left(c+\rho l^2{\omega }^2\right)\frac{{\mathrm{\partial }}^2\bar{U}(\zeta )}{\mathrm{\partial }{\zeta }^2}+h\frac{{\mathrm{\partial }}^2\bar{P}(\zeta )}{\mathrm{\partial }{\zeta }^2}+{{\bar{d}}_1}^2\rho {\omega }^2\bar{U}(\zeta )=0$

$h\frac{{\mathrm{\partial }}^2\bar{U}(\zeta )}{\mathrm{\partial }{\zeta }^2}+b\frac{{\mathrm{\partial }}^2\bar{P}(\zeta )}{\mathrm{\partial }{\zeta }^2}-{{\bar{d}}_1}^2\eta \bar{P}(\zeta )=0$

Appendix

1. Calculation of transfer matrix with flexoelectric effect neglected

With the higher order terms and force-electric coupling caused by flexoelectric effect and micro-inertia effect neglected, the governing equation is simplified as:

$\tilde{c}\frac{d^{2}W}{d^{2}x}+\rho {\omega }^{2}W=0$

The displacement can be expressed as:

$W(x)=e^{i\gamma x}A+e^{-i\gamma x}B$

where $\gamma =\omega \sqrt{\frac{\rho }{\tilde{c}}}$, while the stress is given by:

$\sigma =\tilde{c}\frac{dW}{dz}=i\gamma \tilde{c}{e}^{i\gamma x}A-i\gamma \tilde{c}{e}^{-i\lambda x}B$

In the classical elastic case, the state vector is $v=M_j(x)\{W,\sigma {\}}^T$, where $M_j(x)=\left(\begin{array}{cc}{e}^{i\lambda x}& e^{-i\lambda x}\\ i\lambda \tilde{c}{e}^{i\lambda x}& -i\lambda \tilde{c}{e}^{-i\lambda x}\end{array}\right)$, and $j=1$ and $2$ represent the sublayers A and B of 1D PC. The transfer matrix is then obtained as $T_n=M_2(d)[M_2(d_1){]}^{-1}{M}_1(d_1)[M_1(0){]}^{-1}$.

Appendix

1. Expression of state vector with flexoelectric effect accounted for

The state vectors are: $v=M_{\mathrm{j}}(\zeta )\{\bar{U},\bar{P},{\sigma }^{phys},V^{phys}{\}}^T$

$\begin{array}{rl}& M_j(1,1)=e^{l_{1j}\zeta }{M}_j(1,2)=e^{-l_{1j}\zeta }{M}_j(1,3)=e^{i{l}_{2j}\zeta }{M}_j(1,4)=e^{-i{l}_{2j}\zeta }\\ & M_j(2,1)=\frac{1}{-h_j}(c_j+{\rho }_j{l}_j^2{\omega }^2+\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{1j}^2})e^{l_{1j}\zeta }{M}_j(2,2)=\frac{1}{-h_j}(c_j+{\rho }_j{l}_j^2{\omega }^2+\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{1j}^2})e^{-l_{1j}\zeta }\\ & M_j(2,3)=\frac{1}{-h_j}(c_j-\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{2j}^2})e^{i{l}_{2j}\zeta }{M}_j(2,4)=\frac{1}{-h_j}(c_j-\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{2j}^2})e^{-i{l}_{2j}\zeta }\\ & M_j(3,1)=-(\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{1j}})e^{l_{1j}\zeta }{M}_j(3,2)=(\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{1j}})e^{-l_{1j}\zeta }\\ & M_j(3,3)=i(\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{2j}})e^{i{l}_{2j}\zeta }{M}_j(3,4)=-i(\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{2j}})e^{-i{l}_{2j}\zeta }\\ & M_j(4,1)=l_{1j}{h}_j[1-\frac{b_j}{{h_j}^2}(c_j+{\rho }_j{l}_j^2{\omega }^2+\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{1j}^2})]e^{l_{1j}\zeta }\\ & M_j(4,2)=-l_{1j}{h}_j[1-\frac{b_j}{{h_j}^2}(c_j+{\rho }_j{l}_j^2{\omega }^2+\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{1j}^2})]e^{-l_{1j}\zeta }\\ & M_j(4,3)=i{l}_{2j}{h}_j[1-\frac{b_j}{{h_j}^2}(c_j-\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{2j}^2})]e^{i{l}_{2j}\zeta }\\ & M_j(4,4)=-i{l}_{2j}{h}_j[1-\frac{b_j}{{h_j}^2}(c_j-\frac{{{\bar{d}}_1}^2{\rho }_j{\omega }^2}{l_{2j}^2})]e^{-i{l}_{2j}\zeta }\end{array}$

where $j=1$ and $2$ represent the sublayers A and B of 1D PC, respectively.