Effects of lattice parameters on piezoelectric constants in wurtzite materials: A theoretical study using first-principles and statistical-learning methods

Longitudinal piezoelectric constant (e33) values of wurtzite materials, which are listed in a structure database, are calculated and analyzed by using first-principles and statistical learning methods. It is theoretically shown that wurtzite materials with high e33 generally have small lattice constant ratios (c/a) almost independent of constituent elements, and approximately expressed as e33 ∝ c/a − (c/a)0 with ideal lattice constant ratio (c/a)0. This relation also holds for highly-piezoelectric ternary materials such as ScxAl1−xN. We conducted a search for high-piezoelectric wurtzite materials by identifying materials with smaller c/a values. It is proposed that the piezoelectricity of ZnO can be significantly enhanced by substitutions of Zn with Ca.

I n recent times, owing to their suitable mechanical and semiconducting properties, piezoelectric wurtzite semiconductors such as ZnO and GaN have received a lot of attention in the form of piezotronics and piezo-phototronics device materials. 1) The wurtzite-type piezoelectric materials, especially AlN, have another advantage that these can be used in high-temperature environments, e.g., as sensors in automobile engines, because their noncentrosymmetric structures are stable even at high temperatures. 2,3) However, the piezoelectric constants of wurtzite-type materials are generally much smaller than those of the perovskite-based materials such as Pb(Zr x Ti 1−x )O 3 , known as PZT, by few orders. It remains a challenge to explore better piezoelectric wurtzite materials; there are many reports aiming to enhance piezoelectricity by element doping or codoping into parent materials. [4][5][6][7] Among the wurtzite materials, the highest piezoelectricity has been experimentally discovered for Sc x Al 1−x N (about 25 pC=N for x ∼ 0.5). 2,3) First-principles calculations have shown the microscopic effects of Sc on the piezoelectricity in Sc x Al 1−x N by demonstrating an enhancement in piezoelectricity with increase in x. 8,9) Tasnádi et al. have revealed that the piezoelectric enhancement comes from the large responses of atomic positions to external strain and the significant softening of elastic constant. 8) However, novel materials which are superior to Sc x Al 1−x N have not been synthesized yet as there are no clear and general material-design criteria practically usable for enhancing the piezoelectricity of wurtzite materials. It has been proposed theoretically that there exists a relationship between piezoelectricity and microscopic structures, 10) but the physical mechanism behind the relation with intricate parameters is not fully understood yet. Therefore, better understanding of piezoelectricity in wurtzite materials might lead to simple design criteria making it possible to find a controllable key parameter describing wurtzite piezoelectricity with the help of first-principles and machine learning techniques which have generated a great interest in the recent times. 11,12) In this study, we calculate longitudinal piezoelectric constants (e 33 ) by first-principles method for more than a dozen binary wurtzite materials and investigate the relationships between the piezoelectric constants and several material parameters using the statistical learning methods. We show that a strong correlation exists between e 33 and the lattice parameter ratio (c=a), providing a guiding principle that higher piezoelectricity can be obtained for smaller c=a ratio. This guiding principle could be applied not only for binary materials but also for ternary materials of doped systems. Based on the above relation, we have explored highly piezoelectric wurtzite materials by uncovering small c=a materials.
The wurtzite crystal structure has three independent structure parameters, i.e., the lattice constants (a and c) and the internal atomic coordinate along the c-axis (u). In our study we calculated the parameters for all the binary wurtzite crystals listed in the structure databases. a, c, and u are fully optimized from the first principles for all the calculations. The first-principles calculations are based on the density functional theory (DFT) within the Perdew-Burke-Ernzerhof generalized gradient approximation (GGA) to the exchange and correlation functionals. 13,14) The computations were performed using the HiLAPW code, which employs the allelectron full-potential linearized augmented plane wave (FLAPW) method. [15][16][17] Longitudinal piezoelectric constant along c-axis (e 33 ) was calculated by means of the Berry phase method. [18][19][20] Statistical data analyses were performed using the statistical computing software R with the glmnet-package; this package fits the least absolute shrinkage and selection operator (LASSO) and the elastic-net regularized generalized linear model (GLM) for regression. 21,22) The GLM analysis can identify meaningful variables describing a target variable among a given set of explanatory variables. Table I shows the calculated equilibrium structure parameters (a, c, c=a, and u), piezoelectric constants (e 33 ), Born effective charges ðZ Ã 33 Þ, band gaps (E g ), and elastic constants (C 33 ) of eighteen binary wurtzite materials which are listed in the Inorganic Crystal Structure Database (ICSD). 23 26) the negative piezoelectricity has recently attracted much attention. 29) Among the calculated binary materials, AlN and ZnO are known to exhibit better piezoelectric performances.
From Table I, we can roughly observe a general trend that the value of e 33 increases with decreasing value of c=a indicating that the structural parameter could be a key index describing the piezoelectricity of wurtzite materials. Value of u also increases with decreasing c=a, showing a strong correlation between the structure parameters. The calculation result shows that e 33 does not clearly correlate with the values of the electronic structure parameter E g and energy profile parameter C 33 . Note that E g calculated using DFT-GGA usually underestimates the experimental values. To confirm the above trend, we performed calculations on a hypothetical material MgO with wurtzite structure, which can have much smaller c=a. 30) The calculated values of lattice parameters are consistent with the previous result for wurtzite MgO. 30) Calculated value of e 33 and c=a of wurtzite MgO support the relation of higher e 33 for smaller c=a. Note that the thermodynamically stable structure of MgO is cubic rock-salt type. The structure stabilities of wurtzite MgO have been theoretically studied and it has been concluded that the wurtzite phase transforms to the nonpolar hexagonal structure as an energetically more stable phase. 30,31) Figure 1(a) shows e 33 versus c=a of the calculated wurtzite materials. We can see that e 33 tends to increase with decreasing c=a, although most of the materials for which the calculations have been performed are concentrated in a narrow range of e 33 less than ±0.5 C=m 2 and c=a close to the ideal value of ðc=aÞ 0 ¼ ffiffiffiffiffiffiffi ffi 8=3 p ' 1:633. In order to know where the doped materials fall in this graph, we plotted the calculated values of e 33 for four Sc 0.5 Al 0.5 N models. 9) Interestingly, the ternary Sc 0.5 Al 0.5 N values are seen as an extension of the result of binary materials, showing that there exists a roughly linear relation between e 33 and c=a as e 33 ∝ c=a (we omit intercept terms in the main text). Figure 1(b) shows the calculated values of u versus c=a for the wurtzite materials. We clearly see that there exists a strong linear correlation as u ∝ c=a. Therefore, c=a or u can be the most important key parameters describing e 33 .
To further learn about the calculated e 33 data statistically, we performed the GLM analyses on these data in which the K-fold cross validations (here we set K = 6) were performed to obtain the best penalization parameters giving a minimum error between GLM and DFT results. We set 56 explanatory variables (descriptors) which included materials data (E g and C 33 ), structure data (a, c, volume V, c=a, and u), atomic data of constituent elements [anion's and cation's atomic numbers (N a and N c ), nominal ionic valences (Z v ), ionic radii (R a and R c ), 32) and electronegativities ( χ a and χ c ) 33) ], and many combinations of arithmetic operations among them. The target variable is set to be e 33 , and only five descriptors of c=a, ( χ a − χ c )=χ a χ c , ( χ a − χ c )=( χ a + χ c ), ( χ a − χ c ) 2 , and (N a − N c ) −1 are selected to be nonzero terms. The result shows that c=a is the most dominant descriptor, and others have negligibly tiny contributions with much smaller coefficients than that of c=a (by one or two orders). Therefore, the GLM analysis also shows the relationship of e 33 ∝ c=a.
We should note that, when e 33 =Z v is chosen to be a target variable for the GLM analysis, the minimum error between GLM and DFT results can be slightly reduced compared to   that of the target variable of e 33 . With e 33 =Z v as the target variable, the selected nonzero descriptors are c=a, (R a − R c ) −1 , ( χ a − χ c ) 2 , and (N a − N c ) −1 . In this case, c=a is the dominant descriptor and others have coefficients much smaller (by about two orders of magnitude) than that of c=a; this is roughly described as e 33 =Z v ∝ c=a. This result is quite natural considering the definition of piezoelectric constants, i.e., higher ionic valences can contribute to higher piezoelectric constants of materials. It should also be noted that the descriptors including the electronegativity differences ( χ a − χ c ) are always selected by the GLM analyses, though they have only minor contributions. It has been reported that the value of c=a of wurtzite materials strongly correlates with the square of electronegativity difference, 34) indicating that a smaller value of c=a can be obtained for wurtzite materials with larger ( χ a − χ c ). It has also been pointed out that a larger e 33 might be obtained for materials with larger ( χ a − χ c ). 3) Recent research shows that the relation between piezoelectricity and c=a holds for wurtzite materials, such as codoped AlN films. 35) To the best of our knowledge, it has been not theoretically clarified until now that e 33 correlates strongly with c=a.
Here we discuss the microscopic mechanism behind the relation of e 33 ∝ c=a [or the better description of e 33 ∝ Z v (c=a)]. By the definition of P 3 = e 33 ε 3 , the piezoelectric constants e 33 can be calculated by the first principles as where P 3 and ε 3 are the electric polarization and the external strain along the c-axis, respectively, and q is the electron charge. 9,20,25) For the calculated wurtzite materials, the first (clamped ion) term shows relatively small negative values and does not significantly depend on material. The second (internal strain) term dominantly contributes to the trend of e 33 versus c=a shown in Fig. 1(a); thus, Eq. (1) can be approximately expressed as e 33 / Z Ã 33 ð@u=@" 3 Þ. Z Ã 33 can be roughly approximated with Z v , and now we need to understand how ∂u=∂ε 3 behaves in wurtzite materials. Figure 2 shows the calculated total energies E, u, and uA = ∂u=∂(c=a) as a function of c=a for GaN, AlN, ZnO, and MgO with wurtzite structures. We calculate these values by changing c=a via optimization of u; the volume is fixed at each equilibrium value. GaN, AlN, and ZnO have the energy minimum at each c=a listed in Table I in their wurtzite phases, while the energy minimum of MgO is found to be in the nonpolar hexagonal phase with c=a ∼ 1.2 and u = 0.5; it has a flattened layered structure with an inversion symmetry [see the inset in Fig. 2(b)]. For all of the materials, we can see that u gradually increases as c=a decreases, reaching a value u = 0.5 at c=a smaller than about 1.2, almost independent of materials. uA, which is the slope of the u versus c=a curve, decreases as c=a decreases from 1.9 reaching the minimum uA values at c=a ∼ 1.2-1.3, after which it quickly diminishes to zero as c=a becomes further smaller. In the c=a range higher than about 1.5 [the c=a range shown in Fig. 1(a)], uA is roughly proportional to c=a. Therefore, in this c=a range, we can approximate as ∂u=∂ε 3 ∝ c=a, indicating that Eq. (1) can be approximated as This rough analysis explains our result of the GLM computation. Figure 2(c) implies that the highest value of e 33 can be realized for materials with c=a ∼ 1.2-1.4 (near the minimum value of uA depending on materials) if materials can exist stably in the wurtzite phases. It is also expected that piezoelectricity rapidly diminishes and disappears for materials with smaller c=a values which are less than the minimum uA points, because structures in such cases transform to the nonpolar hexagonal phases. Near the phase boundaries between the polar wurtzite and the nonpolar hexagonal structures, e 33 would exhibit significant nonlinear changes as Fig. 2(c), thus giant piezoelectricity could be obtained beyond the relation of e 33 ∝ (c=a). The result shows that the highest piezoelectricity in wurtzite materials would be generally obtained at c=a near the phase boundary between the polar wurtzite and the nonpolar hexagonal structures, similar to the Sc x Al 1−x N case. 8) Based on the insight above, we have tried to search for novel highly-piezoelectric wurtzite materials by finding materials possessing smaller c=a values. The most practically useful way for enhancing e 33 is element-doping into stable wurtzite materials, as previously conducted. For our study, we chose ZnO as the base material. There have been many reports on the effects of element doping into ZnO, and it has been reported that c=a of ZnO is significantly reduced by Ca doping, 37) compared to other elements. A 16-atom super cell model of ZnO was used to calculate the value of e 33 of the doped systems by substituting Zn with alkaline earth metals. We considered all the possible combinations of substitution sites and the results were averaged over the models in the same manner as in the previous calculations of Sc x Al 1−x N. 9)   Fig. 3(b)]. There is another practical reason for piezoelectricity reduction at a higher x; non piezoelectric rock-salt structures become more stable compared to the wurtzite phases for higher values of x, 36,37) e.g., for x higher than about 0.3 in Ca x Zn 1−x O. 37) Therefore, if the wurtzite structure is retained, Ca x Zn 1−x O with x ∼ 0.5 could give the best piezoelectric performance among the materials for which the calculations were performed. Though the calculated value of e 33 of Ca x Zn 1−x O is smaller than that of Sc x Al 1−x N, we expect that Ca x Zn 1−x O has an advantage concerning materials costs and earth abundances of constituent elements.
In conclusion, we have proposed a possible strategy for maximizing piezoelectricity of wurtzite materials. From the first principles calculations with the help of the statistical learning methods, we showed that e 33 linearly correlates with c=a (in the ranges higher than about 1.  Values are averaged over calculated models as described in the main text.