Piezoresponse in Ferroelectric Materials under Uniform Electric Field of Electrodes

Abstract

The analytical solution for the displacements of an anisotropic piezoelectric material in the uniform electric field is presented for practical use in the “global excitation mode” of piezoresponse force microscopy. The solution is given in the Wolfram Mathematica interactive program code, allowing the derivation of the expression of the piezoresponse both in cases of the anisotropic and isotropic elastic properties. The piezoresponse’s angular dependencies are analyzed using model lithium niobate and barium titanate single crystals as examples. The validity of the isotropic approximation is verified in comparison to the fully anisotropic solution. The approach developed in the paper is important for the quantitative measurements of the piezoelectric response in nanomaterials as well as for the development of novel piezoelectric materials for the sensors/actuators applications.

1. Introduction

Piezoelectric materials are an important class of materials with applications as sensors, actuators, resonators [1,2,3], electric energy harvesters [4,5], in various microelectronic devices [6,7,8], and in the piezoelectric catalysis for wastewater treatment [9]. Evaluating the effective piezoelectric coefficients via macroscopic response and predicting the response based on the known piezoelectric coefficients are vital in designing efficient electromechanical devices based on new materials. Piezoresponse force microscopy (PFM) is a local technique allowing the measurement of the surface piezoresponse with a unique high spatial resolution of 10–20 nm and vertical sensitivity around 100 fm [10,11]. PFM allows measuring the piezoresponse in the micro and nano-objects where macroscopic measurement techniques fail and PFM becomes the only way to quantify the material’s piezoelectric coefficients [12,13,14,15,16,17]. Despite the rapid development of the PFM technique within the last 30 years, many issues still remain, such as the contribution from the electrostatic force [18,19,20] and other parasitic effects [11,21,22,23], as well as the difficulties of the PFM response quantification, i.e., understanding the relationship between the measured surface displacement and components of the piezoelectric tensor [17,24,25,26,27,28].

The most straightforward PFM implementation is in the so-called “global excitation mode”, where an electric field is applied through micro- or nano-sized electrodes, while signal registration is performed utilizing scanning probe microscopy (SPM) probe rastering across the surface [16,29,30,31,32,33,34,35]. Global excitation can be realized via an electric field applied to the conductive tip across the tip-electrode interface and simply by the excitation of the electrode with probe acting as a mechanical sensor. In both cases, the solution of the problem of electromechanical response quantification in the applied electric field is close to the conventional laser interferometry problem of the response under the action of the uniform or close-to-the uniform electric field. Contrary to the highly non-uniform electric field from the SPM probe, where the solution is usually limited to the case of the uniform elastic properties or needs numerical calculations [24,26,27,36,37], the piezoelectric response in the uniform electric field allows an accurate solution. Lefki and Dormans derived the solution [38] for the case of the tetragonal (001) oriented ferroelectric materials bonded to the rigid substrate, as it is imperative to the interpretation of the Doppler vibrometer and laser interferometry data [39,40,41]. More complicated cases considering the flexible substrate contribution, top electrode lateral size, and multilayered structure were studied using sophisticated analytical models and finite element simulations [41,42,43,44,45,46,47].

In this paper, we provide the complete solution of this problem for the arbitrarily oriented ferroelectric sample covered by the top electrode and fixed on the rigid substrate for the cases of the transversally isotropic and fully anisotropic elastic properties. The solutions reveal the contribution of the elastic anisotropy to the piezoresponse, as well they are of interest for the determination of the piezoelectric coefficients from the displacement of the piezoelectric material in the uniform electric field. The results are important for understanding the “global excitation mode” of PFM and conventional PFM response in the case of the measurements in the piezoelectric thin films with the thickness much lower than the tip radius (e.g., for 2D piezoelectric materials) [48]. The solutions are represented in the Wolfram Mathematica interactive code, allowing the visualization and study of the equations in the crystals with arbitrary elastic/piezoelectric matrices and orientation.

2. Theoretical Framework

The direct piezoelectric effect is the appearance of the electric charge at the surface of the crystal with a specific symmetry class under the mechanical force’s action. The polarization vector is linearly connected with the mechanical stress components ${\sigma }_{kl} \left(k,l=\mathrm{1,2,3}\right)$ by the equation:

$$P_i=d_{ikl}{\sigma }_{kl},$$

(1)

where $d_{ikl}$ are the piezoelectric coefficients, representing the 3rd rank tensor.

PFM usually exploits the converse piezoelectric effect representing a change of the piezoelectric crystal dimensions under the external electric field’s action. The deformation tensor ${ϵ}_{ij}$ is linearly dependent on the electric field vector $E_i$ in the mechanically free crystal (${\sigma }_{ij}=0$):

$${ϵ}_{ij}=d_{kij}{E}_k.$$

(2)

In both Equations (1) and (2), the same piezoelectric coefficient appears that can be justified by the thermodynamical analysis. In the general case, the third rank tensor has 27 independent components. ${ϵ}_{ij}={ϵ}_{ji}$, $d_{kij}$ piezoelectric tensor is symmetrical by the two last indexes, which leads to the reduction of the independent component number to 18.

The problem of the crystal deformation under the action of the external electric field can be written in the isothermal state ($\theta =0$) as:

$${\sigma }_{ij}=c_{ijkl}^{E,T}{ϵ}_{kl}-e_{ijm}^T{E}_m,$$

$$D_m={\epsilon }_{mk}^{ϵ,T}{E}_k+e_{mij}^T{ϵ}_{ij},$$

(3)

where D is an electric induction, T is a temperature.

To shorten the notations, upper indexes will be further omitted. The equation of the motion for the piezoelectric media without an accounting of the mass forces can be written as:

$$\frac{\partial {\sigma }_{ij}}{\partial x_j}=\rho \frac{{\partial }^2{u}_i}{\partial t^2}.$$

(4)

Neglecting the magnetic effects and in the quasi-static approximation for the electric field:

$$\frac{\partial D_i}{\partial x_i}=0.$$

(5)

Substituting the expression for the electric field $E_i=-\frac{\partial \phi }{\partial x_i}$ through the scalar potential to the equation of the state (3), the following can be obtained:

$${\sigma }_{ij}=c_{ijkl}\frac{\partial u_k}{\partial x_l}+e_{kij}\frac{\partial \phi }{\partial x_k}$$

$$D_i=e_{ikl}\frac{\partial u_k}{\partial x_l}-{\epsilon }_{ik}\frac{\partial \phi }{\partial x_k}$$

(6)

The system of the linear equation of the piezoelectric media electro-elasticity can be obtained by the substitution equations of the state (6) to the equations of the motion (4) and electrostatics (5):

$$c_{ijkl}\frac{{\partial }^2{u}_k}{\partial x_l\partial x_j}+e_{kij}\frac{{\partial }^2\phi }{\partial x_k\partial x_j}=\rho \frac{{\partial }^2{u}_i}{\partial t^2},$$

(7)

$$e_{ikl}\frac{{\partial }^2{u}_k}{\partial x_l\partial x_i}-{\epsilon }_{ik}\frac{\partial \phi }{\partial x_k\partial x_i}=0.$$

(8)

Let us consider an infinite layer of the piezoelectric material with arbitrary anisotropy of the elastic, dielectric, and piezoelectric properties covered by the electrode from both sides and bonded on the rigid grounded substrate under constant voltage and in the condition of the free surface (Figure 1a). The task is one-dimensional, which means functions changing only along the axis z. In the equilibrium state, it transforms (7) to the system of the uniform differential equation, which can be expressed in the Voigt notation for the tensors as:

$$\left\{\begin{array}{c}\begin{array}{c}{e}_{35}{\phi }^{″}+c_{55}{u}_x^{″}+c_{45}{u}_y^{″}+c_{35}{u}_z^{″}=0\\ e_{34}{\phi }^{″}+c_{45}{u}_x^{″}+c_{44}{u}_y^{″}+c_{34}{u}_z^{″}=0\end{array}\\ \begin{array}{c}{e}_{33}{\phi }^{″}+c_{35}{u}_x^{″}+c_{34}{u}_y^{″}+c_{33}{u}_z^{″}=0\\ e_{35}{u}_x^{″}+e_{34}{u}_y^{″}+e_{33}{u}_z^{″}-{\epsilon }_{33}{\phi }^{″}=0\end{array}\end{array}\right.,$$

(9)

where e_ij are piezoelectric constants.

The boundary conditions for the grounded stationary substrate:

$$\left\{\begin{array}{c}{\left.\phi \right|}_{z=0}=0\\ {\left.u_x\right|}_{z=0}=0\\ {\left.u_y\right|}_{z=0}=0\\ {\left.u_z\right|}_{z=0}=0\end{array}\right..$$

(10)

The boundary conditions for the free electrode surface under the voltage $U$ without applied force can be derived using the condition ${\sigma }_{ik}{n}_k=0$ [49] by the substitution ${\sigma }_{ik}$ from (6):

$$\left\{\begin{array}{c}{\left.\phi \right|}_{z=h}=U\\ {\left.\left(e_{35}{\phi }^{\prime }+c_{55}{u}_x^{\prime }+c_{45}{u}_y^{\prime }+c_{35}{u}_z^{\prime }\right)\right|}_{z=h}=0\\ {\left.\left(e_{34}{\phi }^{\prime }+c_{45}{u}_x^{\prime }+c_{44}{u}_y^{\prime }+c_{34}{u}_z^{\prime }\right)\right|}_{z=h}=0\\ {\left.\left(e_{33}{\phi }^{\prime }+c_{35}{u}_x^{\prime }+c_{34}{u}_y^{\prime }+c_{33}{u}_z^{\prime }\right)\right|}_{z=h}=0\end{array}\right..$$

(11)

The solution of the boundary problem (9)–(11) can be derived as:

$$\left\{\begin{array}{c}\phi \left[z\right]=\frac{Uz}{h}\\ u_x\left[z\right]=-\frac{c_{35}{c}_{44}{e}_{33}-c_{34}{c}_{45}{e}_{33}-c_{34}{c}_{35}{e}_{34}+c_{33}{c}_{45}{e}_{34}+c_{34}^2{e}_{35}-c_{33}{c}_{44}{e}_{35}}{\left(c_{35}^2{c}_{44}-2c_{34}{c}_{35}{c}_{45}+c_{33}{c}_{45}^2+c_{34}^2{c}_{55}-c_{33}{c}_{44}{c}_{55}\right)}·\frac{Uz}{h}\\ u_y\left[z\right]=-\frac{-c_{35}{c}_{45}{e}_{33}+c_{34}{c}_{55}{e}_{33}+c_{35}^2{e}_{34}-c_{33}{c}_{55}{e}_{34}-c_{34}{c}_{35}{e}_{35}+c_{33}{c}_{45}{e}_{35}}{\left(c_{35}^2{c}_{44}-2c_{34}{c}_{35}{c}_{45}+c_{33}{c}_{45}^2+c_{34}^2{c}_{55}-c_{33}{c}_{44}{c}_{55}\right)}·\frac{Uz}{h}\\ u_z\left[z\right]=-\frac{c_{45}^2{e}_{33}-c_{44}{c}_{55}{e}_{33}-c_{35}{c}_{45}{e}_{34}+c_{34}{c}_{55}{e}_{34}+c_{35}{c}_{44}{e}_{35}-c_{34}{c}_{45}{e}_{35}}{\left(c_{35}^2{c}_{44}-2c_{34}{c}_{35}{c}_{45}+c_{33}{c}_{45}^2+c_{34}^2{c}_{55}-c_{33}{c}_{44}{c}_{55}\right)}·\frac{Uz}{h}\end{array}\right.$$

(12)

The complete solution is represented as a Wolfram Mathematica interactive program code (Appendix A) [50], which allows extracting a solution for any given orientation of the piezoelectric plate. In order to study the solution of the boundary problem (9)–(11) illustrating the behavior of the piezoelectric material in the uniform electric field, the solution was simplified for the case of the barium titanate and lithium niobate single crystals, which are model ferroelectric materials with the tabulated piezoelectric properties. These crystals possess different symmetries at room temperature: tetragonal and rhombohedral, respectively.

3. Results and Discussion

3.1. Barium Titanate: Material Tensors

Piezoelectric constant tensor $e_{ij}^0$ and piezoelectric coefficient tensor $d_{ij}^0$ can be written for the tetragonal barium titanate (4 mm point group) in the laboratory system of the coordinates as following [51]:

$$e_{ij}^0=\left(\begin{array}{cc}\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ e_{31}^0& e_{31}^0& e_{33}^0\end{array}& \begin{array}{ccc}0& e_{15}^0& 0\\ e_{15}^0& 0& 0\\ 0& 0& 0\end{array}\end{array}\right)$$

(13)

$$d_{ij}^0=\left(\begin{array}{cc}\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ d_{31}^0& d_{31}^0& d_{33}^0\end{array}& \begin{array}{ccc}0& d_{15}^0& 0\\ d_{15}^0& 0& 0\\ 0& 0& 0\end{array}\end{array}\right)$$

(14)

To write analytical equations for the displacement $\overrightarrow{u}$ through the piezoelectric coefficients, Equation (15) will be used for the relation between the tensors $e_{ijk}$ and $d_{ijk}$, which can lead to the following relationships in the laboratory system of the coordinates:

$$e_{mij}^0=d_{mkl}^0{c}_{klij}^0.$$

(15)

The anisotropic elastic properties for the barium titanate were taken from [52], while isotropic elastic properties: Young modulus and Poisson ratio were calculated by the Voigt averaging of the elastic tensor [53]. Piezoelectric and elastic properties used for the calculations are presented in Appendix B.

3.2. Barium Titanate: Isotropic Elastic Properties

In the approximation of the isotropic elastic properties, an elastic tensor can be reduced to the:

$$c_{ij}^0=c_{iso}=\left(\begin{array}{cc}\begin{array}{ccc}\lambda +2\mu & \lambda & \lambda \\ \lambda & \lambda +2\mu & \lambda \\ \lambda & \lambda & \lambda +2\mu \end{array}& \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{ccc}\mu & 0& 0\\ 0& \mu & 0\\ 0& 0& \mu \end{array}\end{array}\right)$$

(16)

where $\lambda$, $\mu$-Lame elastic modulus, which can be determined through Young’s modulus, $Y$, and Poisson coefficient, $\nu ,$ as follows:

$$\begin{array}{l}\lambda =\frac{\nu Y}{\left(1+\nu \right)\left(1-2\nu \right)}\\ \mu =\frac{Y}{2\left(1+\nu \right)}\end{array}$$

(17)

For the description of the sample rotation in relation to the laboratory coordinate systems, Euler angles $\left\{\varphi ,\theta ,\psi \right\}$ are used, which correspond to the subsequent counterclockwise rotation around the axis $Z_1{X}_2{Z}_3$ (Figure 1b). The rotational matrix used in the calculation can be found in Appendix C. The angle $\varphi$ is responsible for the rotation of the crystal in relation to the crystallographic axis $Z_c$, angle $\theta$-for the tilt of the axis $Z_c$ in relation to the axis $Z$, angle $\psi$-for the rotation of the plate with regard to the axis $Z$.

Equation (12) for the $\overrightarrow{u}$ displacement vector components can be written for the plate of the barium titanate with the thickness, $h$, infinite in the $XY$ plane, and rotated at the random angles $\left\{\varphi ,\theta ,\psi \right\}$ in the approximation of the elastic properties’ isotropy (16) as follows:

$$\left\{\begin{array}{c}{u}_x\left[z,\varphi ,\theta ,\psi \right]=\frac{\mathrm{Sin}\left[\theta \right]\mathrm{Sin}\left[\psi \right]2\left(1+\nu \right)\left(e_{31}^0-e_{33}^0+\left(2e_{15}^0+e_{31}^0-e_{33}^0\right)\mathrm{Cos}\left[2\theta \right]\right)}{Y}·\frac{Uz}{h}\\ u_y\left[z,\varphi ,\theta ,\psi \right]=\frac{\mathrm{Sin}\left[\theta \right]Cos\left[\psi \right]2\left(1+\nu \right)\left(e_{31}^0-e_{33}^0+\left(2e_{15}^0+e_{31}^0-e_{33}^0\right)\mathrm{Cos}\left[2\theta \right]\right)}{Y}·\frac{Uz}{h}\\ u_z\left[z,\varphi ,\theta ,\psi \right]=\frac{\mathrm{Cos}\left[\theta \right]\left(1+\nu \right)\left(2\nu -1\right)\left(2e_{15}^0+e_{31}^0+e_{33}^0-\left(2e_{15}^0+e_{31}^0-e_{33}^0\right)\mathrm{Cos} \left[2\theta \right]\right)}{2Y\left(1-\nu \right)}·\frac{Uz}{h}\end{array}\right.$$

(18)

The solution does not depend on the angle $\varphi$, which is due to the material isotropy in the plane of the crystallographic axis $X_c{Y}_c$. The absolute value of the lateral response ($\sqrt{u_x^2+u_y^2}$) does not depend on the angle of the plate rotation around the Z-axis ($\psi$ angle).

For the case of the isotropic elastic material according to (15), non-zero components of the piezoelectric constant tensor in the laboratory coordinate system can be derived through the piezoelectric coefficients as follows:

$$\begin{array}{l}{e}_{15}^0=\frac{d_{15}^{0}Y}{2+2\nu }\\ e_{31}^0=-\frac{Y\left(d_{31}^0+d_{33}^0\nu \right)}{-1+\nu +2{\nu }^2}\\ e_{33}^0=\frac{Y\left(d_{33}^0\left(-1+\nu \right)-2d_{31}^0\nu \right)}{-1+\nu +2{\nu }^2}\end{array}$$

(19)

Thus, solution (18) can be represented through the piezoelectric coefficients as:

$$\left\{\begin{array}{c}{u}_x\left[z,\varphi ,\theta ,\psi \right]=\mathrm{Sin}\left[\theta \right]\mathrm{Sin}\left[\psi \right]\left(d_{31}^0-d_{33}^0+\left(d_{15}^0+d_{31}^0-d_{33}^0\right)\mathrm{Cos}\left[2\theta \right]\right)·\frac{Uz}{h}\\ u_y\left[z,\varphi ,\theta ,\psi \right]=\mathrm{Sin}\left[\theta \right]\mathrm{Cos}\left[\psi \right]\left(d_{31}^0-d_{33}^0+\left(d_{15}^0+d_{31}^0-d_{33}^0\right)\mathrm{Cos}\left[2\theta \right]\right)·\frac{Uz}{h}\\ u_z\left[z,\varphi ,\theta ,\psi \right]=\mathrm{Cos}\left[\theta \right]\left(\frac{d_{15}^0+d_{31}^0+d_{33}^0-2\nu (d_{15}^0-d_{31}^0)}{4\left(\nu -1\right)}\right.+\hfill \\ \left.+\frac{\left(d_{15}^0+d_{31}^0-d_{33}^0\right)\left(2\nu -1\right)\mathrm{Cos}\left[2\theta \right]}{2\left(\nu -1\right)}\right)·\frac{Uz}{h}\end{array}\right.$$

(20)

It should be noted that the lateral response in the uniform electric field approximation in the form (20) does not depend on the material’s elastic properties and has the only contribution from the piezoelectric coefficients.

3.3. Barium Titanate: Anisotropic Elastic Properties

The tensor of the anisotropic elastic modules for the crystal with the 4 mm symmetry can be written in the laboratory coordinate system as:

$$c_{ij}^0=c_{4mm}=\left(\begin{array}{cc}\begin{array}{ccc}{c}_{11}^0& c_{12}^0& c_{13}^0\\ c_{12}^0& c_{11}^0& c_{13}^0\\ c_{13}^0& c_{13}^0& c_{33}^0\end{array}& \begin{array}{ccc}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array}\\ \begin{array}{ccc}0 & 0& 0\\ 0 & 0& 0\\ 0 & 0& 0\end{array}& \begin{array}{ccc}{c}_{44}^0& 0& 0\\ 0& c_{44}^0& 0\\ 0& 0& \frac{1}{2}\left(c_{11}^0-c_{12}^0\right) \end{array}\end{array}\right)$$

(21)

Equation (12) for the $\overrightarrow{u}$ displacement vector components can be written for the anisotropic barium titanate plate rotated at random angles $\left\{\varphi ,\theta ,\psi \right\}$ as:

$$\left\{\begin{array}{c}{u}_x\left[z,\varphi ,\theta ,\psi \right]=\frac{\mathrm{Sin}\left[\theta \right]\left(\begin{array}{c}(3{\mathrm{c}}_{33}^0\left({\mathrm{e}}_{15}^0+{\mathrm{e}}_{31}^0\right)+\left({\mathrm{c}}_{13}^0+2{\mathrm{c}}_{44}^0\right)\left({\mathrm{e}}_{31}^0-3{\mathrm{e}}_{33}^0\right)-\\ -{\mathrm{c}}_{11}^0\left(3{\mathrm{e}}_{15}^0+{\mathrm{e}}_{33}^0\right)+4({\mathrm{c}}_{11}^0{\mathrm{e}}_{15}^0+{\mathrm{c}}_{33}^0\left({\mathrm{e}}_{15}^0+{\mathrm{e}}_{31}^0\right)-\\ -\left({\mathrm{c}}_{13}^0+2{\mathrm{c}}_{44}^0\right){\mathrm{e}}_{33}^0)\mathrm{Cos}\left[2 \theta \right]-(-{\mathrm{c}}_{33}^0\left({\mathrm{e}}_{15}^0+{\mathrm{e}}_{31}^0\right)+\\ +{\mathrm{c}}_{11}^0\left({\mathrm{e}}_{15}^0-{\mathrm{e}}_{33}^0\right)+\left({\mathrm{c}}_{13}^0+2{\mathrm{c}}_{44}^0\right)\left({\mathrm{e}}_{31}^0+{\mathrm{e}}_{33}^0\right))\mathrm{Cos}\left[4\theta \right])\end{array}\right)}{\left(\begin{array}{c}-{\left(c_{13}^0\right)}^2+c_{11}^0{c}_{33}^0-2c_{13}^0{c}_{44}^0+3\left(c_{11}^0+c_{33}^0\right)c_{44}^0+\\ +4\mathrm{Cos}\left[2\theta \right]\left(-c_{11}^0+c_{33}^0\right)c_{44}^0+\\ +\mathrm{Cos}\left[4\theta \right]\left({\left(c_{13}^0\right)}^2-c_{11}^0{c}_{33}^0+\left(c_{11}^0+2c_{13}^0+c_{33}^0\right)c_{44}^0\right)\end{array}\right)}·\frac{\mathrm{Sin}\left[\psi \right]Uz}{4h}\\ u_x\left[z,\varphi ,\theta ,\psi \right]=\frac{\mathrm{Sin}\left[\theta \right]\left(\begin{array}{c}3c_{33}^0\left({\mathrm{e}}_{15}^0+{\mathrm{e}}_{31}^0\right)+\left({\mathrm{c}}_{13}^0+2{\mathrm{c}}_{44}^0\right)\left({\mathrm{e}}_{31}^0-3{\mathrm{e}}_{33}^0\right)-\\ -{\mathrm{c}}_{11}^0\left(3{\mathrm{e}}_{15}^0+{\mathrm{e}}_{33}^0\right)+4({\mathrm{c}}_{11}^0{\mathrm{e}}_{15}^0+{\mathrm{c}}_{33}^0\left({\mathrm{e}}_{15}^0+{\mathrm{e}}_{31}^0\right)-\\ -\left({\mathrm{c}}_{13}^0+2{\mathrm{c}}_{44}^0\right){\mathrm{e}}_{33}^0)\mathrm{Cos}\left[2\theta \right]-({\mathrm{c}}_{33}^0\left({\mathrm{e}}_{15}^0+{\mathrm{e}}_{31}^0\right)-\\ -{\mathrm{c}}_{11}^0\left({\mathrm{e}}_{15}^0-{\mathrm{e}}_{33}^0\right)-\left({\mathrm{c}}_{13}^0+2{\mathrm{c}}_{44}^0\right)\left({\mathrm{e}}_{31}^0+{\mathrm{e}}_{33}^0\right))\mathrm{Cos}\left[4\theta \right]\end{array}\right)}{\left(\begin{array}{c}-{\left(c_{13}^0\right)}^2+c_{11}^0{c}_{33}^0-2c_{13}^0{c}_{44}^0+3\left(c_{11}^0+c_{33}^0\right)c_{44}^0+\\ +4\mathrm{Cos}\left[2\theta \right]\left(-c_{11}^0+c_{33}^0\right)c_{44}^0+\\ +\mathrm{Cos}\left[4\theta \right]\left({\left(c_{13}^0\right)}^2-c_{11}^0{c}_{33}^0+\left(c_{11}^0+2c_{13}^0+c_{33}^0\right)c_{44}^0\right)\end{array}\right)}·\frac{\mathrm{Cos}\left[\psi \right]Uz}{h}\\ u_z\left[z,\varphi ,\theta ,\psi \right]=\frac{\left(\begin{array}{c}-c_{33}^0{e}_{15}^0-c_{33}^0{e}_{31}^0-2c_{44}^0{e}_{31}^0-2c_{44}^0{e}_{33}^0-\\ \left.-c_{11}^0\left(3e_{15}^0+e_{33}^0\right)\right)+c_{13}^0\left(4e_{15}^0+e_{31}^0+e_{33}^0\right)+\\ +4\left({\mathrm{c}}_{11}^0{\mathrm{e}}_{15}^0-{\mathrm{c}}_{13}^0{\mathrm{e}}_{15}^0+{\mathrm{c}}_{44}^0\left({\mathrm{e}}_{31}^0-{\mathrm{e}}_{33}^0\right)\right)\mathrm{Cos}\left[2\theta \right]+\\ +({\mathrm{c}}_{33}^0\left({\mathrm{e}}_{15}^0+{\mathrm{e}}_{31}^0\right)-{\mathrm{c}}_{11}^0\left({\mathrm{e}}_{15}^0-{\mathrm{e}}_{33}^0\right)-\\ -\left({\mathrm{c}}_{13}^0+2{\mathrm{c}}_{44}^0\right)\left({\mathrm{e}}_{31}^0+{\mathrm{e}}_{33}^0\right))\mathrm{Cos}\left[4\theta \right])\end{array}\right)}{\left(\begin{array}{c}-{\left(c_{13}^0\right)}^2+c_{11}^0{c}_{33}^0-2c_{13}^0{c}_{44}^0+3\left(c_{11}^0+c_{33}^0\right)c_{44}^0+\\ +4\mathrm{Cos}\left[2\theta \right]\left(c_{33}^0-c_{11}^0\right)c_{44}^0+\\ +\mathrm{Cos}\left[4\theta \right]\left({\left(c_{13}^0\right)}^2-c_{11}^0{c}_{33}^0+\left(c_{11}^0+2c_{13}^0+c_{33}^0\right)c_{44}^0\right)\end{array}\right)}·\mathrm{Cos}\left[\theta \right]\frac{Uz}{h}\end{array}\right.$$

(22)

The anisotropic case’s solution becomes significantly more complicated, which is sourced by the diversification of the elastic properties in different crystallographic directions. By analogy with Equation (19), piezoelectric modules relate to the piezoelectric coefficients (21) as:

$$\begin{array}{l}{e}_{15}^0=d_{15}^0{c}_{44}^0\\ e_{31}^0=d_{33}^0{c}_{13}^0+d_{31}^0\left(c_{11}^0+c_{12}^0\right)\\ e_{33}^0=2d_{31}^0{c}_{13}^0+d_{33}^0{c}_{33}^0\end{array}$$

(23)

Substitution of the (23) to the displacement vector $\overrightarrow{u}$ (22) leads to a significant complication of the equations, making it challenging to represent and intuitively interpret it in a general form.

3.4. Analysis of the Solution for the Barium Titanate: Anisotropic Elastic Properties

To illustrate the behavior of the solutions (18) and (22), they were analyzed for three orientations (cuts) of barium titanate crystals: (001)_c, (011)_c, and (100)_c, which correspond to the (0, 0, 0), (180°, 45°,0), and (90°, 90°, 0) values of the Euler angles, respectively. The solutions (18), (20), (22), and (22) through (23) for the displacement vector of the differently oriented crystals (chosen angles $\theta$, $\psi =0,$ and $\forall \varphi$) are summarized in

Table 1. The expressions of the displacement vector for the case of the differently oriented crystals of the barium titanate derived from (18), (20), (22), (22) via (23).

$\theta$ $\overrightarrow{u}$	${\left(001\right)}_c$ $\left\{\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right\}$	${\left(100\right)}_c$ $\left\{\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right\}$	${\left(011\right)}_c$ $\left\{\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right\}$
(18)	$\left\{\begin{array}{c}0\\ 0\\ \frac{U\left(1+\nu \right)\left(1-2\nu \right)e_{33}^0}{Y\left(\nu -1\right)}\end{array}\right\}$	$\left\{\begin{array}{c}0\\ \frac{2U\left(1+\nu \right)e_{15}^0}{-Y}\\ 0\end{array}\right\}$	$\left\{\begin{array}{c}0\\ \frac{U(\nu +1)(e_{31}^0 - e_{33}^0)}{\sqrt{2}Y}\\ \frac{U\left(1+\nu \right)\left(1-2\nu \right)\left(2e_{15}^0+e_{31}^0+e_{33}^0\right)}{2\sqrt{2}Y\left(\nu -1\right)}\end{array}\right\}$
(20)	$\left\{\begin{array}{c}0\\ 0\\ \frac{2U\nu d_{31}^0}{\nu -1}-U{d}_{33}^0\end{array}\right\}$	$\left\{\begin{array}{c}0\\ -U{d}_{15}^0\\ 0\end{array}\right\}$	$\left\{\begin{array}{c}0\\ \frac{\left({\mathrm{d}}_{31}^0-{\mathrm{d}}_{33}^0\right)U}{\sqrt{2}}\\ \frac{U\left(\left(1-2\nu \right)d_{15}+\left(1+2\nu \right)d_{31}+d_{33}\right)}{2\sqrt{2}\left(\nu -1\right)}\end{array}\right\}$
(22)	$\left\{\begin{array}{c}0\\ 0\\ -\frac{U{e}_{33}^0}{c_{33}^0}\end{array}\right\}$	$\left\{\begin{array}{c}0\\ -\frac{U{e}_{15}^0}{c_{44}^0}\\ 0\end{array}\right\}$	$\left\{\begin{array}{c}\frac{\begin{array}{c}-U\left(\left(c_{13}^0+2c_{44}^0\right)\left(e_{31}^0-e_{33}^0\right)-c_{33}^0\left(e_{15}^0\left.-e_{31}^0\right)+\right.\right.\\ \left.+c_{11}^0\left(e_{15}^0+e_{33}^0\right)\right)\end{array}}{\sqrt{2}\left({\left(c_{13}^0\right)}^2+2c_{13}^0{c}_{44}^0-c_{33}^0{c}_{44}^0-c_{11}^0\left(c_{33}^0+c_{44}^0\right)\right)}\\ \frac{\begin{array}{c}U\left(c_{33}^0\left(e_{15}^0+e_{31}^0\right)\right.-\\ \left.\left.-c_{13}^0\left(2e_{15}^0\right.+e_{31}^0+e_{33}^0\right)+c_{11}^0\left(e_{15}^0+e_{33}^0\right)\right)\end{array}}{\sqrt{2}\left({\left(c_{13}^0\right)}^2+2c_{13}^0{c}_{44}^0-c_{33}^0{c}_{44}^0-c_{11}^0\left(c_{33}^0+c_{44}^0\right)\right)}\end{array}\right\}$
(22)	$\left\{\begin{array}{c}0\\ 0\\ -\frac{2c_{13}^{0}U{d}_{31}^0}{c_{33}^0}-U{d}_{33}^0\end{array}\right\}$	$\left\{\begin{array}{c}0\\ -U{d}_{15}^0\\ 0\end{array}\right\}$	$\left\{\begin{array}{c}\frac{\begin{array}{c}-U\left(-{\left(c_{33}^0\right)}^2{d}_{33}^0+4\left(c_{13}^0+2c_{44}^0\right)\left(c_{12}^0{d}_{31}^0-\right.\right.\\ \left.-2c_{13}^0{d}_{31}^0+c_{13}^0{d}_{33}^0\right)+c_{11}^0\left(-c_{44}^0\left(d_{15}^0-8d_{31}^0\right)+\right.\\ \left.+\left(c_{12}^0-2c_{13}^0+3c_{33}^0\right)d_{31}^0+\left(c_{13}^0-3c_{33}^0\right)d_{33}^0\right)+\\ +c_{33}^0\left(c_{44}^0{d}_{15}^0+3c_{12}^0{d}_{31}^0-d_{33}^0\left(c_{13}^0+8c_{44}^0\right)-\right.\\ \left.\left.-2c_{13}^0{d}_{31}^0\right)+{\left(c_{11}^0\right)}^2{d}_{31}^0\right)\end{array}}{2\sqrt{2}\left({\left(c_{13}^0\right)}^2+2c_{13}^0{c}_{44}^0-c_{33}^0{c}_{44}^0-c_{11}^0\left(c_{33}^0+c_{44}^0\right)\right)}\\ \frac{\begin{array}{c}U\left(c_{13}^0\left(2c_{44}^0{d}_{15}^0\right.\right.\left.+c_{12}^0{d}_{31}^0\right)+\\ +{\left(c_{13}^0\right)}^2\left(2d_{31}^0+d_{33}^0\right)+c_{33}^0\left(c_{44}^0{d}_{15}^0+c_{12}^0{d}_{31}^0\right)+\\ +c_{11}^0\left(c_{44}^0{d}_{15}^0+\left(c_{13}^0+c_{33}^0\right)d_{31}^0+c_{33}^0{d}_{33}^0\right.)\end{array}}{\sqrt{2}\left({\left(c_{13}^0\right)}^2+2c_{13}^0{c}_{44}^0-c_{33}^0{c}_{44}^0-c_{11}^0\left(c_{33}^0+c_{44}^0\right)\right)}\end{array}\right\}$

. Equation (22) is well-known as the Lefki and Dormans solution [38]. The expression has the most simple shape for the case of the crystal cut along the possible polarization directions, while it becomes much more complicated for the case of (011)_c crystal orientation, where a combination of the elastic modulus contributes.

Further, the response was analyzed in dependence on the angles $\theta$ and $\psi$(Figure 2), which was the equivalent in the “global excitation mode” PFM experiments to the choice of the specific cut of the crystal at the angle $\theta$, or rotation of the sample under SPM probe at angle $\psi$, respectively. (100)_c cut of the crystal most effectively illustrates piezoresponse behavior with the crystal rotation (Figure 2a,b). The same as it is in the isotropic approximation, the rotation of the (100)_c oriented plate around the axis of the normal to the surface conserves the absolute value of the lateral piezoresponse, while $u_x$ and $u_y$ components transform by the cosine and sinus laws, respectively. It means that PFM in the “global excitation mode” should be gradually redistributed from the twisting to the buckling component of the cantilever motion [10]. As expected, normal piezoresponse, i.e., vertical PFM signal, was absent on this crystal cut. Comparison between anisotropic and isotropic cases revealed that the piezoresponse in the anisotropic case was around 20 percent lower. However, the general trend of piezoresponse behavior conserves for both cases (Figure 2a,b).

As expected, the normal piezoresponse in (001)_c had a contribution not only from the d₃₃ coefficient but from the d₃₁ as well (Table 1). In the case of the anisotropic elastic properties, contribution to the displacement vector from the d₃₁ coefficient was proportional to the ratio between c₃₃ and c₁₃ elastic coefficients, which was about 0.92. In contrast, in the isotropic case, the displacement vector was proportional to the ν/(ν − 1), which was about 0.67. Thus, the deviation of the response due to d₃₁ was more significant for the case of the anisotropic elastic properties.

The plot of the surface displacement $\overrightarrow{u}$ vector components for the (010)_c cut of the crystal ($\theta =45°$) depending on the angle $\psi$ is presented in Figure 2c,d. In this case, the lateral response conserves the behavior, and a normal response appears, which was independent of the rotational angle. On this cut, the piezoresponse for the anisotropic crystal properties was more than twice as larger as the response calculated in the isotropic approximation. Other informative plots represent the dependencies of the lateral and vertical surface displacements on the angle θ (the angle between the outer normal to the plate and the crystallographic axis $Z_c$), which revealed a significant difference of the dependencies of the vertical and lateral piezoresponses in the isotropic and anisotropic cases (Figure 2e,f). Notably, not only the value of the piezoresponse changes after the transition from the anisotropic case to isotropic approximation but also the profile of the piezoresponse signal changes significantly.

3.5. Lithium Niobate Material Tensors

Lithium niobate is a rhombohedral material at room temperature, which belongs to 3m point group [54]. The tensor of piezoelectric constants $e_{ij}^0$ and coefficients $d_{ij}^0$ in the laboratory coordinate system has the form:

$$e_{ij}^0=\left(\begin{array}{cc}\begin{array}{ccc}0& 0& 0\\ -e_{22}^0& e_{22}^0& 0\\ e_{31}^0& e_{31}^0& e_{33}^0\end{array}& \begin{array}{ccc}0& e_{15}^0& -e_{22}^0\\ e_{15}^0& 0& 0\\ 0& 0& 0\end{array}\end{array}\right)$$

(24)

$$d_{ij}^0=\left(\begin{array}{cc}\begin{array}{ccc}0& 0& 0\\ -d_{22}^0& d_{22}^0& 0\\ d_{31}^0& d_{31}^0& d_{33}^0\end{array}& \begin{array}{ccc}0& d_{15}^0& -d_{22}^0\\ d_{15}^0& 0& 0\\ 0& 0& 0\end{array}\end{array}\right)$$

(25)

The anisotropic elastic properties for the lithium niobate were taken from [54], while isotropic elastic properties, Young’s modulus and Poisson ratio, were calculated using Voigt averaging the elastic tensor [53]. Piezoelectric and elastic properties used for the calculations are presented in Appendix B.

3.6. Lithium Niobate: Isotropic Elastic Properties

The displacement vector $\overrightarrow{u}$ for XY-infinite lithium niobate plate of the thickness h rotated at arbitrary Euler angles $\left\{\varphi ,\theta ,\psi \right\}$ can be written in the approximation of elastic properties isotropy (50) as:

$$\left\{\begin{array}{c}{u}_x\left[z,\varphi ,\theta ,\psi \right]=\left(1+\nu \right)\mathrm{Sin}\left[\theta \right](-\mathrm{Sin}\left[\theta \right](\mathrm{Cos}\left[\psi \right]\left(\mathrm{Sin}\left[\varphi \right]+3\mathrm{Sin}\left[3\varphi \right]\right)+\mathrm{Cos}\left[\theta \right](\mathrm{Cos}\left[\varphi \right]+\\ +3\mathrm{Cos}\left[3\varphi \right])\mathrm{Sin}\left[\psi \right])e_{22}^0+2\mathrm{Sin}\left[\psi \right]\left(e_{31}^0+\mathrm{Cos}\left[2\theta \right]\left(e_{15}^0+e_{31}^0-e_{33}^0\right)-e_{33}^0\right))·\frac{Uz}{Yh}\\ u_y\left[z,\varphi ,\theta ,\psi \right]=\mathrm{Cos}\left[\psi \right]\left(e_{31}^0-e_{33}^0+\left(2e_{15}^0+e_{31}^0-e_{33}^0\right)\mathrm{Cos}\left[2\theta \right]-e_{22}^0\mathrm{Cos}\left[3\varphi \left]\mathrm{Sin}\right[2\theta \right]\right) \\ ·\left(e_{22}^0\mathrm{Cos}\left[3\varphi \left]\mathrm{Sin}\right[2\theta \right]\right)+2{\mathrm{e}}_{22}^0\mathrm{Sin}\left[\theta \right]\mathrm{Sin}\left[3\varphi \right]\mathrm{Sin}\left[\psi \right])\frac{Uz\left(1+\nu \right)}{Yh}\\ u_z\left[z,\varphi ,\theta ,\psi \right]=\left(1+\nu \right)\left(2\nu -1\right)(\left(-2e_{15}^0-e_{31}^0-e_{33}^0\right)\mathrm{Cos}\left[\theta \right]+\left(2{\mathrm{e}}_{15}^0+{\mathrm{e}}_{31}^0-{\mathrm{e}}_{33}^0\right)\mathrm{Cos}\left[\theta \right]\mathrm{Cos}\left[2\theta \right]+\\ +2e_{22}^0\mathrm{Cos}\left[3\varphi \right]\mathrm{Sin}{\left[\theta \right]}^3)·\frac{Uz}{2Yh\left(\nu -1\right)}\end{array}\right.$$

(26)

Similar to the case of the barium titanate, the absolute value of the lateral response ($\sqrt{u_x^2+u_y^2}$) does not depend on the angle of plate rotation around the Z-axis (angle $\psi$) but depends on the other angles. According to (15), the non-zero components of the piezoelectric constant tensor in the elastically isotropic material can be expressed in the laboratory coordinate system through the components of the elastic moduli tensor as follows:

$$\begin{array}{l}\begin{array}{l}{e}_{15}^0=\frac{Y{d}_{15}^0}{2+2\nu }\\ e_{22}^0=\frac{Y{d}_{22}^0}{1+\nu } \end{array} \\ \begin{array}{l}{e}_{31}^0=-\frac{Y\left(d_{31}^0+\nu d_{33}^0\right)}{-1+\nu +2{\nu }^2} \\ e_{33}^0=\frac{Y\left(-2\nu d_{31}^0+\left(-1+\nu \right)d_{33}^0\right)}{-1+\nu +2{\nu }^2}\end{array}\end{array}$$

(27)

Therefore, the solution (26) is expressed through (27) as:

$$\left\{\begin{array}{c}{u}_x\left[z,\varphi ,\theta ,\psi \right]=-\mathrm{Sin}\left[\theta \right](\mathrm{Cos}\left[\psi \right]\mathrm{Sin}\left[\theta \right]\left(\mathrm{Sin}\left[\varphi \right]+3\mathrm{Sin}\left[3\varphi \right]\right)d_{22}^0+\mathrm{Sin}\left[\psi \right](\mathrm{Cos}\left[\theta \right]· \\ ·\left(\mathrm{Cos}\left[\varphi \right]+3\mathrm{Cos}\left[3\varphi \right]\right)\mathrm{Sin}\left[\theta \right]d_{22}^0-2d_{31}^0-\mathrm{Cos}\left[2\theta \right]\left(d_{15}^0+2d_{31}^0-2d_{33}^0\right)+2d_{33}^0))·\frac{Uz}{h}\\ u_y\left[z,\varphi ,\theta ,\psi \right]=(\mathrm{Cos}\left[\psi \right]\left(d_{31}^0-d_{33}^0+\left(d_{15}^0+d_{31}^0-d_{33}^0\right)\mathrm{Cos}\left[2\theta \right]-d_{22}^0\mathrm{Cos}\left[3\varphi \left]\mathrm{Sin}\right[2\theta \right]\right)+\\ +2d_{22}^0\mathrm{Sin}\left[\theta \right]\mathrm{Sin}\left[3\varphi \right]\mathrm{Sin}\left[\psi \right]))·\frac{Uz}{2h}\\ u_z\left[z,\varphi ,\theta ,\psi \right]=(\mathrm{Cos}\left[\theta \right]\left(d_{15}^0+d_{31}^0+d_{33}^0-2d_{15}^0\nu +2d_{31}^0\nu +\left(d_{15}^0+d_{31}^0-d_{33}^0\right)\left(-1+2\nu \right)\mathrm{Cos}\left[2\theta \right]\right)\\ +2d_{22}^0\left(-1+2\nu \right)\mathrm{Cos}\left[3\varphi \right]\mathrm{Sin}{\left[\theta \right]}^3)\frac{Uz}{2h\left(\nu -1\right)}\end{array}\right.$$

(28)

The lateral response in the form (28), as well as for barium titanate, does not depend on the elastic properties of the material.

3.7. Lithium Niobate: Anisotropic Elastic Properties

The tensor of elastic moduli for crystals of similar symmetry (3m) in the laboratory coordinate system:

$$c_{ij}^0=c_{3m}^0=\left(\begin{array}{cc}\begin{array}{ccc}{c}_{11}^0& c_{12}^0& c_{13}^0\\ c_{12}^0& c_{11}^0& c_{13}^0\\ c_{13}^0& c_{13}^0& c_{33}^0\end{array}& \begin{array}{ccc}{c}_{14}^0& 0& 0\\ -c_{14}^0& 0& 0\\ 0& 0& 0\end{array}\\ \begin{array}{ccc}{c}_{14}^0& -c_{14}^0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{ccc}{c}_{44}^0& 0& 0\\ 0& c_{44}^0& c_{14}^0\\ 0& c_{14}^0& \frac{1}{2}\left(c_{11}^0-c_{12}^0\right)\end{array}\end{array}\right)$$

(29)

The piezoelectric constants are expressed through the piezoelectric coefficients as:

$$\begin{array}{l}\begin{array}{l}{e}_{15}^0=c_{44}^0{d}_{15}^0-c_{14}^0{d}_{22}^0 \\ e_{22}^0=-c_{14}^0{d}_{15}^0+\left(c_{11}^0-c_{12}^0\right)d_{22}^0\end{array}\\ \begin{array}{l}{e}_{31}^0=c_{11}^0{d}_{31}^0+c_{12}^0{d}_{31}^0+c_{13}^0{d}_{33}^0\\ e_{33}^0=2c_{13}^0{d}_{31}^0+c_{33}^0{d}_{33}^0 \end{array}\end{array}$$

(30)

The expression for the displacement vector $\overrightarrow{u}$ in the fully anisotropic case is too tedious to present here in the analytical form. Thereby, the solution is given in the main text only for some specific chosen Euler angles.

3.8. Analysis of the Solution for the Lithium Niobate: Anisotropic Elastic Properties

Analytical expressions for the displacement vector of points of the plate surface ($z=h$) at various angles $\theta$ (at $\psi =0$ and $\forall \varphi$) for the lithium niobate are summarized in

Table 2. The expressions of the displacement vector for the case of the differently oriented crystals of the lithium niobate derived from (26), (27), (12) via (29), (12) via (29), and (30).

$\theta$ .. $\overrightarrow{u}$	${\left(001\right)}_{\mathrm{c}}$ $\left\{\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right\}$	${\left(100\right)}_c$ $\left\{\begin{array}{c}{u}_x\\ u_y\\ u_z\end{array}\right\}$
(26)	$\left\{\begin{array}{c}0\\ 0\\ \frac{\left(1+\nu \right)\left(1-2\nu \right)e_{33}^{0}U}{Y\left(\nu -1\right)}\end{array}\right\}$	$\left\{\begin{array}{c}\frac{2e_{22}^{0}U\left(1+\nu \right)}{Y}\\ -\frac{2e_{15}^{0}U\left(1+\nu \right)}{Y}\\ 0\end{array}\right\}$
(27)	$\left\{\begin{array}{c}0\\ 0\\ \frac{2U\nu d_{31}^0}{\nu -1}-U{d}_{33}^0\end{array}\right\}$	$\left\{\begin{array}{c}2d_{22}^{0}U\\ -d_{15}^{0}U\\ 0\end{array}\right\}$
(12) via (29)	$\left\{\begin{array}{c}0\\ 0\\ -\frac{U{e}_{33}^0}{c_{33}^0}\end{array}\right\}$	$\left\{\begin{array}{c}-\frac{2\left({\mathrm{c}}_{14}^0{\mathrm{e}}_{15}^0+{\mathrm{c}}_{44}^0{\mathrm{e}}_{22}^0\right)U}{2{\mathrm{c}}_{14}^0{}^2+\left({\mathrm{c}}_{12}^0-{\mathrm{c}}_{11}^0\right){\mathrm{c}}_{44}^0}\\ \frac{\left({\mathrm{c}}_{11}^0{\mathrm{e}}_{15}^0-{\mathrm{c}}_{12}^0{\mathrm{e}}_{15}^0+2{\mathrm{c}}_{14}^0{\mathrm{e}}_{22}^0\right)U}{2{\mathrm{c}}_{14}^0{}^2+\left({\mathrm{c}}_{12}^0-{\mathrm{c}}_{11}^0\right){\mathrm{c}}_{44}^0}\\ 0\end{array}\right\}$
(12) via (29) and (30)	$\left\{\begin{array}{c}0\\ 0\\ -\frac{2c_{13}^{0}U{d}_{31}^0}{c_{33}^0}-U{d}_{33}^0\end{array}\right\}$	$\left\{\begin{array}{c}\frac{2\left({\mathrm{c}}_{14}^0{}^2+\left(-{\mathrm{c}}_{11}^0+{\mathrm{c}}_{12}^0\right){\mathrm{c}}_{44}^0\right){\mathrm{d}}_{22}^{0}U}{2{\mathrm{c}}_{14}^0{}^2+\left({\mathrm{c}}_{12}^0-{\mathrm{c}}_{11}^0\right){\mathrm{c}}_{44}^0}\\ \frac{\left(-2{\mathrm{c}}_{14}^0{}^2{\mathrm{d}}_{15}^0+\left({\mathrm{c}}_{11}^0-{\mathrm{c}}_{12}^0\right){\mathrm{c}}_{44}^0{\mathrm{d}}_{15}^0+\left({\mathrm{c}}_{11}^0-{\mathrm{c}}_{12}^0\right){\mathrm{c}}_{14}^0{\mathrm{d}}_{22}^0\right)U}{2{\mathrm{c}}_{14}^0{}^2+\left({\mathrm{c}}_{12}^0-{\mathrm{c}}_{11}^0\right){\mathrm{c}}_{44}^0}\\ 0\end{array}\right\}$

. Analysis of the piezoresponse dependencies on the angle $\psi$ revealed a similar trend as in barium titanate, but without such a large difference between isotropic and anisotropic cases. It did not exceed 3% for various orientations. It was also important that (100)_c and (010)_c oriented crystals possessed equivalent $\psi$-angle dependencies in the tetragonal barium titanate, while the dependencies were different for the lithium niobate, where piezoresponse did not nullify on (100)_c cut at the angle $\psi =0$ due to rhombohedral symmetry of the crystalline lattice (Figure 3a,b). The dependencies on the cut angle $\theta$ were also significantly different for the case of lithium niobate in comparison to barium titanate (Figure 2e,f and Figure 3e,f). In the lithium niobate case, a negligible difference between isotropic and anisotropic cases was observed, which demonstrated that the isotropic approximation’s applicability is highly dependent on the symmetry and elastic properties of the crystal.

Table 3. The ratios of the contributions to the piezoresponse in the lithium niobate from each of the piezoelectric tensor components.

	(001)c				(011)c
	d15	d22	d31	d33	d15	d22	d31	d33
ux	0	0	0	0	0	0	0	0
uy	0	0	0	0	0	157%	−4%	−53%
uz	0	0	−5%	105%	60%	26%	−2%	16%
	(010)c				(100)c
	d15	d22	d31	d33	d15	d22	d31	d33
ux	0	0	0	0	0	100%	0	0
uy	101%	−1%	0	0	96%	4%	0	0
uz	0	100%	0	0	0	0	0	0

shows the ratios between the contributions to the piezoelectric response from each component of the piezoelectric coefficient tensor calculated for lithium niobate crystals with the different orientations. The total measured piezoresponse was taken as 100%. Negative values, thus, indicate that the given coefficient reduced the total piezoelectric response. For example, a major contribution to the piezoresponse in (001)_c cut lithium niobate comes from the d₃₃ coefficient, while the contribution from the d₃₁ coefficient reduces the total electromechanical response. d₁₅ and d₂₂ coefficients determine the response in the (100)_c− and (010)_c nonpolar cuts. In contrast, the (011)_c cut of the crystal has a valuable contribution from all piezoelectric coefficients of the lithium niobate piezoelectric tensor, which makes difficult the recovery of piezoelectric coefficients in this orientation from the values of the measured piezoresponse.

4. Verification of the Analytical Model

4.1. Verification of the Analytical Model by the Finite Element Modeling (FEM)

The analytical results were verified by the quantitative comparison with the numerical solution of the same problem with FEM in the COMSOL Multiphysics software. The model was realized using “Electrostatics” and “Structural mechanics” modules. The piezoelectric sample of a finite-size was clamped to the rigid non-deformable steel plate. The constant potential was applied to the top surface, while the bottom surface of the piezoelectric plate was grounded. The rotation of the plate was performed with the same Euler matrix as used in the analytical theory (Appendix C). The comparison of the results was made for the (010)_c cut of the lithium niobate single crystal. The same properties of the crystal were used in the FEM and analytical model.

The image illustrating the piezoelectric displacement in the electric field for the (010)_c cut of the lithium niobate single crystal is presented in Figure 4a. It is seen that shear deformation due to the d₁₅ shear coefficient dominates against other contributions. The top surface moved along x-direction in the sample Cartesian coordinate system (polarization direction in (010)_C crystal cut). At the same time, the bottom part of the plate was constrained to the substrate. The uniform electric field led to the uniform displacement of the piezoelectric surface almost in the whole area of the sample. The non-uniformity of the strain could be observed only immediately in the vicinity of the sample edges, which is related to the contribution of the sample edges in the fixed-size sample (enhanced electric field on the sample edges). This non-uniformity is not crucial if measurements are done far from the sample edges (Figure 4b). However, the electrode deposition in the experimental realization is also usually imperfect near the boundaries. Thus, the measurements near the sample edges should be avoided anyhow. Thereby, FEM verified well the validity of the one-dimensional approximation in the solution of the problem. The piezoresponse angular dependencies calculated with the analytical solution and FEM (data from the middle part of the sample) revealed perfect coincidence with each other (Figure 5).

4.2. Experimental Verification of the Analytical Model

Experimental verification of the theoretical results was performed in (010)_c cut crystal of the lithium niobate (Yamaju Ceramics, Owariasahi, Japan). 50-nanometer-thick copper electrodes were deposited on both crystal surfaces using magnetron sputtering. The sample was glued onto the metal disc by the silver paint (Figure 6a). “Global-excitation” mode PFM measurements were realized with the NTEGRA Aura scanning probe microscope (NT-MDT, Russia) using an external Zurich Instruments HF2LI lock-in amplifier. Commercial uncoated Multi−75EG (Budget Sensors, Bulgaria) cantilevers were utilized for the piezoresponse signal detection. The scheme illustrating experimental conditions is presented in Figure 6a. 5 V amplitude AC was applied to the bottom electrode while the top electrode and AFM tip were grounded. PFM measurements were done at low frequency (400 Hz) on the resonance-free plateau. The phase offset of the lock-in was kept at zero for all the measurements, and X = Rcosφ was measured in dependence on the rotational angle ψ. To calibrate vertical and lateral piezoresponse “inverted optical level sensitivity” calibration factors were extracted from the vertical and lateral force-distance curves’ measurements [55,56]. The values of the piezoresponse were divided by the amplitude of the applied AC voltage and by the respective vertical and lateral shape factors of the cantilever extracted using the procedure described in [56,57]. Lateral response angular dependence was measured and quantified (Figure 6b) to compare theoretical results with the experiments, while the vertical signal was captured in the 180-degree angular position, where the contribution of the cantilever buckling was minimal [22]. A pure vertical displacement was evaluated to be −14.5 pm/V. Both vertical and lateral measured responses matched well the results of the analytical model.

5. Conclusions

In this work, piezoelectric response in the uniform electric field created by the top electrode was analyzed for the case of the arbitrarily oriented piezoelectric material bonded to the rigid substrate, which presents “global excitation mode” PFM measurements. As the measurements of the piezoresponse angular dependencies often become crucial for interpreting the material piezoelectric properties [12], the particular focus in this paper was on the analysis of the angular dependencies. Based on the results of differently oriented barium titanate and lithium niobate crystals, the piezoelectric response was found to contain several contributions to the piezoelectric coefficients, which are often ignored in the measurements. The displacements were calculated for the case of fully coupled conditions both for the anisotropic and isotropic elastic properties. The anisotropic case was compared with the often considered transversally isotropic approximation. The difference between these two cases was shown to be significant for the case of barium titanate crystals and insignificant for lithium niobate. Thus, the isotropic approximation can be safely used only in specific crystals. The derived analytical theory was verified by finite element modeling and the global excitation mode PFM experiment demonstrated the very well coincidence.

Additionally, the solutions were represented in the Wolfram Mathematica interactive code, allowing easy access for researchers to the developed theoretical methodology. The code was shared in the Wolfram Notebook Archive cloud service and is available for the researcher community. The analytical model and interactive program code can be used to predict the piezoresponse angular behavior in various piezoelectric materials. This study is important for the further development of the quantitative “global excitation mode” PFM measurements, especially in low-dimensional materials, such as 2D and 1D piezoelectric nanostructures.