Data-driven discovery of high performance layered van der Waals piezoelectric NbOI2

Using high-throughput first-principles calculations to search for layered van der Waals materials with the largest piezoelectric stress coefficients, we discover NbOI2 to be the one among 2940 monolayers screened. The piezoelectric performance of NbOI2 is independent of thickness, and its electromechanical coupling factor of near unity is a hallmark of optimal interconversion between electrical and mechanical energy. Laser scanning vibrometer studies on bulk and few-layer NbOI2 crystals verify their huge piezoelectric responses, which exceed internal references such as In2Se3 and CuInP2S6. Furthermore, we provide insights into the atomic origins of anti-correlated piezoelectric and ferroelectric responses in NbOX2 (X = Cl, Br, I), based on bond covalency and structural distortions in these materials. Our discovery that NbOI2 has the largest piezoelectric stress coefficients among 2D materials calls for the development of NbOI2-based flexible nanoscale piezoelectric devices.

Dynamical charge ( ,1 * ) and 1 of each atom in NbOX2.------------21 Fig. S17 High-throughput calculation results for maximum out-of-plane sheet piezoelectric tensor elements.  Table S15 Table of    Here we note the disparity in the 22 values for As2X3 (X=S, Se) obtained in our study and those by Gao et al. 1 and attribute the disparity to the use of different exchange-correlation functionals.
In the figures below, we verify the dynamical stability of NbOI2 by making sure there are no imaginary phonon frequencies for all q in the phonon dispersion spectra as well as by ensuring that there is no bond breaking in the molecular dynamics calculation.
The phonon dispersion of NbOI2 is calculated using Quantum Espresso 7,8,9 working with a Projector Augmented Wave (PAW) approach 10 with the Perdew, Becke and Ernzerhof (PBE) Generalized Gradient Approximation (GGA) of the exchange-correlation functional 11 . We use a kinetic energy cut-off of 60 Ry for the plane wave basis set and a Monkhorst-Pack k point meshes of 12×6×1 is used. The atomic coordinates are fully relaxed using the conjugate gradient scheme until the maximum energy difference between iterations is less than 10 −14 Ry and the residual force is less than 0.0001 Ry/Bohr. confirming that the structure is dynamically stable.
Ab initio molecular dynamics of NbOI2 is performed using the VASP code 12 using the same set up as the targeted study on NbOX2. Here, we use a 4×2×1 supercell and a Monkhorst-Pack k point mesh 13 with a density of 3×3×1. We use the NVT ensemble at 298 K with the Nose-Hoover thermostat for 5ps and a timestep of 1fs.   and Compliance tensor elements S11, S12, S22, S66 in m N -1 .
Headings of the first 3 columns are the exact database dictionary keys used in this database.
Heading of the last column is in the format MP_{dictionary_key_in_materials_project}.          The in-plane cofficients of α-In2Se3 are also measured and found to be 11 = 1.7 pm/V and 22 = 6.0 pm/V.

Formalisms
In this section, refers to the three-dimensional (3D) relaxed-ion piezoelectric stress tensor, which is defined by the relation 29 where is the electric polarization, is the homogeneous strain, is the mechanical stress, and is the homogeneous electric field. = { , , }) and = {1 … 6} as in Voigt notation.
consists of two parts: clamped-ion contributions ( ) as well as ionic contributions ( ) 30 .
is a second-derivative response function tensor of energy (E) with respect to homogeneous electric field and homogeneous strain.
It is considered a clamped-ion quantity because the ionic coordinates ( ) are not relaxed when the homogeneous electric field and strain are applied, hence it represents the piezoelectric contributions from the electrons alone.
On the other hand, is defined, in an implied sum notation, as where Ω 0 is the cell volume before deformation and * is the Born effective charge (m is a composite label for atom and displacement directions, ranging from 1 to 3N). represents the ionic contribution to due to the relaxation of ionic positions after the application of strain.
is computed using density functional perturbation theory (DFPT) 30  The sheet defined in the main text is obtained by multiplying the 3D by the cell height. Table S7 is obtained from 3D by setting elements related to the z-direction to 0 and dividing the rest of the elements by cell height 33 .

The sheet defined in
Using both sheet and sheet for calculation of is equivalent to using their bulk counterparts because the scaling by height in both terms are cancelled off.
Note that when strain and electric field are simultaneously present, a more accurate formulation, as presented by Wu et al. 30 , needs to be invoked. For the sake of readability, we present the "improper" formulation here while understanding that the "proper" terms as presented by Wu et al. are used in the DFPT calculations. Also, we note that the "improper" and "proper" terms are equivalent when = {1, 2, 3}.

Supplementary Note 1: Polarization switching in NbOX2
We study the polarization switching in bulk NbOI2 microscopically via polarization versus electric field (P-E) measurements 34,35 and locally through spectroscopic PFM characterizations.
The P-E ferroelectric hysteresis loop of NbOI2 when an external electric field is applied to its polar axis is displayed in Fig. S18b. When the applied electric field is in the vicinity of the coercive field (~ 8.5 kV/cm), the polarization shows a drastic variation due to domain reversal, i.e., the field is large enough to switch domains with the unfavorable direction of polarization. hysteresis disappears when the field is applied along the nonpolar axis of NbOI2 (Fig. S18c).

Supplementary Note 2: Ferroelectric-paraelectric phase transition in NbOI2
The ferroelectric-paraelectric phase transition in NbOI2 was confirmed by temperature- This corresponds to a transition from ferroelectric phase with non-centrosymmetric C2 (No. 5) symmetry to paraelectric (PE) phase with centrosymmetric C2/m (No. 12) symmetry. SHG is highly sensitive to the inversion-symmetry breaking that accompanies a ferroelectricparaelectric phase transition; only non-centrosymmetric structures are capable of emitting SHG light. We found that the SHG signal vanishes once is exceeded and emerges again upon cooling (Fig. S20b).

Supplementary Note 4: Strain-assisted Ferroelectric Switching
As discussed in the main text, the magnitude of the applied electric field can be much reduced if one applies compressive strain to the materials. This effect is especially large for NbOI2 which has the largest piezoelectric effect.

Fig. S22 | Effect of strain on the Nb-O bond lengths in NbOX2.
The black and red lines denote the longer and shorter Nb-O bonds respectively. A centred symmetric structure is obtained when compressive strain of over -3.5%, -3.0% and -2.5% along the -direction is applied to NbOCl2, NbOBr2 and NbOI2 respectively.
The large response of Nb atoms in response to strain in the -direction also enables both the magnitude of polarization and the ferroelectric switching barriers to be controlled by strain. A modest compressive strain of -2.5% in NbOI2 results in a centred symmetric structure ( Supplementary Fig. S22), allowing the polarization direction to be set with an electric field of small magnitude. The polarization magnitude can then be enhanced by applying tensile strain. Table S15 | Table of