Density Functional Perturbation Theory to Predict Piezoelectric Properties

Abstract

Among the various computational methods in materials science, only first-principles calculation based on the density functional theory has predictability for unknown material. Especially, density functional perturbation theory (DFPT) can effectively calculate the second derivative of the total energy with respect to the atomic displacement. By using DFPT method, we can predict piezoelectric constants, dielectric constants, elastic constants, and phonon dispersion relationship of any given crystal structure. Recently, we established the computational technique to decompose piezoelectric constants into each atomic contribution, which enable us to gain deeper insights to understand the piezoelectricity of material. Therefore, in this chapter, we will introduce the computational framework to predict piezoelectric properties of polar material by means of DFPT and details of decomposition technique of piezoelectric constants. Then, we will show some case studies to predict and discover new piezoelectric material.

1. Introduction

In this chapter, we will introduce how recent computational techniques can successfully predict response properties, represented as piezoelectricity, by means of perturbation method. Piezoelectricity is the polarization change in response to external mechanical force. Inversely, if electrical field is applied to piezoelectric material, mechanical strain is induced (inverse piezoelectric effect). Therefore, piezoelectric materials are widely used as vibrational censors, surface acoustic wave devices, and actuators. Only the material having no inversion symmetry shows piezoelectricity. For example, Figure 1 shows schematic illustration of the piezoelectric effect. Positions of positively charged ion (cation) and negatively charged ion (anion) are represented as plus and minus symbols. Figure 1a shows the paraelectric phase, where ions are orderly located with inversion symmetry. On the other hand, ions are slightly displaced by δ with respect to those in paraelectric phase, as shown in Figure 1b. Such small displacement induces microscopic polarization P_s along the ionic displaced direction.

Figure 1.

Ionic configuration of (a) paraelectric phase and (b) ferroelectric phase. (c) Ionic displacement according to the external force.

Because ferroelectric phase is energetically more stable than paraelectric phase under low temperature, P_s is frequently referred as the spontaneous polarization. Above Curie temperature, ferroelectric properties are disappeared since paraelectric phase becomes more stable than ferroelectric one. Figure 1c shows the schematic illustration of the principal of piezoelectricity, where external stress (red-colored arrows) increases the ionic displacement and resultant polarization. In this case, external stress increases the spontaneous polarization by ΔP_s = P_s^’−P_s. Therefore, piezoelectric constant is defined as the derivative of the spontaneous polarization with respect to the external field. More detailed and comprehensive description of piezoelectricity is reviewed by Martin [1].

First-principles calculation based on density functional theory (DFT [2, 3]) has been widely utilized as the computational method to predict the electronic properties of material under the ground state. Ideally, required information to conduct the first-principles calculation is only the crystal structure, including atomic species and position of periodic/nonperiodic structure unit. The most significant advantage of first-principles calculation is its predictability. Since King-Smith and Vanderbilt showed the theoretical methodology to calculate change in polarization per unit volume ΔP [4], dielectric and piezoelectric properties of wide range of materials in which electronic correlations are not too strong [5, 6, 7] have been accurately predicted. The derivative of total energy determines various properties. For example, determined forces, stresses, dipole moment (first-order derivatives), dynamical matrix, elastic constants, dielectric and piezoelectric constants (second-order derivative), nonlinear dielectric susceptibility, phonon–phonon interaction and Grüneisen parameters (third-order derivative), and so on. Thus, first-principles calculation has been made use of calculating the perturbed total energy of materials because of its accuracy. Although perturbations were made by hand up to the early 1980s, sophisticated methodology of density functional perturbation theory (DFPT) was proposed in 1987 by Baroni et al. [5]. They showed general formulation of total energy change with respect to atomic displacement and opened the way to efficiently compute the energy derivative with respect to the perturbation [5]. DFPT can compute response properties directly arising from the perturbations of strain, atomic displacement, and electric field by making use of linear response theory [8, 9, 10, 11]. Numbers of ferroelectric materials are theoretically investigated on the origin of their ferroelectric properties (including piezoelectricity and dielectric properties) by using DFPT. Because of technological importance, such theoretical researches have been focused on Pb-based perovskite material (e.g., PbTiO₃, PbZrO₃, and their solid solution [12, 13, 14, 15]) because they have excellent piezoelectric properties and are widely applied for actuators. However, due to the restriction of hazardous substance (RoHS) directive, researches on lead-free ferroelectric materials have gathered great attraction. By taking advantage of the predictability of DFPT, various lead-free ferroelectric oxide and nitride materials were theoretically investigated on their piezoelectric properties [16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. Moreover, DFPT calculations showed that piezoelectricity can be greatly enhanced by imposing isotropic stress for PbTiO₃ [26, 27], uniaxial stress for SrHfO₃ [28], uniaxial and biaxial strain for AlN-GaN solid solution alloy [29], and two-dimensional epitaxial strain for doped ZnO [30]. As latterly explained, those enhancements of piezoelectric constant are thought to be closely related to the phase transition. In the next section, we will show the definition of piezoelectric constants within the framework of DFPT.

2. Formulation of piezoelectric constants

Formulation and calculation methodologies to obtain response properties of materials in the framework of DFPT have been developed in a step-by-step manner, because degrees of freedom by perturbations of atomic displacement, homogeneous electric fields, and strain are often strongly coupled. For example, piezoelectricity affects elastic and dielectric properties. Therefore, special care must be paid for the calculation of coupled properties. In 2005, Hamann et al. demonstrated that elastic and piezoelectric tensors can be efficiently calculated by treating homogeneous strain within the framework of DFPT [31]. At the same time, Wu et al. systematically formulated response properties with respect to displacement, strain, and electric fields [32]. In this section, we will briefly introduce how piezoelectric properties are formulated in the framework of DFPT. In each formulation, Einstein implied-sum notation is used. Cartesian directions {x, y, z} are represented as α and β. Subscription of j and k = 1, …, 6 is the standard Voigt notation (represents directions of xx, yy, zz, yz, zx, and xy). The subscripts m and n are the degrees of freedom in the cell. They range from 1 to 3i, where i is the number of irreducible atoms because each atom has three degree of freedom along x, y, and z directions.

Total energy of material under perturbation of atomic displacement u, electric field σ, and strain η, E(u,σ,η), is defined as follows:

$$E\left(u,\sigma ,\eta \right)=\frac{1}{{\mathrm{\Omega }}_0}\left[E^0-\mathrm{\Omega }\sigma \bullet \mathrm{P}\right]$$E1

where E⁰ is the total energy of material under the ground state, Ω₀ is volume of the unit cell (smallest repeat unit of crystal), Ω is deformed volume of the unit cell, and P is the electric polarization. Following response functional tensor can be obtained by second-order differential of Eq. (1):

$$\text{Force constant matrix}:K_{\mathit{mn}}={\left.{\mathrm{\Omega }}_0\frac{{\partial }^{2}E}{\partial u_m\partial u_n}\right|}_{\sigma ,\eta }$$E2

$$\text{Clamped}-\mathrm{ion}\text{term of electric susceptibility}:{\overline{\mathrm{\chi }}}_{\mathit{\alpha \beta }}=-{\left.\frac{{\partial }^{2}E}{\partial {\sigma }_{\alpha }\partial {\sigma }_{\beta }}\right|}_{u,\eta }$$E3

$$\text{Clamped}-\mathrm{ion}\text{term of elastic tensor}:{\overline{C}}_{\mathit{jk}}=-{\left.\frac{{\partial }^{2}E}{\partial {\eta }_j\partial {\eta }_k}\right|}_{u,\sigma }$$E4

$$\text{Born effective charge tensor}:Z_{\mathit{m\alpha }}={\left.-{\mathrm{\Omega }}_0\frac{{\partial }^{2}E}{\partial u_m\partial {\sigma }_{\alpha }}\right|}_{\eta }$$E5

$$\text{Force}-\text{response internal}-\text{strain tensor}:{\Lambda }_{\mathit{mj}}={\left.-{\mathrm{\Omega }}_0\frac{{\partial }^{2}E}{\partial u_m\partial {\eta }_j}\right|}_{\sigma }$$E6

$$\text{Clamped}-\mathrm{ion}\text{piezoelectric tensor}:{\overline{e}}_{\mathit{\alpha j}}={\left.\frac{{\partial }^{2}E}{\partial {\sigma }_{\alpha }\partial {\eta }_j}\right|}_u$$E7

Clamped-ion term is a frozen quantity, which indicates that atomic coordinates are not allowed to relax as the homogeneous electric field or strain. Therefore, dynamical term should be added into the clamped-ion term in order to obtain proper response properties.

Simplest and physically well-understandable piezoelectric constant can be expressed as follows:

$$e_{\mathit{\alpha j}}=\frac{\partial P_{\alpha }}{\partial {\eta }_j}$$E8

In this expression, it is easily understood that piezoelectric e constant e_αj is a measure of the change in polarization induced by the external strain. As the atomic positions are changed according to the strain, change of the polarization includes both electronic contribution (clamped-ion term) and dynamical contribution (internal-strain term). The internal-strain term of piezoelectric constant is represented as follows:

$${\widehat{e}}_{\mathit{\alpha j}}=\frac{1}{{\mathrm{\Omega }}_0}{Z}_{\mathit{m\alpha }}{\left(K^{-1}\right)}_{\mathit{mn}}{\Lambda }_{\mathit{nj}}$$E9

Thus, proper piezoelectric constant can be obtained by Eqs. (7) and (9):

$$e_{\mathit{\alpha j}}={\left.\frac{{\partial }^{2}E}{\partial {\sigma }_{\alpha }\partial {\eta }_j}\right|}_u+\frac{1}{{\mathrm{\Omega }}_0}{Z}_{\mathit{m\alpha }}{\left(K^{-1}\right)}_{\mathit{mn}}{\Lambda }_{\mathit{nj}}$$E10

Here, the first and second terms on the right-hand side in Eq. (10) are the clamped-ion term and internal-strain term, respectively. The former shows the electronic contribution ignoring the atomic relaxation effect, and the latter shows the ionic contribution including the response of the atomic displacement to the strain. The Born effective charge Z_mα, force-constant matrix K_mn, and internal-strain tensor Λ_nj are the second derivatives of the energy with respect to the displacement and electric field, pairs of displacements, and displacement and strain, respectively. The internal-strain term of the piezoelectric stress constants can be further decomposed into the individual atomic contributions when the above second-derivative tensors are fully obtained.

On the other hand, the internal-strain term of the piezoelectric stress constant e_αj is frequently described by the following equation, using the Born effective charge Z_αβ and displacement u_β of each atom in the calculation cell:

$${\widehat{e}}_{\mathit{\alpha j}}=Z_{\mathit{\alpha \beta }}\frac{\partial u_{\beta }}{\partial {\eta }_j}$$E11

where $\partial u_{\beta }/\partial {\eta }_j$ shows the response of the first-order atomic displacement to the first-order strain. In this expression, the meaning of the piezoelectric stress constant, i.e., e_j is a measure of the change in polarization induced by the external strain, is much more visible than in Eq. (9). In the DFPT formalism, $\partial u_{\beta }/\partial {\eta }_j$ is implicitly calculated as a displacement-response internal-strain tensor Γ as follows [32]:

$${\Gamma }_{\mathit{nj}}={\Lambda }_{\mathit{mj}}{\left(K^{-1}\right)}_{\mathit{mn}}$$E12

Because the subscript n in Γ_nj indicates the degrees of freedom, Γ_nj can be decomposed into the individual atomic components, which also enables to calculate individual contribution of each atom for total piezoelectric constant.

Here, piezoelectric e constant defined as Eq. (9) is frequently referred as “piezoelectric strain constant.” On the other hand, it is much more natural and easy to control the stress (electric field) than to control the strain in any case. In this case, the piezoelectric strain constant d_αj is usually measured. It can be obtained from piezoelectric strain constant e_αj using the following relation:

$$d_{\mathit{\alpha j}}=s_{\mathit{jk}}{e}_{\mathit{\alpha k}}$$E13

where s_jk is the elastic compliance, which is given by the inverse matrix of the elastic constants C_jk.

Those formulations are implemented in specific first-principles simulation packages such as ABINIT [33] and Vienna ab initio simulation package (VASP) [34], and piezoelectric constants can be calculated on a daily basis. From the next section, we will show how DFPT calculation precisely gives piezoelectric properties of ferroelectric materials.

3. Introduction of target material

In this chapter, we selected LiNbO₃ as a target material to show the predictability of DFPT calculation. LiNbO₃ is one of ferroelectric materials and widely used as surface acoustic wave (SAW) and optical waveguide elements. Crystal structure of LiNbO₃, which belongs to the space group of R3c, is frequently referred as “strained perovskite structure.” Schematic illustrations of crystal structure of ABO₃ perovskite and LiNbO₃ are shown in Figure 2.

Figure 2.

Crystal structures of (a) ferroelectric phase with R3c symmetry and (b) paraelectric phase with R$\overline{\mathit{3}}$c symmetry LiNbO₃. Yellow green-, green-, and red-colored balls represent Li, Nb, and oxygen atoms, respectively. Bonding structures between Nb and surrounding oxygen atoms are represented as green-colored polyhedron. Two orthogonal crystallographic directions are shown as both a- and c-axes.

Crystal structures shown in the present chapter was visualized by using VESTA software [35]. Curie temperature of LiNbO₃ is quite high and ranges from 1140 [36] to 1210°C dependent on the quality of sample (variation of Li/Nb relation can shift Curie temperature [37]). Below Curie temperature, ferroelectric phase with R3c symmetry (crystal structure can be classified into 230 types of space group according to the symmetry group) shown in Figure 2a is stable. Paraelectric phase with R$\overline{\mathit{3}}$c symmetry, shown in Figure 2b, becomes stable above Curie temperature. In the paraelectric phase, it can be seen that both Li and oxygen is positioned with the same height along c-axis, and the position of Nb is just the center between two oxygen layers along c-axis. On the other hand, both Li and Nb are shifted in ferroelectric R3c phase along downward direction of c-axis with respect to those in paraelectric R$\overline{\mathit{3}}$c phase.

Due to the different bonding nature between Li-O and Nb-O, atomic positions of Li and Nb are off-centered within oxygen layers along c-axis. This structural characteristic is the ferroelectric nature of LiNbO₃. One of the notable properties of LiNbO₃ is its high-curie temperature (~1400 K). However, piezoelectric properties of LiNbO₃ are not so much superior as compared with Pb-based perovskites. Crystal structure of piezoelectric ABO₃ perovskite is based on the cubic structure (of Pm3m symmetry), shown in Figure 3a.

Figure 3.

(a) Crystal structure of cubic ABO₃ perovskite and possible polarization directions. (b) Representative unstable vibrational mode of cubic ABO₃ perovskite showing as arrows.

Cubic lattice is symmetric and usually high-temperature phase, same as LiNbO₃. The “strained perovskite structure” expression for LiNbO₃ means that LiO₆ and NbO₆ polyhedron are largely rotated with respect to the cubic perovskite structure. However, because of the simple atomic configuration of cubic structure, atoms can be displaced along various directions and change crystalline symmetry as shown in Figure 3a. Crystalline lattice is vibrated (referred as phonon) under finite temperature. Some lattice vibrations along specific directions are unstable. This specific phonon is called as soft mode with imaginary frequency. In such case, atoms are displaced along unstable phonon mode to lower the total energy. For example, cooperative atomic displacement along [001] direction shown in Figure 3b (referred as Γ₁₅ mode) changes symmetry from cubic to tetragonal (of P4mm symmetry), which leads polarization along [001] direction. Thus, polarization direction of perovskite is not restricted and allowed to be changed. This characteristic rotational polarization direction is favorable for piezoelectricity because grains in polycrystalline material are oriented along various directions. Thus, careful controlling of crystal structure is essential to obtain superior piezoelectric properties.

The most convenient way to control and drastically change the crystal structure is imposing high pressure. Many compounds have found to be possible to form LiNbO₃-type structure under the high-pressure synthesis [38], and some of them were quenchable phase. For example, LiNbO₃-type structured ZnSbO₃ was successfully synthesized [39] under high pressure, and improvement of the spontaneous polarization is suggested by enhancement of the covalency of Sn site from first-principles simulation [40]. Moreover, high-pressure synthesized research on LiNbO₃-type structure is now extended to more complex compounds such as oxynitrides [41, 42].

The crystal structure of ABO₃ compound is generally determined by the balance between the ionic radius of A and B element, which is frequently referred as tolerance factor. Due to the small size of the Li ion with respect to the tolerance factor of LiNbO₃, this compound cannot form stably the popular perovskite structure under the ambient condition. On the other hand, we predicted the crystal structures of high-pressure phase of LiNbO₃ [43], which were not completely elucidated by experimental study [44]. Revealed structures are NaIO₃-type structure (Pnma) as room temperature high-pressure phase and apatite-like structure (P6₃/m) as high-temperature high-pressure phase. It should be noted that the NaIO₃-type structure is closely related with the popular GdFeO₃-type perovskite structure. The only difference between these structures is that A-site position and B-site position are inter-exchanged. Therefore, there seems to be a possible way to connect the perovskite structure and LiNbO₃-type structure.

In our previous theoretical study on high-pressure phase, analysis was mainly concerned with phase transition mechanism only from the viewpoint of subgroup symmetry and energy barrier [43]. It will be instructive to deal with this phase transition phenomenon from the viewpoint of lattice instability as discussed in the field of the ferroelectric instability analysis. In the following section, we will show investigation on the potential piezoelectric properties of LiNbO₃ with various hypothetical crystal structures by the method of DFPT, and possible phase transition mechanism will be discussed from the viewpoint of soft mode.

4. Computational methodology

First-principles calculation was performed by using VASP code [34]. Interactions between ion and electron were treated by projector augmented wave (PAW) method [45]. PBEsol functional [46] was used to approximate exchanges and correlate interactions of electrons, which can be used to reproduce the lattice constants of various materials [45]. Precise calculation on the lattice constant is essential to predict piezoelectric properties because they depend on volume of unit cell Ω as shown in Eq. (1). The kinetic energy cutoff for plane waves was set at 500 eV, and the k-point mesh was set at ~0.03/Å intervals to obtain the converged total energy at less than 0.1 meV/atom. Before calculating the piezoelectric constants, atomic positions and cell parameters were optimized until the forces on each atom and cell converged at below 5 × 10⁻⁴ eV/Å.

Since VASP does not directly calculate Eq. (12), we added routine to calculate displacement-response internal-strain tensor Γ_nj and decompose piezoelectric constants into each atom. The sum of the decomposed piezoelectric constants was confirmed to accurately reproduce the total piezoelectric constants. Careful convergence tests with a higher energy cutoff and denser k-point mesh showed that the numerical accuracy of the calculated Γ_nj was less than 0.01. It was confirmed that this error does not influence our discussion and conclusion. Moreover, it was confirmed that the values of Γ_nj obtained by the DFPT method were consistent with those calculated by the direct method, in which the strain-displacement relation of each ion was explicitly calculated.

On the basis of cubic Pm3m phase, lattice instability analysis was performed by phonon calculation utilizing phonopy code [47]. Force constant matrix shown in Eq. (2) was constructed by DFPT calculation implemented in VASP code combined with supercell approach. Supercell was constructed by using unit cell so that orthogonal three axes of the supercell exceed 10 Å. Note that although supercell is not required in DFPT approach, the present VASP code implements perturbation at the zone center.

5. Calculated piezoelectric properties of LiNbO₃

Calculated piezoelectric properties of LiNbO₃ in ferroelectric phase are summarized in Table 1. Some experimentally measured values are also shown in Table 1. All properties are confirmed to be well reproduced by calculation. In a technological importance, 33 components are the most important because C-axis of LiNbO₃ is polarization direction. Calculated values of e₃₃, C₃₃, and ε₃₃ are especially well reproduced. It should be mentioned here that chemical composition of LiNbO₃ used for experiment is congruent and includes Li vacancy. On the other hand, calculation was performed by using stoichiometric LiNbO₃.

	Calculated value	Experimental value
Piezoelectric stress constant (C/m^2)
e15	3.73	3.655 ± 0.022 [48], 3.7 [49]
e22	2.51	2.407 ± 0.015 [48], 2.5 [49]
e31	0.21	0.328 ± 0.032 [48], 0.2 [49]
e33	1.69	1.894 ± 0.054 [48], 1.3 [49]
Elastic constant (GPa)
C11	190.7	198.86 ± 0.033 [48], 203 [49]
C12	58.3	54.67 ± 0.04 [48], 53 [49]
C13	62.4	67.99 ± 0.55 [48], 75 [49]
C14	13.5	7.83 ± 0.02 [48], 9 [49]
C33	220.0	234.18 ± 0.75 [48], 245 [49]
C44	49.2	59.85 ± 0.01 [48], 60 [49]
Dielectric constant
ε11	40.6	44.9 ± 0.4 [48], 44 [49]
ε33	24.1	26.7 ± 0.3 [48], 29 [49]

Table 1.

Piezoelectric constant, elastic constant, and dielectric constant calculated by DFPT and experimentally measured values.

Thus, Li vacancy is considered to have negligible influence on the piezoelectric properties. Decomposed ionic contribution of piezoelectric strain constant e₃₃ is summarized in Table 2. Although the Born effective charge of Nb is larger than its formal charge +5e, displacement-response internal-strain constant of Nb is negative value. This indicates that piezoelectricity of LiNbO₃ is mainly dominated by displacement of Li. Born effective charge indicates a degree of polarization induced by atomic displacement and dominated by the change in the orbital hybridization. Although anomalously large Born effective charge is crucial for superior piezoelectric properties of perovskite ABO₃ materials [50], the present study of decomposition of piezoelectric constant shows that coupling degree between external strain and atomic displacement is also indispensable to understand the piezoelectric properties.

Decomposed e33’ (C/m^2)			Born effective charge Z33 (e)			Displacement-response internal-strain constant Γ33
Li	Nb	O	Li	Nb	O	Li	Nb	O
0.1	0.05	0.16	1.03	6.77	−2.6	0.67	−0.05	−0.21

Table 2.

Decomposed piezoelectric constants of LiNbO₃.

6. Piezoelectric properties of perovskite-LiNbO₃

Next, we will show how piezoelectric properties are affected by crystal structure, while chemical composition is kept as LiNbO₃. Various hypothetical crystal structures common for perovskite-type structure were constructed, and their energetic stabilities were examined by calculating enthalpy H = U + PV (U is total energy obtained by first-principles calculation, P is external pressure, and V is equilibrium volume under pressure P) as a function of external pressure. Imposing high pressure is most convenient method to modify crystal structure and find unexpected stable phase. The following eight types of phases were considered:

Cubic, Pm-3 m; tetragonal, P4mm; and rhombohedral, R-3 m.

LiNbO₃-ferroelectric phase, R3c, and LiNbO₃-paraelectric phase, R-3c.

Orthorhombic, Amm2 and Cmmm, and high-pressure phase, P63/m.

where names of space groups are used to distinguish each structure. Crystal structure of each phase is shown in Figure 4a. Polyhedra shown in Figure 4a correspond to Nb-centered bonding structure of Nb-O bondings. Figure 4b shows the enthalpy difference of each phases as a function of external pressure. Here, external pressure is assumed to be isotropic. Standard of enthalpy was set to be the enthalpy of most stable R3c phase under ambient condition. At the positive (compressive) pressure region, P6₃/m phase becomes stable above 21 GPa, which is close to the experimental phase transition pressure of 25 GPa [43]. Details of the phase transition behavior under high pressure are theoretically investigated in our previous work [44]. Unfortunately, P6₃/m phase is highly symmetric and shows no piezoelectricity. At the negative (expansive) pressure region, enthalpy difference becomes smaller as there is an increase of negative pressure except for R$\overline{\mathit{3}}$c and P6₃/m phases.

Figure 4.

(a) Schematic illustration of eight types of perovskite-structured LiNbO₃ and their space groups. (b) Enthalpy differences of each phase measured from the enthalpy of R3c phase as a function of pressure.

Imposing negative pressure can be achieved by solid solution with parent phase of larger lattice constant. At −6 GPa, P4mm phase becomes stable, while R3m and Amm2 phases become stable at −9 GPa. However, bond breaking occurs in Nb-O bonding above −6 GPa for P4mm phase. The same bond breaking occurs in Amm2 and R3m phases at −11GPa and −14 GPa, respectively. Thus, those phase transitions occur just before bond breaking.

Within the eight phases shown in Figure 4a, only P4mm, R3m, R3c, and Amm2 phases show piezoelectricity. Piezoelectric stress constant, elastic constant, and dielectric constant of P4mm, R3m, and Amm2 phases are compared with those of R3c phase in Table 3. Various piezoelectric properties are observed by each phase. Especially for P4mm and Amm2 phases, high e₃₃ and relatively low C₃₃ values are predicted, which are advantageous for large piezoelectric strain constant d₃₃. On the other hand, R3m phase was found to be unstable because following mechanical stability conditions of rhombohedral symmetry:

$$C_{11}+C_{12}>0,C_{33}>0,{\left(C_{11}+C_{12}\right)}^{\ast }{C}_{33}>2C_{13},C_{11}-C_{12}>0,C_{44}>0,{\left(C_{11}-C_{12}\right)}^{\ast }{C}_{44}>2C_{14}$$E14

are broken because of C₄₄ < 0.

	R3c	P4mm	R3m	Amm2
Piezoelectric stress constant (C/m^2)
e15	3. 73	1.14	5.10	0.64
e22	2.51	—	1.19	—
e31	0.21	0.46	0.24	0.80
e33	1.69	3.28	1.92	2.99
Elastic constant (GPa)
C11	190.7	297.2	203.0	321.8
C12	58.3	48.9	169.0	91.6
C13	62.4	77.7	90.6	92.8
C14	13.5	—	−42.1	—
C33	220.0	157.7	206.6	176.4
C44	49.2	39.5	−29.6	32.5
Dielectric constant
ε11	40.6	56.2	36.6	28.2
ε33	24.1	12.3	16.2	13.3

Table 3.

Piezoelectric constant, elastic constant, and dielectric constant of R3c, P4mm, R3m, and Amm2 phases calculated by DFPT.

Figure 5a and b show piezoelectric properties of P4mm phase as a function of pressure and corresponding volume of unit cell. Dotted lines indicate zero pressure states. Piezoelectric stress constant e₃₃ of P4mm phase shows parabolic behavior and maximum value at zero pressure state. On the other hand, elastic constant C₃₃ of P4mm phase continuously decreases as volume increases, because orbital hybridization of Nb-O bonding along polarization direction decreases as bond length increases. At the pressure of −6 GPa, C₃₃ of P4mm phase shows almost zero value. This indicates that Nb-O bonding is broken. Piezoelectric stress constant d₃₃ shown in Figure 5b increases as volume, because of increase of elastic compliance. Especially at the pressure of −5 GPa just before bond breaking, d₃₃ shows maximum value or approximately 1000 pC/N.

Figure 5.

(a) Piezoelectric stress constant e₃₃ and elastic constant C₃₃ and (b) piezoelectric strain constant d₃₃ of P4mm phase as a function of pressure and corresponding volume of unit cell.

This giant piezoelectric constant is almost comparable to that of PZT material [51]. Giant piezoelectric constant is understood as a result of phase instability in morphotropic phase boundary [52]. The same as P4mm phase of LiNbO₃, we revealed that ZnO also showed anomalously large piezoelectric constant just before phase transition [30].

Finally, we would like to show phase transition path between cubic perovskite structure and LiNbO₃ structure. Figure 6a shows the energy change of Pm3m phase as a function of Li displacement along <001>, <011>, and <111> directions. Pm3m phase is paraelectric phase. Because ferroelectricity and piezoelectricity of LiNbO₃ are dominated by off-centering and displacement of Li, respectively, phase transition from Pm3m phase is also expected to be occurred by Li displacement. Li displacement along <001>, <011>, and <111> directions induces tetragonal, orthorhombic, and rhombohedral phase transition from cubic phase. Figure 6a clearly shows that tetragonal phase transition from Pm3m phase to P4mm phase is the most energetically advantageous. Figure 6b shows the phonon dispersion curve of cubic Pm3m phase of LiNbO₃. Horizontal axis corresponds to sampling path along high symmetric reciprocal point (q-point). Within the whole Brillouin zone of reciprocal space, unstable phonon modes with imaginary phonon frequencies are observed. Here, imaginary phonon frequency is represented as negative value for convenience. Therefore, cubic Pm3m phase of LiNbO₃ is thermodynamically unstable and considered to show phase transition in accordance with specific phonon mode of imaginary frequency (referred as soft mode). Thus, modulated structures were constructed by imposing atomic displacement along normal modes at each symmetric q-points. Modulated structures were structurally relaxed, and their space group and energy change from cubic P4mm phase were investigated. Summary of such modulated structures are shown in Table 4. At Γ point, tetragonal phase transition along with Γ₁₅ soft mode of cubic phase shown in Figure 3b shows energy gain of −0.422 eV/formula unit (f.u.). On the other hand, it was found that modulation at R point gives more stable energy gain of −0.682 eV/f.u. In this case, R₂₅ soft mode induces phase transition from cubic Pm3m phase to R$\overline{\mathit{3}}$c phase shown in Figure 2b.

Figure 6.

(a) Energy change of Pm3m phase as a function of li displacement along <001>, <011>, and <111> directions. (b) Phonon dispersion curve of Pm3m phase.

q-point	Frequency (THz)	Structure	Space group	Energy gain (eV/f.u.)
Γ	−7.22	Tetragonal	P4mm	−0.422
X	−5.38	Orthorhombic	Pmma	−0.220
M	−6.56	Orthorhombic	Pmma	−0.125
R	−5.51	Rhombohedral	R-3c	−0.682

Table 4.

Summary of imaginary phonon frequency at each symmetric q-point in P4mm phase of LiNbO₃ and corresponding structural phase transition.

“Structure” indicates Bravais lattice of modulated structure from P4mm phase. Space group of the relaxed modulated structure and energy gain is also shown.

Figure 7 shows schematic illustration of phase transition mechanism from cubic perovskite structure to LiNbO₃ structure. On the contrary to the result of Figure 6a, R₂₅ soft mode is represented as rotation of NbO₆ polyhedra. Then, Γ₁₅ soft mode of R$\overline{\mathit{3}}$c phase leads R3c phase, which is ground state of LiNbO₃. Although the present study shows that perovskite-structured LiNbO₃ is thermodynamically unstable while its piezoelectricity is excellent, it can be possible to control phase transition behavior by dopant substitution.

Figure 7.

Phase transition path between cubic perovskite structure (Pm3m) and LiNbO₃-structure. Unit cell of LiNbO₃ is enclosed with blue lines.

7. Summary and conclusion

In this chapter, we briefly introduced sophisticated method of density functional perturbation theory. DFPT can effectively calculate the second derivative of the total energy with respect to the atomic displacement within the framework of first-principles calculation. By using DFPT method, we can predict piezoelectric constants, dielectric constants, elastic constants, and phonon dispersion relationship of any given crystal structure. Moreover, we showed our established computational technique to decompose piezoelectric constants into each atomic contribution, which enable us to gain deeper insights to understand the piezoelectricity of material. By using LiNbO₃ as a model material, we showed the predictability of DFPT for piezoelectric properties. In addition, we showed that superior piezoelectric properties are hidden in perovskite-structured LiNbO₃. Structural relationship and possible phase transition path between LiNbO₃ structure and perovskite structure were discussed and concluded that perovskite-structured LiNbO₃ is thermodynamically unstable. Further studies are expected to control relative phase stability between perovskite and LiNbO₃ structure by dopant substation and solid solution.