Twist piezoelectricity: giant electromechanical coupling in magic-angle twisted bilayer LiNbO3

Twisted a pair of stacked two-dimensional materials exhibit many exotic electronic and photonic properties, leading to the emergence of flat-band superconductivity, moiré engineering and topological polaritons. These remarkable discoveries make twistronics the focus point of tremendous interest, but mostly limited to the concept of electrons, phonons or photons. Here, we present twist piezoelectricity as a fascinating paradigm to modulate polarization and electromechanical coupling by twisting precisely the stacked lithium niobate slabs due to the interlayer coupling effect. Particularly, the inversed and twisted bilayer lithium niobate is constructed to overcome the intrinsic mutual limitation of single crystals and giant effective electromechanical coupling coefficient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${k}_{t}^{2}$$\end{document}kt2 is unveiled at magic angle of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$111^\circ$$\end{document}111∘, reaching 85.5%. Theoretical analysis based on mutual energy integrals shows well agreements with numerical and experimental results. Our work opens new venues to flexibly control multi-physics with magic angle, stimulating progress in wideband acoustic-electric, and acoustic-optic components, which has great potential in wireless communication, timing, sensing, and hybrid integrated photonics.

In a piezoelectric crystal, the acoustic fields equations: where  is stress tensor,  ⃗ ⃗ is the particle displacement,  is strain tensor, and the electromagnetic fields equations: where  ⃗⃗ is electric fields,  ⃗⃗ is magnetic induction intensity fields,  ⃗⃗⃗ is magnetic strength fields and  ⃗⃗ is electric displacement, are coupled through the piezoelectric constitutive relations, which can be given in strain-charge form: ⃗⃗ =  T ⋅  ⃗⃗ + :  (S5) where  T is permittivity at zero or constant stress,  is piezoelectric strain constants,   is compliance constants, as well as in stress-charge form: ⃗⃗ =  S ⋅  ⃗⃗ + :  (S7) where  S is permittivity at zero or constant strain,  is piezoelectric stress constants,   is elastic constants at zero or constant electric field.Combine the equations S1-S8, we obtain the piezoelectrically coupled wave equations in piezoelectric crystals as given as equation (1) and equation (2) in main text of this letter.The matrix form of equation (1) and equation ( 2) can be expressed as: and where As can be seen that there are three solutions to equation (S16) while one corresponding to pure longitudinal acoustic waves and the other two correspond to the ordinary light and the unusual light in LN crystal.Thus, the acoustic modes were called quasi-longitudinal wave, which coupled with shear polarized quasi-electromagnetic waves.
The dispersion of another two SH polarized waves were solved as: Since the forms of piezoelectric coupling apply an additional coefficient calculated by  and  , which is mathematically different from coupling form of quasi-waves, these kinds of piezoelectric acoustic waves are called as piezoelectrically stiffened waves.The longitudinal polarized coupled electric fields indicates that the curl of the electric fields maintaining zero even in time harmonic vibration and thus there is no time harmonic vibrating magnetic fields.Furthermore, the electric displacement fields should be zero indicated by equation (S4).Therefore, the coupled longitudinal electric fields can be calculated as: The irrotational electric fields allows to replace electric fields vector with a scalar potential Φ as  ⃗ = −∇Φ, which is known as the quasi-static approximation and used primarily in FEM simulations and applications.
Another important characteristic of SH stiffened waves is its distribution of stress, derived by equation (S1) and equation (S2), can be expressed as: with  1 equals to zero as well as  2 ,  3 and  4 .Thus, the corresponding stress fields exhibit the same polarization as its particle velocity/displacement polarization, with the phase on the contrary because of the partial derivatives with respect to space involves an additional factor − to its amplitude.

B. Definition of "∇" operator
According to the variable type of the object of the operator "∇", there are three types of calculations need to be defined: where the operator acted on a symmetric second order tensor (like  or ) and the result should be a vector (like  ⃗ ⃗ or where the operator acted on a vector and the result should also be a vector, where the operator acted on a vector and the result should be a symmetric second-order tensor.

A. Partial wave methods
The dispersion characteristics and vibrating fields distributions of acoustical resonant eigenmodes in ITBLN, acoustic fields and related boundary conditions were solved by superposition of partial waves.All possible plane wave solutions which can be also called waveguide modes should be firstly solved for sake of superimposition with amplitudes and then be determined by related boundary conditions, which take form as: where  ̂ is the unit vector normal to the boundary.In this letter, the used partial waves were the two SH polarized waves solved in Supplementary I which can be expressed as: where T is refraction coefficients and Γ is reflection coefficients, solved by the boundary conditions as equation (S23) and equation (S24).To be noted that the superscript of coefficients indicates the medium from which the incident wave propagates, more specifically says, 1 indicates the upper medium and 2 indicates the lower medium.
Additionally, the boundary conditions on the upper and lower sides of cavity can be expressed as a total reflected equation: where  2×2 is second order identity matrix.Simultaneous solving equation (S17), equation (S18) and equation (S19), the determinant of dispersion characteristic and eigenfrequency can be obtained:

B. Symmetric distribution of stress
Supplementary Figure 1.Schematics of inversion and rotation of ITBLN, the front view of interface between upper and lower layers shows the defined coordinate angular bisector systems which used for derivation of symmetry of distributions of stress fields.
Using the angular bisector defined coordinate system (see Supplementary Figure 1), the component of particle where  , is acoustic impedance defined as / , , which can be solved by equation (S17).If the thickness of the two LN layers is the same, marked as ℎ, two series of symmetry can be involved as  , = ± , as well as  , = ± , .It is evident that the two series exhibite the opposite symmetric of stress and particle velocity distributions as indicated by equations (S29-S36).
Supplementary Figure 2. Electrical responses of ITBLN at different thickness ratios while the angle were set as a, 116° and b, 297°.
To be noted that the equations (S29-S36) are valid only on the boundary between the two LN layers if the origin of the X-axis is set at this boundary.As for the distribution characteristics along X direction in the whole cavity, the rotation of polarization would be expected to be remarkable as the result of different phase velocities/wavenumbers of  and  type SH partial waves for superposition of total fields.
Once the thickness of two piezoelectric layers are unequal, the symmetry of fields distribution would be broken.Two extra parasitic modes can be observed, as shown in Supplementary Figure 2. Despite the coupling coefficients of these modes are quite small, the parasitic modes would degrade the performance of acoustic devices in potential applications.

C. Calculation of refraction 𝐓 and reflection 𝚪 coefficients
In the first case, we consider an -type polarized SH wave incident from one LN layer into the other LN layer, as seen in Fig. S3(a), at the interface between two layers, there would be two reflected waves(-type polarized,  ⃗  and -type polarized,  ⃗  ) and two refracted waves (-type polarized,  ⃗ T and -type polarized,  ⃗ T ).As a ratio of amplitudes of these waves, the refraction and reflection coefficients (T and Γ) can be involved in the expression of waves at the interface: To simplify the derivation for a more intuitive and concise solutions with physical pictures, the piezoelectricity was neglected here.Thus, the stress of all plane waves presupposed above can be calculated by the combination of equation ( S14) and the following equation: Supplementary Figure 3. Schematic of reflection and refraction of a, -polarized incident waves and c, -polarized incident waves with the b, coordinated system used in derivation of coefficients calculations.
Before the calculation of refraction and reflection coefficients, the influence of piezoelectricity term −    in equation (S42) on the acoustic impedance needs to be analyzed.Firstly, since either  and  exhibits shear-horizontal polarization (  and   ) and merely couples with longitudinal electric fields (  ), substituting Eq(S11), (S12) and (S22) into (S42) would lead to nonzero value of  5 and  6 as well as the other components remain zero, given as: Here we marked the two solutions to Eq(451) or ( 452) as  and , which corresponds to the  and  modes solved by Eq(448).Obviously, both  and  must satisfy the equations above.
Thirdly, with help of the polarization equations derived above, the stress and displacement of  and  modes can be related.Considering a wave propagating alone +x direction with  or  polarization type, its displacement field can be written as: where the harmonic term exp( − ) is omitted.Substituting  into Eq(S53) gives: ⃗ ⃗ =   ( ̂+ ) (S54) Meanwhile, according to the projection of stress along x axis: ⋅  ̂=  1  ̂+  6  ̂+  5 ̂(S55) Combination with Eq(S47) and (S48), the projection of stress field can be written as: where relation Eq(S52) was used during the derivation.To be noted that the complex form of displacement and stress field involves an imaginary part in expressions, which defines the phase.On the other hand, the definition of acoustic impedance   gives: where  equals either to the two solutions ( or ).Therefore, the projection stress of  or  mode can be related with its displacement field: ⋅  ̂= −   ⃗ ⃗ = −   ⃗ ⃗ (S59) Furthermore, if the wave propagates along -x direction, its harmonic term should be exp( + ).Under this [ ] is a 6 × 6 sized matrix, representing the elastic constants of LN and [] is a 6 × 3 sized matrix, representing the piezoelectric constants of LN and [ϵ S ] is a 3 × 3 sized matrix, representing the permittivity constants −  ⋅    −  ⋅ {    +  ⋅    +  ⋅ in lower LN layer.The subscript  and  indicates different polarization types.The reflection and refraction coefficients at the boundary between the two piezoelectric mediums require that: −   = (̂cos   +  ̂sin   ) −   (S37)  ⃗  =  ̂1Γ   +   = Γ  (̂cos   +  ̂sin   ) +   (S38)  ⃗  = ̂1Γ   +   = Γ  (̂cos   −  ̂sin   ) +   (S39)  ⃗ T =  ̂2T   −   = T  (̂cos   −  ̂sin   ) −   (S40)  ⃗ T = ̂2T   −   = T  (̂cos   +  ̂sin   ) −   (S41)