Large Piezoelectric Response of van der Waals Layered Solids
The bulk piezoelectric response, as measured by the piezoelectric modulus tensor (d), is determined by a combination of charge redistribution due to strain and the amount of strain produced by the application of stress (stiffness). Motivated by the notion that less stiff materials could exhibit large piezoelectric responses, herein we investigate the piezoelectric modulus of van der Waals-bonded quasi-2D ionic compounds using first-principles calculations. From a pool of 869 known binary and ternary quasi-2D materials, we have identified 135 non-centrosymmetric crystals of which 48 systems are found to have d components larger than the longitudinal piezoelectric modulus of AlN (a common piezoelectric for resonators), and three systems with the response greater than that of PbTiO, which is among the materials with largest known piezoelectric modulus. None of the identified materials have previously been considered for piezoelectric applications. Furthermore, we find that large d components always couple to the deformations (shearing or axial) of van der Waals “gaps” between the layers and are indeed enabled by the weak intra-layer interactions.
Coupling between the mechanical degrees of freedom and the electric polarization in a solid, a hallmark of piezoelectric materials, has found use in applications that range from sensors, resonators, motors, and actuators, to high-resolution ultrasound devices and miniature filters for cellular communications.jaffe2012a; uchino1996a; liu2018double Interestingly, only about 10 piezoelectric materials, including SiO (quartz), LiTaO, LiNbO, PZT (lead zirconate titanate)-based, BaTiO-based, (K,Na)NbO-based, BiTiO-based, AlN, and ZnO, are technologically relevant and cover virtually all of these applications.ieeeweigel; wu2016piezotronics; uchino1998piezoelectric Expanding the pool of known piezoelectric materials would help broaden the range of applications and allow earth-abundant and non-toxicsaito2004lead replacements for materials that are presently in use. Furthermore, it would offer more cost- and performance-effective choices beyond a relatively limited set of materials that are at present used for their piezoelectric properties.
The need for new piezoelectrics was recognized before and was the main motivation behind recent computational efforts in high-throughput screening of inorganic materials for piezoelectric performance. These include the work on perovskite alloys,armiento_PRB:2014; armiento2011screening and creation of a database of piezoelectric properties of compounds.de2015database Materials with reduced dimensionality, such as 2D mono- and multi-layers have also been investigated recently.acsnano2017; blonsky2015ab; Li2015; yin2017giant; duerloo2012intrinsic; hinchet2018piezoelectric In virtually all of these works, the quantity of interest was the piezoelectric coefficient tensor , which relates the polarization to strain ,nye1985physical; eijdij (in Voigt notation,voigt2014lehrbuch ) and is typically obtained from first-principles calculations.
In this work, we also apply first-principles calculations to screen for candidate piezoelectrics, but instead of the piezoelectric coefficient tensor we consider the piezoelectric modulus tensor , which relates the induced polarization to the applied stress, .nye1985physical; eijdij The advantage of using is in that it represents an important figure of merit in a wide range of technological applicationsjaffe2012a; uchino1996a and is a more commonly measured piezoelectric property, especially for materials in the thin film form. Because the induced polarization depends on the applied stress through a combination of the charge redistribution due to strain and the amount of strain that is produced by the applied stress, the piezoelectric modulus tensor depends on both the piezoelectric coefficient tensor and the elastic tensor as:
From equation (1), we can infer that large piezoelectric modulus can be expected in materials with large large , but also in systems with low stiffness . Softer (less stiff) materials can exhibit large piezoelectric response, as measured by , compared to stiffer materials with similar . This is illustrated in Figure 2, where we notice less stiff materials overall exhibit larger piezoelectric responses in the corresponding direction.
For this reason, we concentrate our investigation on layered (quasi-2D) van der Waals bonded solids, which can be expected to have relatively low stiffness in the out-of-plane direction due to the presence of weak van der Waals (vdW) interactions between the layers. The questions we are addressing are: (a) whether materials that belong to this class and exhibit strong piezoelectric response () can be found, and (b) if this is true, what is the role of vdW interactions in enabling the strong response. To our knowledge, these questions have not been previously addressed in a systematic way.
Our results confirm the expectations. The search has revealed a number of quasi-2D materials with relatively large components. Out of 869 considered binary and ternary layered VdW systems, we have identified 135 non-centrosymmetric crystals. Out of those we find more than one third (48 compounds) exhibit piezoelectric moduli greater than that of AlN (=5.5 pC/N), a commonly used piezoelectric material in resonator applications. In addition, we found three layered systems with piezoelectric moduli even larger than that of PbTiO (=119 pC/N ), another established piezoelectric material known for its very large response. It is important to note that none of these vdW systems have been considered previously for piezoelectric applications.
After performing a thorough analysis of the coupling between various stress and components, we find that in all of these systems, large piezoelectric response is always coupled to the stress components that imply deformations (axial or shear) of the van der Waals “gaps” between the layers. This is consistent with the fact that the softest elastic constants are related to these deformations. Ultimately, our results point to the layered vdW systems as a rich chemical space for finding new piezoelectric materials, and introduce elastic properties as additional design criteria for finding materials with large piezoelectric modulus.
Ii Computational Methodology
In this section, we discuss the details of the computational methodology used in our calculations, broadly divided into three main subsections. The first subsection describes the procedure for identifying quasi-2D structures from the Inorganic Crystal Structure Database (ICSD).bergerhoff1983inorganic; belsky2002new Next, we describe the drawbacks of GGA exchange correlation functionals for predicting the properties of layered materials and our approach to overcome this issue. Then, we provide detailed descriptions of our workflow for evaluating the piezoelectric modulus tensors (d) of quasi-2D solids.
ii.1 Automated Identification of Quasi-2D Materials
An essential component of this work is the identification of layered (quasi-2D) vdW materials from ICSD. To accomplish this, we extend the applications of our previously-developed proceduregorai2016computational for automated identification of the binary quasi-2D materials from the ICSD to include ternary chemistries. Similar algorithms have been developed and used for the purpose of broad identification of vdW bonded layered systems by others, including the work of Ashton et al.,ashton2017topology Mounet et al.,mounet2016novel and Cheon et al.cheon2017data
Our procedure relies on a slab cutting routine and bond counting. In the first step, we cut out stoichiometric slabs of a certain thickness for all symmetry inequivalent sets of Miller indices within a certain range (). Next, for each slab we find the terminations of its surfaces that minimize the number of broken bonds by translating surface atoms from one side of the slab to the other using appropriate lattice vectors. We then count the (minimal) number of broken bonds, i.e., the under-coordination, of the surface atoms. The condition of quasi low-dimensional crystals then implies the existence of directions for which the corresponding slabs do not have any under-coordinated atom relative to their bulk coordination in the first shell. If there is exactly one such , the material is a layered material with relatively large spatial gaps separating individual layers. If the number of directions is larger than one, then the corresponding systems are of lower dimensionality, quasi-1D for two such directions and molecular crystals for larger than two. If there are no such directions, the structure is then a connected 3D structure without large spatial gaps.
In Ref. gorai2016computational we demonstrated the success of our algorithm in searching complex quasi-2D materials, including those with layer stacking in oblique directions, materials with corrugated, accordion-like layers, and those with individual layers composed of multiple atomic layers. Using the described automated algorithm, we have in this work considered 3500 binary and 8000 additional ternary compounds from the ICSD and classified them into layered (quasi-2D) and not layered materials. We restricted our search to stoichiometric and ordered systems that do not contain rare earth elements and have 50 of less atoms in the unit cell. From our calculations, we have identified 426 binary and 443 ternary quasi-2D compounds. A full list of these materials can be found in supplementary information with their ICSD id.
ii.2 Calculating Piezoelectric Properties of Quasi-2D Materials
In quasi-2D structures considered in this work, the individual layers are held together by relatively weak vdW interactions. The standard exchange-correlation functionals typically employed in density functional theory (DFT) calculations, including the calculations of elastic and piezoelectric properties, are known to fail to describe the vdW interactions. This is evident from relatively large errors in the out-of-plane lattice constants and the associated elastic properties.lebegue2013two To overcome this issue, we employed a vdW-corrected functional (optB86) as implemented in VASP (Vienna Ab-initio Simulation Package) codekresse1993ab; kresse1996efficient to calculate the lattice parameters, elastic, and piezoelectric properties of quasi-2D materials.PhysRevB.83.195131; klimevs2009chemical To evaluate the piezoelectric coefficient tensors, we utilize the VASP implementation of the density functional perturbation theory (DFPT)baroni2001phonons; baroni1987green; gonze1995adiabatic calculations. A relatively large plane wave cutoff energy of 540 eV is used for structural relaxation, calculation of elastic tensors, and piezoelectric coefficient tensors. A dense k-point grid, defined by , where is number of atoms in the primitive cell and is the number of k-points, is employed. In all our calculations, a very high tolerance of 10 eV for energy convergence is used, which is an important consideration for conducting DFPT calculations.de2015database For calculation of elastic tensors we use a finite difference method. Here, the full elastic tensor is calculated by conducting six finite distortions of the lattice and obtaining elastic constants () from the stress-strain relationship.le2002symmetry; PRB2010elastic
The importance of incorporating vdW-corrections is illustrated in Figure 3, where we notice significant improvement in predicting elastic constants, particularly and , with vdW-corrected functional (optB86) PhysRevB.83.195131; klimevs2009chemical compared to standard GGA-PBE functional.perdew1996generalized The GGA-PBE predicted elastic constants are sourced from Ref.jain2013commentary. A more detailed analysis of the data presented in Figure 3 reveals that the GGA-PBE is still better in predicting in-plane elastic coefficients and , but the error in reproducing and is a factor of 10 or larger. This is due the fact that these particular two elastic constants are directly related to the deformations of the spatial gaps between the layers and the failure of GGA-PBE in reproducing relatively weak vdW interactions. The comparison of predicted properties with GGA-PBE and vdW-corrected functional is limited only to the elastic constants (Figure 3); the lack of experimental data on piezoelectric properties of quasi-2D materials prevented us from making similar comparisons for predicted piezoelectric properties. In the supplemental information Table S1, we have provided the experimental details (e.g. measurement techniques, measuring temperatures, etc.) for each compound shown in Figure 3. The comparison of calculated and experimental values for the piezoelectric modulus of a few commercially important piezoelectric materials (including AlN and PbTiO) are shown supplementary information. Predicted piezoelectric moduli are found to be in good agreement with experimental values.
ii.3 Workflow for identifying quasi-2D piezoelectrics
A complete workflow we developed for identifying promising quasi-2D piezoelectric materials is illustrated in Figure 4. The binary and ternary crystal structures from the ICSD database are first screened using the automated algorithm for identifying quasi-2D structures. Then, we filter out all centrosymmetric structures based on the space group assigned in ICSD. Out of 11500 binary and ternary materials we find 869 layered systems, out of which 135 are identified as having non-centrosymmetric structures. Next, the non-centrosymmetric, structures are relaxed using the previously described first-principles calculations employing a vdW-corrected functional. As the piezoelectric materials need to have sizable band gaps for their properties to not be screened by the existence of free charge carriers we next employ a band-gap filter. As suggested in the previous works,gorai2016computational; marom2009describing for electronic structure calculation, we perform self-consistent GGA-PBE calculations on the vdW-relaxed structures using dense k-point grids. Because of the known band gap error in DFT calculations we use a relatively generous band gap cutoff of 0.1 eV. Fifty additional materials with their band gap smaller than 0.1 eV are discarded as a result. Finite difference calculations are performed with vdW-corrected functional to obtain the elastic constants (C). Five materials with elastic tensors with negative eigenvalues are also discarded. According to Born stability criteria,born1940stability the elastically stable materials always have positive eigenvalues of stiffness matrix. This means that an elastically stable materials always have positive elastic energy for arbitrary homogeneous deformation by an infinitesimal strain.mouhat2014necessary With the remaining 80 candidates, we proceed our calculations by performing Density Functional Perturbation Theory (DFPT) baroni2001phonons; baroni1987green; gonze1995adiabatic calculations to assess their piezoelectric coefficient tensor, e. Out of the remaining 80 candidates, 38 materials contain transition metal elements. For these, we perform a limited search of the magnetic ground state by enumerating all magnetic configurations in the unit cell and calculating their total energies. The elastic and piezo calculations are then performed only on the lowest energy spin state. This step is necessary because of the strong dependence of the electronic structure and other properties on spin configuration discussed in more details in Ref. gorai2016thermoelectricity. The automated DFT calculations including initial file generation, calculating properties, data extraction and data handling are performed with the help of PyLada,pylada a Python framework for high-throughput first-principles calculations.
Iii Results and Discussions
iii.1 Promising Quasi-2D Piezoelectric Materials
Based on the calculated piezoelectric modulus tensor, a number of candidate materials with relatively large components have emerged. They are shown in Figure 5 with the full list together with the corresponding , and values provided in the supplementary information. In addition, the top 20 most promising systems, based on the largest component, are listed in Table 1.
The piezoelectric modulus tensor (d) is a third rank tensor, and any isotropic averaging scheme will yield zero.li2000effective To rank these materials based on the merit of their piezoelectric response, we define as the largest element of the absolute matrix. Then the is plotted against the space group number in Figure 5. The two reference lines have been drawn for categorizing these candidates – one representing calculated of PbTiO and the other representing calculated of AlN (see supplementary information for the benchmark agains experiments).
Please note that PbTiO is also a ferroelectric and here we are only using the value for its piezoelectric response. More precisely, the PbTiO reference line corresponds to the bulk piezoelectric modulus corresponding to the single-crystal single-domain samples.
The materials shown in the left panel of Figure 5 can broadly be divided in three categories. The first category is comprised of quasi-2D compounds with larger than the longitudinal piezoelectric modulus of PbTiO (=119 pC/N)jaffe2012a – the key end member of most commercial high-strain piezoelectrics. We found three materials (InTe, PbS, and GeTe) in this category. Among them PbS is not in its ground state rocksalt phase, but in the hypothetical distorted NiAs structure which has found its way into the ICSD.PhysRevB.84.045206 The other two compounds have previously been experimentally synthesized,chattopadhyay1987neutron; sutherland1976indium but their piezoelectric moduli have not been reported so far.
In the second category we group all compounds which have larger than the longitudinal piezoelectric modulus of AlN (=5.5 pC/N)dubois1999properties and lower than the longitudinal piezoelectric modulus of PbTiO. The majority (48 compounds) of the piezoelectric candidates from our study fall in this category revealing that overall the vdW bonded quasi-2D systems indeed exhibit a propensity toward large piezo-response. This group is composed from oxides and other chalcogenides such as CuVO, SnO, BiInO, MoVO, SnS, InSe, CsTe and other (total of 30); halides such as KSnF, AgI, MgCl, and PbI (total of 4); pnictides such as NaSnN, NaSnP, KSnAs, and NaN (total of 4). Also, a number of materials in this group (10) are the mixed anion systems, e.g., NaIO, SbFCl. Not surprisingly, the more ionic systems like oxides, halides and nitrides are more frequently found closer to the top of the range. Finally, two compounds CsTe and BiWO are found in calculations to be dynamically unstable, but both have been experimentally synthesized (likely high-temperature phases).chuntonov1982synthesis; wolfe1969crystal
The last group is composed of materials with the piezoelectric response lower than the longitudinal piezoelectric modulus of AlN. Though these candidates exhibit low piezoelectric response, they could still be useful as the calculated moduli are comparable with that of quartz (=2.27 pC/N).bottom1970measurement We found a total of 12 compounds which fall in this category. Examples include: WS, RhF, ZrCl, and GaInS. The moduli d, e and C with other informations such as band gap, space group of these compounds are provided in the supplementary information.
We also observe that the piezoelectric compounds in the quasi-2D family of solids are clustered mainly in three specific space groups, i.e., space group no. 36 (Cmc2), 160 (R3m), and 186 (P6mc). This is mainly a reflection of the population bias, as these are the three most frequently occurring non-centrosymmetric space groups in the quasi-2D family of crystals.
Table 1 shows that in 9 of the high-response quasi-2D piezoelectric compounds, the component appears as . Note that because of the freedom in choosing the in-plane axes, and are virtually indistinguishable (see the discussion section). This component corresponds to the thickness shearing deformation where the material shears like a deck of cards in the in-plane direction, with no change in the other dimension. Materials with large can be used in a variety of applications including: sensors, actuators, accelerometer, material testing structural health monitoring, non-destructive testing (NDT), and non-destructive evaluation (NDE).boivin2017torsional The components of appearing as usually coincide with the components of appearing as . The distribution of different components of (and ) appearing as (and ) are discussed in more detail in the next section.
Another quantity that could influence the piezoelectric response is the band gap of the material. In this work the band gaps are calculated at the DFT level and are also shown in Table 1 and supplementary information. Only about 1/3 of the studied materials are found to have DFT band gaps below 1 eV. Given the well known underestimation of band gaps in DFT based methods we do not think that materials with DFT band gaps larger than 1 eV would suffer from problems related to the existence of free charge carriers due to thermal fluctuations. However, for those with smaller gaps thermal fluctuation may cause sufficient number of free charge carriers, which may lower the polarization upon straining these materials despite having large piezoelectric moduli. Of course, provided that the real band gap is sufficiently close to the DFT one. For these materials, a more accurate assessment of the electronic structure might be needed before they are considered for applications.
In relation to the chemical composition and toxicity it is also important to not that currently, the most widely used piezoelectric material is lead zirconate titanate (PbZrTiO or PZT).shrout2007lead However, PZT causes significant environmental problems because of its high lead content.shrout2007lead Hence, significant efforts have been made to develop lead-free piezoelectric materials. shrout2007lead; takenaka2005current; baettig2005theoretical In our work, we have identified 44 candidates that do not contain any toxic elements including Pb. Out of 44 candidates, 33 of them have their piezoelectric modulus larger than AlN.
iii.2 Role of van der Waals Interactions
In order to understand the role of van der Waals interactions on the piezoelectric response of quasi-2D materials, we analyze the relationships of and to the corresponding strain and stress components, respectively. A histogram showing the number of compounds with a given component appearing as is shown in Figure 6 (a). We observe that the most frequent is . This indicates that in the majority of quasi-2D materials the largest piezoelectric response, as measured by the , is along the layer stacking direction. The reason for this behavior is that the relatively weak vdW interactions allow large charge redistribution in the layering direction upon straining the system.
On the other hand, the piezoelectric modulus tensor d relates relates stress to polarization and combines two types of effects: (1) amount of strain due to application of stress, and (2) amount of charge redistribution (polarization) due to resultant strain produced by the applied stress. A similar histogram of the components appearing as is shown in Figure 6(b). We divide the components into three groups depending on the deformation types (stress components) and the polarization direction. The schematics of the polarization directions and the associated stress components is shown in Figure 6(c) together with the components connecting the two. Every component of represents a separate piezoelectric operating mode. Schematics of all possible piezoelectric modes are provided in supplementary information.
Group I: Applied stress deforms van der Waals bonds and the measured polarization coincides with the direction of the deformations. The longitudinal mode () and shear modes ( and ) fall in this class. As evident from the histogram in Figure 6(b) these are the most frequently appearing components. In our considered materials, the modes and are indistinguishable because of the arbitrariness of the choice of axes ‘1’ and ‘2’, while the axis ‘3’ is fixed by layer stacking directions. In both of these modes, the same stress component () is responsible for the deformation, which implies shearing of the van der Waals gaps. On the other hand, in the , the applied stress axially deforms (stretches or compresses) van der Waals gaps. Hence, the large piezoelectric responses are achieved by deforming (axial or shear) the relatively soft van der Waals bonds. The bar-heights of and in the histogram of in Figure 6(b) are larger compared to . This is because the shearing resistance values () of quasi-2D materials are lower compared to their axial resistance values () (refer to Figure 3(b)).
Group II: Applied stress deforms van der Waals bonds, but the measuring polarization directions are different from their deformation directions. The face shear modes ( and ) and the thickness-extension modes ( and ) fall in this class. The schematics of and are provided in the supplementary information. Usually, the direction of polarization is facilitated by the direction of deformation. In such modes, the deformation directions are different from their measured polarization directions. This is the reason behind their low occurrence in the histogram of though the van der Waals bond are deformed by the applied stress.
Group III: The van der Waals bond does not deform by the applied stress in these modes. This is why these modes do not appear frequently in the histogram of . This includes length or width extension modes ( or ) and shearing modes of type and .
The above discussion implies that the large piezoelectric response is always accompanied and caused by the stress that deforms the “soft” van der Waals gaps either through stretching, compression, or shearing. In addition, the analysis of the maximal piezoelectric moduli and the associated stress components provided guidance to experimentalists of how thin films should be grown to utilize large . It also describes what kind of mechanical actuation is necessary to achieve large piezoelectric response. This also helps in the design of new devices to take advantage of large . Finally, these findings can open up a wide variety of devices based on their operation modes or based on new materials for non-conventional modes. For example, materials with large face-shear () mode response will be attractive for for torsional applications, like novel gyroscopic sensors or high-precision torsional MEMS actuators.boivin2017torsional
iii.3 Axial Piezoelectric and Elastic Anisotropy
In addition to revealing new candidate piezoelectric materials and explaining the origin of strong piezoelectric response, we have also investigated the anisotropy in axial piezoelectric response in relation to the axial elastic anisotropy. We analyze how the response in the out-of-plane direction compares to the in-plane responses. The in-plane and out-of-plane directions in quasi-2D materials are trivial to define and are illustrated in Figure 1. We define the axial anisotropy in both elastic and piezoelectric responses by the ratio of the out-of-plane component to the in-plane component. Here, the in-plane component is defined by the arithmetic mean of the ‘11’ and ‘22’-components (invariant to the choice of the in-plane axes), whereas the out-of-plane response is solely defined by the ‘33’-component. Hence, the axial anisotropy of d, e, and C can be expressed as , , and respectively. If the axial anisotropy equals or nearly equals 1 then the material is considered to be isotropic in the corresponding quantity responses with respect to the in-plane and out-of-plane directions. If the axial anisotropy is greater (or lower) than 1 the response of a material to axial deformation is dominated in out-of-plane (in-plane) direction.
All possible correlations among the axial anisotropy of d, e, and C have been investigated. The correlation between axial anisotropy of d with respect to the axial anisotropy of e and C are shown in Figure 7 (a) and (b) respectively. From the comparative studies between Figure 7 (a) and (b), we see that the axial anisotropy of d is mainly dictated by the axial anisotropy of e not by the axial anisotropy of C. The plot between axial anisotropy of C and e is provided in the supplementary information.
From Figure 7 (a) and (b), we observe that most quasi-2D materials have axial piezoelectric anisotropy parameter (both in d and e) 1 and axial elastic anisotropy parameter 1. This implies that the quasi-2D piezoelectric materials are dominant in out-of-plane piezoelectric responses but elastically they are dominant in-plane directions. This result corroborates the correlation between elastic softness and large piezoelectric response, i.e., the large piezoelectric responses are observed in elastically softer directions. We found quasi-2D piezoelectric compounds such as CuVO and LiSbO are dominated by in-plane piezoelectric response. Compounds such as BiWO, AlZnS, GeZnO, and GeTe are nearly isotropic in their axial elastic responses but highly anisotropic in axial piezoelectric responses.
In conclusion, we performed a large-scale computational (first-principles) assessment of the bulk piezoelectric properties of layered (quasi-2D), vdW bonded materials. In our study we concentrate on the piezoelectric modulus as the measure of the piezoelectric response, which relates mechanical stress and electric polarization and depends on a combination of charge redistribution due to strain and the amount of strain produced by the stress. Overall, out of 135 non-centrosymmetric quasi-2D binary and ternary structures from ICSD we have discovered 51 materials with piezoelectric response larger than that of AlN, a well-known piezoelectric materials used in applications. Out of these 51 systems, we find three with the piezoelectric modulus even larger than that of PbTiO that has the piezoelectric modulus among the largest known. More importantly, 33 out of the 51 layered compounds do not contain any toxic elements including Pb. Our results also reveal that the large piezoelectric modulus in vdW systems is directly enabled by the vdW interactions between layers as in majority of compounds the large components of the piezoelectric modulus tensor couple to the stress components that imply deformations (both shear and axial) of the “soft” vdW bonds between layers. Our results suggest that quasi-2D layered materials are a rich structural space for discovering new piezoelectric materials.
Acknowledgement The authors gratefully acknowledge the support of the National Science Foundation through Grant No. DMREF-1534503. The calculations were performed using the high-performance computing facilities at National Renewable Energy Laboratory (NREL) and at Colorado School of Mines (the Golden Energy Computing Organization) .