Data-driven methods for discovery of next-generation electrostrictive materials

All dielectrics exhibit electrostriction, i


INTRODUCTION
The development of novel materials exhibiting exceptional strain response in the presence of an external electric field is imperative to the development of tunable electromechanical devices [1][2][3][4][5][6][7][8] .Promising applications involve next-generation transducers 9 , ultra-efficient actuators 1,2,10 , sonars 11 , and low power logic devices 12 with application spaces in medicine, astronomy, and consumer electronics.Recent discoveries of 'giant' electrostrictors, for example, Gd-doped ceria 8 , and (Nb, Y) doping in bismuth oxide 13 and lanthanum molybdenum oxide 7 , soft nano-composite based on carbon nanotubes (CNT) 14 and liquid crystalline graphene 6 have opened promising applications for superior electromechanical properties, and also as a replacement of lead-based Pb(Ti 1-x Zr x )O 3 (PZT) solid solutions that are widely used in current applications.Of particular interest is the possibility to induce an effective piezoelectric behavior under bias electric field, generating extraordinarily large piezoelectric coefficients (d 33

~
pC N -1 ) 15 .These 'giant' electrostrictors have been reviewed recently 16 .There is a clear advantage in utilizing electrostrictive materials that exhibit a quadratic strain in response to an electric field.This may allow for a more significant strain response for a given applied external field, as well as take advantage of the anhysteretic nature of electrostriction for highprecision displacements.
This effect is mediated via an applied electric field ( ) and the dielectric polarization ( ) vectors as described by: ( ( where are the strain components, and are piezoelectric tensors, and and are fourth-rank tensors describing the electrostrictive response.In Equation (1), the total strain (at constant stress) is expressed in terms of an applied electric field, whereas in Equation (2) it is given as a function of the polarization.The M tensor is of particular interest for actuator applications whereas Q is used as a design criterion in sensing applications.
The search for novel electrostrictors has historically appeared bounded by the universal empirical relation connecting the electrostrictive coefficient to the ratio of the elastic compressibility over the permittivity 5 until the discovery of 'giant' electrostrictors that are several orders of magnitude greater than the empirically expected values 16 .These include 15 , lead halide perovskites ( m 2 V -2 and 1266 -1417 m 4 C - 2 ) 17 , and Ce 0.8 Gd 0.2 O 1.9 ( m 2 V -2 ) [6][7][8] .The electrostriction coefficients of these materials are three orders larger than typical perovskites, including BaTiO 3 , PbTiO 3 , SrTiO 3 , and Pb(Ti 1-x Zr x )O 3 for which is in the range 0.02 -0.07 m 4 C - 2 5 .
Using large open materials databases composed of accurately determined materials properties, the programmatic search for materials exhibiting specific properties becomes possible.In line with data-driven methods, advanced modeling approaches across length scales, and the successful insertion of the concept of materials genomics that allows for rapid development and deployment of novel materials 18,19 , we employ here a combination of data mining and first-principles modeling in the discovery of materials exhibiting exceptional electrostrictive coefficients.Our results show several non-oxide inorganic compounds with electrostrictive properties of one or two orders of magnitude higher than existing oxide ceramics and with comparable or higher electrostriction coefficients of polymers.

Assessment of Electrostrictive Properties
The tensorial notation that has been employed in Equations ( 1) and ( 2) can be reduced to the following criteria to assess the magnitude of the electrostrictive properties of a wide variety of materials (see Supplementary Information), which is also defined as the Newnham's proxy , based on the approach presented in Refs. 5,10.can be derived from as below: (3) Here is the average eigenvalue of the dielectric tensor with , the relative permittivity of the material, and the permittivity of vacuum. in Equation ( 3) is the Reuss elastic compressibility given as 20,21 : ( The Supplementary Information provides a more detailed treatment to define constitutive relations for the electrostrictive coefficients and .Equations ( 3) and (4) follow from specific electrical and mechanical boundary conditions corresponding to hydrostatic pressure/expansion and a single electric field component (uniaxial applied electric field in a short-circuit configuration) as discussed in the Supplementary Information.We note that Equation (3) describes the criterion used by Newnham et al. in their landmark paper that classified electrostriction for a range of materials over an extensive response space 5 .We add here an additional criterion for M h 16,22 .Although closely related, the quantities M h and Q h are not interchangeable as converse effects (Supplementary Information).Subtle differences in how these are defined give rise to slightly different criteria for the proxy derivation used to approximate M h and Q h .M h is readily accessible physical parameter in experiments.
In addition to the static empirical relations in Equation (3) and Equation ( 4), it is advantageous to use the differential forms of Q h and M h as described in the Supplementary Information, which is essentially the dielectric response with change in external hydrostatic pressure ( ).In the following, Q h and M h derived from our density functional theory (DFT) calculations are defined as and , respectively.Thus, It must be noted that Equations ( 6) and ( 7) capture the mechanical response as differential in pressure ) and the equations are independent of Equation (3) and Equation (4).
DFT has been used to study electrostriction properties of compounds having high structural symmetry 23,24 .Of particular interest in this regard is the study by Tanner et al. wherein the electrostrictive coefficients of cubic dielectric materials were computed using the differentials of the lattice parameters, which is similar to differentials of hydrostatic pressure for cubic crystal system 25 .We have adopted the latter method using Equation ( 6) and Equation (7).
In our computations, the trace of external stress tensor, defined as pstress, is treated as the hydrostatic pressure with the conventional notation that positive and negative pstress values mean hydrostatic expansion and compression, respectively.This approach allows for uniform evaluation of electrostrictive coefficients over all seven crystal systems independent of their crystal structure.

Data Mining
To assign proxies for selecting materials exhibiting high electrostriction, we utilize a set of relationships between the dielectric response, the elastic coefficients, the electric field-induced strain, and the electrostrictive properties of materials as described in the previous section.A similar methodology developed by Newnham et al. has shown a semi-empirical connection between the elastic compliance, the dielectric response, and the hydrostatic electrostrictive Existing databases of calculated functional properties provide easy access to data that can be used to shortlist materials with exceptional properties via a simple database query of related proxies 26,27 .The Materials Project (MP) serves as an ideal platform for novel materials discovery, given the large number of calculated material properties and optimized crystal structures contained in the database 28 .The database contains 124,515 inorganic compounds with properties calculated via medium / high accuracy DFT methods.The MP team has calculated the elastic tensors for 13,751 and the dielectric tensors for 4,892 inorganic materials so far 28 and made these publicly available 20 .In addition, there are currently 1,974 materials for which MP has calculated both the dielectric and elastic properties.This set of materials is ideal for data mining and training machine learning models to predict electrostrictive response.Starting from this set, we down-select using a set of simple property filters to find new promising candidates for electrostriction applications.We limit the materials by considering non-elemental (> 1 element) dielectrics with a bandgap > 0.5 eV and seek compounds that should exhibit relative stability as determined via the energy above the hull (< 0.1 eV).The schematic representation of the workflow is shown in Fig. (1).
These materials are then separated into centrosymmetric and non-centrosymmetric datasets based on their point group symmetries.This discretization serves the purpose of allowing us to consider materials lacking an inversion symmetry (non-centrosymmetric), which exhibit a spontaneous polarization via a piezoelectric response resulting from the dielectric and electrostrictive contributions separately from those that solely display an electrostrictive response (see Fig. 1).Taking these latter materials, it is then possible to consider the ratio of for each material and rank these according to the proxies described in Equations ( 3) and (4).To reduce the dielectric and elastic tensors to a scalar value for each data point, the inverse of the Reuss average bulk modulus was used in place of the elastic compressibility ( ) and the averaged dielectric tensor eigenvalue in place of , thus is considered a replacement for and a proxy for predicting .Alternatively, is the expression used for querying for high .These proxies are utilized throughout the rest of the study in combination with DFT calculated values for shortlisting materials which should exhibit exceptional electrostrictive strain response to an applied external electric field.

First Principles (DFT) Validation
The accuracy of DFT for materials properties that are determined at the electronic and atomic levels has been well tested; it provides a platform to bridge the gap between fundamental theory, high-throughput computational data science, and experiments [29][30][31][32][33][34][35][36][37] .Here, we perform first-principles calculations with DFT using the plane-wave pseudopotential method 38 .The structural units of the compounds are taken from the MP database.We have used the PBEsol exchange-correlation functional 39 in place of the standard General Gradient Approximation (GGA) in the MP.PBEsol is a revised Perdew-Burke-Ernzerhof GGA chosen due to the enhanced prediction of equilibrium properties of densely packed solids, such as the dielectric and elastic tensors for insulating materials compared to GGA 25 .Geometric optimization was carried out to obtain the equilibrium lattice parameters and corresponding atomic positions.The k-points grid size was set to 9 × 9 × 9, but models with atoms larger than 20 were treated with 7 × 7 × 7 grid size.The kinetic energy cutoff for the plane waves is 800 eV.The symmetry tag was switched on for the calculations so that the system did not deviate from the initialized structures.
Non-spin polarized calculations were performed as most compounds were not expected to have a magnetic ground state.Geometrical optimization for the models was carried out with the tolerance for total energy convergence set to 10 −5 eV and was performed in five cycles.The relaxed model from one cycle was then fed as the initial model for the next cycle to ensure a thorough relaxation.The calculations were performed with the Vienna ab initio simulation package (VASP) [40][41][42] .
The elastic compliance tensor was calculated by performing six finite distortions of the supercells and constructing the Hessian matrix from total energies obtained from density functional perturbation theory (DFPT) 43 .The static ion-clamped dielectric matrix was computed using DFPT via the linear Sternheimer equations 44,45 .This method bypasses the complexity of computing the dielectric function through relaxation under finite electric fields.Equations ( 6) and ( 7) rely on the change in dielectric response as a function of hydrostatic pressure 46 .The DFT computation of the dielectric function was done at 25 kbar pstress above and 25 kbar pstress below the equilibrium volume for which the pstress is zero.The mechanical perturbation of the order ± 25 kbar is well within the harmonic approximation and provides minimal deviation of the unit cell symmetry.We found slight variations in the equilibrium lattice parameters when compared to the MP.This is due to differences in the pseudopotentials that were utilized in this study in the relaxation procedures.

DISCUSSION
Applying the selection criteria developed in Section III has allowed us to shortlist a set of materials that may potentially exhibit a significant electrostrictive strain response.5]47 in Fig. , respectively).We also plot a comparison of with and with in Fig. 3(a) and Fig. 3(b), respectively.There are significant deviations from the linear correlation that follows from the Newnham semi-empirical correlation between and .The result ascertains that Newnham's semi-empirical relation, is not a universal measure that is applicable across all materials.
The most striking aspect of the top compounds is that they are halides, chlorides, bromides, and iodides.We anticipated this due to recent observations of high electrostrictive properties in fluorides such as CaF 2 , BaF 2 , SrF 2 , LiF, and NaF with electrostrictive responses around 0.36 -0.46 m 4 C -2 3-5,47 .It is interesting to note that none of these compounds made our top ten obtained either via data mining of MP or from further refinement of MP results through high-precision DFT.The materials that are discovered through our approach are mainly binary non-cubic semi-metal halides.We note that three materials appear common to the top ten and top ten datasets: GaBr 3 , WCl 6 , and SbCl 3 .These compounds have not yet been explored with respect to their electrostrictive properties.We certainly hope that our results will generate more experimental work to verify our predictions.We note that not all the compounds in Table 1 are viable materials, some are liquid at room temperature (e.g., SiCl 4 , PCl 3 , PBr 3 O), some are gases or volatile at room pressure (e.g., BCl 3 , SF 6 , or WCl 6 ).But others such as GaBr 3 with of 560.53 m 4 C -2 show significant promise.
To explain these results and to identify the origin of the giant electrostrictive properties of materials displayed in Tables 1 and 2, we carried out a detailed analysis of bonding, charge distribution, and hybridization for selected compounds.For example, we have examined GaBr within these bands, the DOS is composed of narrow peaks for GaBr 3 whereas Ga 2 O 3 exhibits a more continuous dispersion of its DOS.
The consequence of DOS differences is fundamentally connected with difference in elastic properties and thus on the compressibility , which for GaBr 3 is 22.67 GPa -1 .This value is three orders of magnitude larger than of Ga 2 O 3 (0.06 GPa -1 ).The plot of versus is shown to gain an insight into which physical factors play more important role in the derived properties.Fig. ( 5) is a four-dimensional plot with and as the absicissa and ordinate.The color coding differentiates the crystal system of the compounds, and the circle diameter scales according to the value of .Note that the values of and are independently computed using DFT.The figure is divided into four quadrants.Quadrant-II covers the region of large but small , and has large .Conversely, the quadrant-IV covers the region of small but large , and has smaller .These seem consistent with the semi-empirical relation of Equation (1).
However, we observe that quadrant-III is filled with data points showing that there are several materials with low and low that make up the chemical space.On the other hand, quadrant-I is almost empty.The result does not reflect the same trend as shown in Refs. 5,10, indicating that electrostrictive properties cannot be generalized using Newnham's proxy defined through Equations ( 3) and ( 4).The and derived from the Materials Project data are shown in Supplementary Figure 1 for all the datamined compounds and is provided in the Supplementary Information.
As a cross-check of the quality of data obtained from our study, we calculated the effective piezoelectric coefficients and from the derived and .The electrostriction coefficients are related to piezoelectric charge (d) and voltage (g) coefficients as: Trujillo et al., Discovery of Novel Electrostrictors, page 13 (8)   where is the breakdown voltage.This quantity essentially determines the maximum amount of strain that can be obtained through an applied electric field which is also true for electrostrictors.
For stress sensing, the piezoelectric charge coefficient (d) relates the applied stress ( ) to the dielectric displacement generated ( ) in short-circuit conditions whereas the piezoelectric voltage coefficient (g) relates the applied stress to the generated electric field ( ).
There is no direct way of computing .We used the breakdown voltage expression from Refs. 48,49: ( where is the maximum phonon frequency at the -point that were calculated from the linear response approach and is the DFT obtained band gap.Although the band gap is underestimated in DFT, we do not expect it to affect the results of because Equation ( 10) is a machine-learned relationship that already accounts for the underestimated band gaps.7), Eq. ( 8), and Eq (9), respectively.The unit of , , and are MV/m, m/V and, Vm/N, respectively.The data for of top-ten materials are tabulated.All calculations were carried out at 0K and therefore on materials in their solid states.

4 C 2 .
3 , one of the top candidates with = 560.53m 4 C -2 and = m 2 V -2 .Its oxide counterpart Ga 2 O 3 has significantly lower electrostriction response factor, = 0.02 m As such, Ga 2 O 3 did not make it to our lists even within the top 50 materials.The crystal structure and the atomic charge densities of GaBr 3 and Ga 2 O 3 are shown in Fig. 4 (a) and Fig. 4 (b), respectively.The atom resolved density of states (DOS) of the respective materials are plotted in Fig. 4 (c) and Fig. 4 (d).Despite each having 24 electrons per formula units in their valence band, the DOS is quite different.The bandgap observed from the DOS of the two compounds are 3.36 eV for GaBr 3 and 2.04 eV for Ga 2 O 3 , which is greater than 0.5 eV, a criterion used for our data mining.This corroborates with the charge separation observed between the cation and anion sites shown in Fig. 4 (a) and Fig. 4 (b) for GaBr 3 and Ga 2 O 3 , respectively.Owing to the large band gaps, the two materials have similar order of magnitude permittivities ( 2.86 for GaBr 3 and 9.28 for Ga 2 O 3 ).The overall energy spread of both the valence and conduction bands of GaBr 3 is smaller than for Ga 2 O 3 .In addition, Trujillo et al., Discovery of Novel Electrostrictors, page 12

Figure 1 .
Figure 1.Flow chart showing data curating and first-principles computation on datamined compounds.Data-mining Materials Project database and DFT computation workflow for search of high electrostriction materials.

Figure 3 .
Figure 3.Comparison of electrostrictive coefficients of and of obtained from proxy and DFT computation.Comparison of (a) versus and (b) versus obtained from direct application of empirical formula to the Materials Project data and our DFT computations.

Figure 4 .
Figure 4. Electronic structure of selected compounds.The crystal structure and the charge densities of GaBr 3 and Ga 2 O 3 are shown in (a) and (b).The atom projected electronic density of states for GaBr 3 and Ga 2 O 3 are shown in (c) and (d).

Figure 5 .
Figure 5. Plot of permittivity versus compressibility obtained from DFT calculations.Electrostrictive coefficient as in Equation (3) falls along the quadrants I and III with larger values appearing is quadrant-III and small values in quadrant-I.
ten high materials appear in both proxy and DFT data sets.These are: SiCl 4 , BCl 3 , PCl 3 , GaBr 3 , and WCl 6 .On the other hand, three out of the top ten materials appear in both the proxy and DFT data sets.These are SbSI, SbCl 2 , and WCl 6 .Materials in the top ten list of WCl 6 and PBr 3 are 22.83 m 4 C -2 and 20.99 m 4 C -2 , respectively, while for these compounds is significantly larger: 560.53 m 4 C -2 for WCl 6 and 393.32 m 4 C -2 for PBr 3 .Other materials such as GaCl 3 have similar values from the empirical and DFT relations (20.97 m 4 C -2 and 17.78 m 4 C -2 DFT calculations are listed in Tables1 and 2, respectively.It is encouraging to find that five out Trujillo et al., Discovery of Novel Electrostrictors, page 10 of our top

Table 3
common in the top ten list of and .The data signifies further that fluorides, chlorides,

Table 1 .
The list of top ten materials predicted from Newnham's proxy relation ( ) and our DFT calculations ( ).The unit of is m 4 /C 2 .All calculations were carried out at 0K and therefore on materials in their solid states.

Table 2 .
The list of top ten materials predicted from and our DFT calculations ( ).The unit of is m 2 /V 2 .All calculations were carried out at 0K and therefore on materials in their solid states.

Table 3 .
Values of dielectric breakdown, , piezoelectric coefficients , and computed from calculated from Eq. (