High-throughput ab initio calculations on dielectric constant and band gap of non-oxide dielectrics

High-k dielectrics, materials having a large band gap (Eg) and high dielectric constant (k) simultaneously, constitute critical components in microelectronic devices. Because of the inverse relationship between Eg and k, materials with large values in both properties are rare. Therefore, massive databases on Eg and k will be useful in identifying optimal high-k materials. While experimental and theoretical data on Eg and k of oxides are accumulating, corresponding information is scarce for non-oxide dielectrics with anions such as C, N, F, P, S, and Cl. To identify promising high-k dielectrics among these material groups, we screen 869 compounds of binary carbides, nitrides, sulfides, phosphides, chlorides, and fluorides, through automated ab initio calculations. Among these compounds, fluorides exhibit an Eg-k relation that is comparable to that of oxides. By further screening over ternary fluorides, we identify fluorides such as BiF3, LaF3, and BaBeF4 that could serve as useful high-k dielectrics.

growing computational speed, it becomes possible to conduct massive calculations on dielectric and electronic properties of crystals. In a previous study 20 , we carried out high-throughput DFT screening over ~1,800 oxides and identified new candidate high-k oxides such as c-BeO whose figure of merit far exceeding that of industry-standard HfO 2 . However, discussions in the above imply that non-oxide dielectrics would be also valuable in view of exploiting diverse chemistry. Recently, Petousis et al. performed high-throughput screening on ~1,000 inorganic compounds including both oxide and non-oxide compounds 21 . However, the search space in ref. 21 . was limited to stable or metastable phases with hull energies in the phase diagram less than 20 meV/atom, and relatively small primitive cells containing less than 20 atoms. However, such restrictive conditions could miss promising high-k materials. For example, the industry-standard HfO 2 thin film includes a significant portion of tetragonal or cubic phases that exhibit high k values 22 . However, these are high-temperature phases and their energies are higher than for monoclinic HfO 2 by more than 50 meV/atom 23 . In addition, c-BeO which was suggested to be a promising high-k candidate in ref. 20 is also a metastable phase with the energy of 483 meV/atom with respect to the stable wurtzite BeO. Therefore, a more extensive table of E g -k relations for non-oxide compounds is in demand. We note that the open material database such as Materials Project 24 or AFLOW 25 provides E g computed by the semilocal functional. However, the semilocal functional severely underestimates E g such that it is not appropriate in screening dielectric materials, which requires accurate band gaps. In addition, these open databases do not provide dielectric constants. (Materials Project seems to open dielectric constants only for materials studied in ref. 21 .).
In this study, by conducting high-throughput ab initio calculations on accurate E g and k, we screen high-k non-oxide materials such as carbides, nitrides, sulfides, phosphides, chloride and fluorides. Since the number of these compounds amounts to ~30,000 entries in the current the Inorganic Crystal Structure Database (ICSD) 26,27 , the computation on the whole materials is not feasible within our computational resource. Therefore, we first limit the screening to binary phases (869 structures) and compare E g -k relations depending on the anion species. We find that the E g -k relation of binary fluorides looks most promising and so extend the screening space to ternary fluorides (415 structures). Consequently, we identify candidate fluorides that are suitable for high-k dielectrics.

Results
Automation workflow. Figure 1 shows the overall workflow of the automated computations. First, from ICSD (the 2015 version), we garner structural information on ordered binary crystals that were experimentally identified and contain only one of carbon, nitrogen, fluorine, phosphorus, sulfur or chlorine atoms (non-oxide groups hereafter). We exclude compounds including 3d transition metal elements with partially occupied d orbitals (V ~ Cu) because their band gaps are usually smaller than 3 eV, and so they are not suitable for high-k applications. In addition, DFT + U methods that are necessary for 3d orbitals can significantly underestimate dielectric constants by hardening phonon modes 28 . We also omit large primitive cells that contain more than 50 atoms in the unit cell due to a sheer computational cost. (The total number of such structures is 53, and their dielectric constants are typically small.) After these pre-screening steps, we perform the structural relaxation and calculate E g and k using the in-house automation package (Automated Ab initio Modeling of Materials Property Package (AMP 2 )) 29 for high-throughput calculations. (See the Methods section for further computational details.).

Validation of automatic calculations.
In the present work, the band gap is calculated within the hybrid functional (HSE06) with structural parameters (lattice vectors and atomic coordinates) fixed to those obtained using the generalized gradient approximation (GGA) functional. To reduce the computational cost, we employ the HSE@GGA scheme in which the HSE calculation is performed on the band edge points identified by GGA. In ref. 20 , the HSE@GGA scheme was validated by comparing with experimental and other theoretical band gaps for selected oxides. Similarly, Fig. 2(a) tests this approach against various non-oxide compounds. The estimated band gaps are in good agreement with the experiment except for compounds with E g > 8 eV that show noticeable discrepancies 30 . This is caused by the fixed fraction of the exact exchange term; it is known that materials with large E g require higher fractions of the Fock term due to the weak electronic screening 31 . For a comparison purpose, we also present in Fig. 2(a) results with G 0 W 0 @HSE in which one-shot GW calculation is performed based on the HSE result 32 . It is seen that E g 's from G 0 W 0 @HSE are in very good agreement with experimental values. However, the GW method is too expensive to be used in the high-throughput screening.
For dielectric constants, we compare GGA and LDA results of selected compounds with experimental values (see Fig. 2(b)). The estimated dielectric constants are in good agreements with experiment regardless of the functional, though results show that GGA tends to give higher dielectric constant compared to LDA. The mean average deviation (MAD) is slightly lower with LDA than GGA (0.70 and 0.98, respectively). Considering that GGA tends to overestimate dielectric constants for high-k materials 20 , we employ the LDA scheme in evaluating dielectric constants. To note, we do not consider the HSE functional for evaluating k because HSE tends to underestimate the optical dielectric constant 32 . Furthermore, HSE is not superior to LDA for the static dielectric constant, as was demonstrated for TiO 2

.
High-throughput screening for binary non-oxides. For binary dielectrics, we consider 869 compounds (76 carbides, 123 nitrides, 132 fluorides, 205 sulfides, 194 phosphides and 139 chlorides). Among them, 20 carbides, 49 nitrides, 86 fluorides, 82 phosphides, 114 sulfides, and 98 chlorides are found to be insulators. The E g -k relation of each non-oxide group is presented in Fig. 3. All the numerical data are provided in Tables S1-S6 in Supplementary Information. For the comparison purpose, E g -k of oxides from ref. 20 are also plotted as gray points in Fig. 3(a). The metallic compounds are not displayed since their dielectric constants are ill defined. We also exclude unstable structures that show imaginary optical modes in the unit cell calculation.
In Fig. 3, it is seen that the inverse relationship between E g and k persist for every compound. However, the detailed distributions are distinct between the material groups. Especially, it is noticeable that fluorides have wide band gaps than other non-oxide groups. The band gap is related to the energy splitting between bonding and anti-bonding orbitals. Therefore, a large difference in electronegativity, small ionic size, and a high degree of orbital overlap contribute to strong bond strength and large band gap. Fluorine has the highest electronegativity in the Periodic Table and its ionic size is the smallest among considered anions. This results in a much broader E g range than those of other non-oxides. The maximum band gap follows the order of fluorides > oxides > chlorides > nitrides > sulfides > carbides > phosphides, which is also in line with the order of electronegativity. The correlation between E g and electronegativity was also discussed previously 33,34 . Distribution of dielectric constants. In Fig. 3, the inverse relationship between E g and k consistently appear, which puts a fundamental limitation on the existence of ideal high-k dielectrics that have large E g and k simultaneously. For detailed analysis on this, we divide the dielectric constant into electronic and ionic contributions (k el and k ion , respectively), and plot E g -k el and E g -k ion relations separately for all the data on binary compounds (see Figs. 4(a) and 4(b), respectively.) The density of data points is drawn in contours by representing each point with a Gaussian. In Fig. 4(a), a clear inverse relation is found between E g and k el . This can be rationalized by the fact that the band gap reflects the bonding-antibonding separation and so a small E g means higher electronic polarizability associated with the facile excitation into antibonding states, leading to higher k el . On the other hand, it is seen in Fig. 4(b) that the data points in E g -k ion are more scattered such that the inverse relationship between E g and k ion is weaker than for between E g and k el . The ionic dielectric constant is dictated by off-centering of the cations with respect to the anions (and vice versa) under electric fields, the degree of which depends on the bond strength. Compared to the band gap, the bond strength or phonon frequency is highly sensitive to the bond length. Therefore, various bond lengths among similar compositions result in a wide variation of k ion . Since the dielectric constant of high-k materials is mostly contributed by k ion , this weak correlation between E g and k ion increases a chance of finding new high-k materials by expanding the search space. High-throughput screening for ternary fluorides. Figure 3 indicates that band gaps of carbides, nitrides, phosphides and sulfides are distributed mostly over medium to small values and therefore, these material groups might not be appropriate for high-k applications. Both fluorides and chlorides follow similar E g -k relations but fluorides show a broader distribution. Considering the above discussion on a loose relation between E g and k ion , it might be worthwhile to extend the search space to ternary fluorides. There are 644 ternary fluorides reported in ICSD and the calculated E g -k relations for 415 ternary fluorides with finite gaps are provided in Fig. 5, together with the binary fluorides presented in Fig. 3(c). The numerical data are compiled in Table S7 in Supplementary Information. To rank the candidate materials, we assign E g ·k as the figure of merit (FOM) because E g and k are approximately proportional to the logarithm of the leakage current density 20,35,36 . Data points in Fig. 5 are colored according to FOM.
In Fig. 5, we identify candidate fluorides that merit consideration as high-k dielectrics. For the candidate materials, we confirm the dynamical stability with phonon analysis and neglect dynamically unstable structures. (See the Methods section for details on the phonon calculations.) As a reference, monoclinic and tetragonal HfO 2 (m-HfO 2 and t-HfO 2 , respectively) are also marked. It is seen that no fluorides outperform tetragonal HfO 2 (t-HfO 2 ) that is currently industry-standard high-k dielectrics. However, t-HfO 2 is a high-temperature phase and hence its stabilization at the room temperature requires strain engineering or dopants during device fabrication 37,38 . In contrast, many candidate materials in Fig. 5 are stable phases (9 stable and 3 metstable) and so synthesis would be more straightforward than HfO 2 .

Discussions
In Table 1, we enlist candidate fluorides that were identified in Fig. 5. Besides E g and k values, relative energies with respect to the most stable phase at ambient conditions are also provided. Some materials in Table 1 were also noted in ref. 21 but the dielectric constants were larger than the present results because of the functional difference (see above). Among the binary phases, BiF 3 (see Fig. 6(a)) looks promising with E g and k values close to those of t-HfO 2 . BiF 3 was used as dielectric buffer layer for surface plasmon resonance 39 . However, unlike t-HfO 2 that is metastable, BiF 3 is the stable phase and therefore we expect that the fabrication would be easier, and the film quality would be more uniform than HfO 2 . The polymorphs of LaF 3 are also intriguing as they possess large band gaps of 9~10 eV. LaF 3 was used as dielectric buffer layer 39 or UV coating 40 . There are several polymorphs in LaF 3 and their dielectric constants range over 11~17 with metastable phases showing larger values (see Table S3). The Pmmn structure with the largest k is shown in Fig. 6(b).
There are several Ge-F compounds such as Ge 3 F 8 , Ge 5 F 12 , GeF 2 , and GeF 4 . The band gap of these compounds range over 5~7.5 eV. (See Table S3) This is in contrast with GeO 2 whose band gap is only 2.8 eV. The small band gaps of GeO 2 and suboxides result in high-density defect states, degrading the interface quality in Ge transistors 41 . In this respect, the large band gap of germanium fluorides may contribute to forming stable interfaces in Ge devices by playing as the passivation layer that removes interface defects.
Among the ternary fluorides, BaBeF 4 is noticeable as the material possess a large band gap of 9.7 eV and dielectric constant of 17, surpassing that of LaF 3 . (See Fig. 6(c).) We note that tetragonal TlAlF 4 , which has the highest FOM in ref. 21 , is found to be dynamically unstable in the present calculation (both LDA and GGA) while dynamically stable monoclinic TlAlF 4 (see Fig. 6(d)) has k of 27. Nevertheless, the dielectric property of monoclinic TlAlF 4 approaches that of t-HfO 2 .
In summary, we conducted high-throughput calculations of E g and k for 449 binary non-oxides and 415 ternary fluorides. We confirmed that inverse relationships between E g and k are present in non-oxide compounds like in oxide compounds. Among the different anion groups, binary fluorides are the most promising as they show a wide distribution of E g . By further screening over ternary fluorides, we identified fluorides such as BiF 3 , LaF 3 , Figure 5. E g -k map for 86 binary and 415 ternary fluorides. Each material is color coded according to the figure of merit that is the product of E g and k. The candidate fluorides for high-k dielectrics are marked in red (stable phase) or green (metastable phase) circles. As a reference, monoclinic and tetragonal HfO 2 (m-HfO 2 and t-HfO 2 , respectively) are also noted. and BaBeF 4 that could serve as useful high-k dielectrics. We believe that the suggested fluoride compounds may contribute to resolving various issues in microelectronic devices that is caused by using only oxide dielectrics.

Methods
Computational details. The DFT calculations were performed using Vienna ab initio simulation package (VASP) 42-44 based on the projector augmented wave (PAW) pseudopotential 45,46 . For the exchange-correlation functional, we employ the generalized gradient approximation (GGA) in the form of the Perdew-Burke-Ernzerhof (PBE) 47 . Compounds bearing Zn, La, and Ce are performed with the GGA + U method 48 . The effective U parameter of 7.5 eV is used for Zn d and La f orbitals 49,50 while 4.5 eV is used for Ce f orbitals 51 . We also carry out the hybrid functional (HSE06) calculation 52 to overcome underestimation of the band gap in GGA. The k-point meshes are selected by ensuring that the total energy, stress tensor components, and all forces are converged within 10 meV/atom, 10 kbar, and 0.02 eV Å −1 , respectively. The atomic positions and lattice parameters are relaxed until the total energy, atomic forces, and stress tensors are reduced to within the same criteria.
Band gap. We adopt an automated approach in estimating the band gap as was detailed in our former high-throughput study 20 . To be brief, we first identify k points corresponding to the valence band maximum and conduction band minimum by sweeping along the lines connecting high-symmetry points 53 using GGA(+U) functional. Then, we perform the one-shot hybrid functional calculation on these band-edge points.
Dielectric constant. The density-functional perturbation theory (DFPT) 54,55 implemented in VASP is used to estimate Born effective charges, phonon frequencies, and dielectric constants. In calculating the dielectric constants, the k-point density is doubled because DFPT is sensitive to this computational parameter.
Phonon analysis. Phonon bands are calculated using the PHONOPY package 56 .

Data Availability Statement
All data sets used in this work are available from the corresponding author on reasonable request.