High‐Throughput Electronic Structures and Ferroelectric Interfaces of HfO2 by GGA+U(d,p) Calculations

The electronic structure, vacancy symmetry, defect levels, ferroelectric phases, and interface properties of HfO2 are studied using a GGA + U(d,p) approach, a simplified version of the ACBN0 method. Introducing an on‐site Coulomb interaction to both Hf 5d orbitals and O 2p orbitals reproduces the experimental bandgap, gives band energies similar to those of hybrid functionals, gives the correct symmetry for the oxygen vacancy, and describes the Schottky barriers at the metallic contacts like TiN correctly. The energetics of phase energies and strain arising from different ferroelectric–electrode interfaces are tested. The GGA + U(d,p) approach is a useful tool to study various HfO2 configurations by rapid ab initio molecular dynamics calculations.


Introduction
The demonstration of ferroelectricity (FE) in HfO 2 -based films has triggered intense research in next-generation electronic devices such as ferroelectric memories and negative capacitance field-effect transistors. [1][2][3][4] Modeling such HfO 2 films with grain boundaries at the atomic scale requires the ability to carry out various ab initio molecular dynamics (MD) calculations to understand their properties, device models, and reliability. The calculations should be able to handle large supercells, both normal to and along the HfO 2 -electrode interface, and have the correct bandgap, phase energies, and description of gap states due to grain boundaries and oxygen vacancies.
It is well known that density functional theory (DFT) with local or semilocal exchange-correlation functional underestimates the bandgap of semiconductors and insulators and can also cause problems for defects and domain-wall states. There are many methods to correct this shortcoming. At one extreme, costly methods like the GW are unable to handle large supercells. On the other hand, hybrid functional calculations, which mix a fraction of nonlocal Hartree-Fock (HF) exchange like the HeydÀScuseriaÀ Ernzerhof (HSE) [5] or screened-exchange (SX) [6] methods, correct the bandgap error but at a two-orders-higher computation cost than simple DFT functionals. Nevertheless, the added fraction of HF exchange in hybrid functionals can cause a divergence at the Fermi energy with metal electrodes, such as for Schottky barriers or metallic grain boundaries. [7] It is therefore interesting to study recent pseudohybrid functional or Agapito Curtarolo Buongiorno Nardelli ab initio DFT functional (ACBN0) methods. [8][9][10] These are essentially a type of DFT þ U with speeds of traditional DFT but with an accuracy more like hybrid functionals in terms of bandgaps, localization, and defect levels and do not suffer from a divergence at the Fermi energy.
The DFT þ U method first arose from the Anderson model for localized magnetic impurities in free-electron metals, where a Hubbard U potential was added to open-shell magnetic impurity d orbitals. [11] It was widely used for transition metals and their oxides. It was then used empirically to widen the bandgap of semiconducting oxides like ZnO, adding a U potential to closed-shell Zn 3d states, [12] and was extended further to systems like TiO 2 and Cu 2 O. [13,14] The method was widely applied but often needs unphysically large U values. It was then noted that adding a U potential to both metal d states and filled oxygen 2p states can open the gap correctly without needing such large U values. [15][16][17] Recently, the pseudohybrid or ACBN0 or GGA þ U þ V methods were introduced. [8][9][10] Here, the U term denotes on-site interactions, whereas the V term is for twocentered interactions. These studies give the methods a more firm basis and explain the effects of the various parameters.
Here we give a detailed analysis of the GGA þ U(d,p) method applied to HfO 2 to discuss its band structures, defect levels, electron affinities, and the phase energies. The calculations are carried out using the plane-wave code CASTEP. [18] Ultrasoft pseudopotentials are used for oxygen, whereas the Hf pseudopotential is generated by the on-the-fly (OTF) scheme provided with CASTEP, [19] which better describes the Hf electronic properties and needs a smaller cutoff energy of only 340 eV. The total energy is converged to 10 À6 eV/atom with a residual force of 0.02 eV Å À1 . A 7 Â 7 Â 7 k-point mesh is used for cubic HfO 2 primitive cell and a 5 Â 5 Â 5 k-point mesh is adopted for the other three phases of HfO 2 primitive cells. Different U value combinations are applied on Hf d and O p orbitals and are denoted as U ¼ (m,n), where (m,n) is U ¼ m eV on Hf d and U ¼ n eV on O p orbitals, respectively. For each (m,n) combination, the lattice constants are correspondingly relaxed.

Method
The parameters are selected by minimizing the error between the calculated GGA þ U band energy and its calculated SX value for the Γ 15 , Γ 1 , Γ 12 , and X 1 conduction band (CB) states, referred to as the Γ 15 valence band (VB) maximum, plus the O 2p VB width at X. Figure 1 shows the band structures of cubic-phase HfO 2 calculated with the Perdew-Burke-Ernzerhof (PBE) functional, the SX functional, and the GGA þ U(d,p) method. The PBE calculation clearly underestimates the bandgap to only 3.67 eV, as shown in Figure 1a. The SX band structure in Figure 1b reproduces our previous result [6] and is close to the bands found by the GW method, particularly the conduction bands such as the Γ 15 energy level. [20] The GGA þ U(d,p) band structure in Figure 1b is much closer to the SX band compared with the PBE functional, either the calculated gap value or the Γ 15 energy level. The detailed band energy comparison is shown in Table 1 and Figure 2.
The success of using U(d) terms in ZnO suggests applying U(d) to the empty Hf 5d states. However, sizeable U(d) terms can reduce the Γ 12 ÀΓ 15 CB crystal field splitting and eventually invert it. This can be a serious error as the phase energies are related to the Γ 12 ÀΓ 15 crystal field. In contrast, applying a U(p) term to O 2p levels mostly shifts the VB energies downward and leaves the CB relatively unchanged. The U(p) also widens the VB to near the correct value, but U(d) wrongly narrows it.
The error minimization also includes the phase energies for the cubic (c-), tetragonal (t-), monoclinic (m-), and orthorhombic (o-, with space group Pca2 1 ) phases in the GGA. These are expected to be related to the CB energies in that these reflect the crystal field splitting. Overall, the U(6,6) value gives the lowest mean error and is chosen.
We note that the HfO 2 phase energies among the different symmetries are in the correct order by the potentials in the VASP code, which, however, has a rather large cut-off energy of 500 eV, whereas the order for the norm-conserving potentials in CASTEP is incorrect. [21] This can be corrected using the CASTEP OTF potentials, as here, which have a lower cutoff energy of only 340 eV. We also note that a correct phase order is critical for MD simulations; otherwise, HfO 2 may relax to the wrong phases for ferroelectric interfaces. Table 2 shows the relative total energies for the various phases of HfO 2 with various U combinations. Our OTF potentials give the correct phase-energy order, with the c-, t-, o-, and finally the m-HfO 2 phases. This phase-energy order agrees with experiments. [22,23] However, while fitting these phase energies, care must be taken that other factors such as the electron affinity (EA) are not too strongly disturbed. A larger U value only on O p, such as U ¼ (0,9), provides a phase energy closer to the PBE results, as shown in Table 2, but eventually at the expense of a worse EA of 4.7 eV, which affects the interface calculations. The U ¼ (6,6) gives overall better results. A disadvantage of the U(d,p) approach is that a larger Hf U(d) value reduces the energy difference between t-and m-HfO 2 phases.  The ionization potential (IP) of bulk cubic HfO 2 was calculated using a nonpolar (110) HfO 2 slab with a thick vacuum layer. Together with the bandgap of 6.04 eV by optimized GGA þ U(d,p), we obtain the IP of 9.60 eV and EA of 3.56 eV. Note that although the IP/EA is higher than the experimental value, [24] they are much better than traditional PBE calculations, as shown in Table 3. The accuracy of the charge neutrality level (CNL), the branch point of the valence and conduction band states, [24] is also improved. Note that a reasonable CNL value allows correct Schottky barrier heights (SBHs) and defect-free interfaces. It is somewhat anomalous that density functionals are often judged on their bandgaps, with little attention to their electron affinities or IPs. The correct EA is also needed to get the correct energy levels with respect to the vacuum level.

Results
The oxygen vacancy defect calculation uses a supercell of 32 cubic HfO 2 formula units, with the 3 Â 3 Â 3 k-point grid. The optimized U ¼ (6,6) combination is always used. As noted previously, one should do both the geometry relaxation and partial density of states (PDOS) by hybrid functional or by GGA þ U, not a GGA relaxation followed by a PDOS calculation by HSE. [25] Here we show that we can get acceptable defect levels using the optimized GGA þ U(d,p) method. For each charged state, the supercell is fully relaxed at constant volume. The defect formation energy (E f ) is obtained by the well-known formula where E tot ðV O q Þ and E tot ðHfO 2 Þ is the total energy of the charged HfO 2 supercell with charged V o defect and the ideal HfO 2 supercell. The E corr is handled as previously. [25] The use of U(d,p) also provides a reasonable description of the defect levels of the O vacancy V O , the most important of the point defects in HfO 2 . [25][26][27] We see in Figure 3 that GGA þ U(d,p) gives a good defect formation energy of 6.2 eV for the neutral state of V O under O-rich limit. It also gives a reasonable description of the transition states, with (þ2/0) level at 4.06 eV and (0/À2) level at 5.09 eV, as shown in Figure 3, compared with the transition levels calculated in the SX functional (4.0 and 5.5 eV, respectively). [25] Note that the þ1 charge state is not stable in both SX and GGA þ U(d,p) methods, with the negative-U behavior, while is stable by HSE. [27] The most critical aspect of the O vacancy is state localization. Pacchioni and co-workers noted that simpler GGA methods often find defects in insulators such as TiO 2 to have higher symmetry geometry than in experiment. [28] This is a general feature of GGA. Hybrid functional methods correct this error, and Figure 3 shows that this version of GGA þ U also corrects this error for defects in HfO 2 .
The nearest Hf─Hf bond length of cubic HfO 2 in the ideal case is 3.80 Å. Once the oxygen vacancy is induced, its adjacent Hf atoms will relax into the energy-minimum positions. For the neutral vacancy, the four neighboring Hf atoms do not relax  Table 2. Phase-energy comparison with different functionals by firstprinciples calculation. The previous work [23] calculated by GGA functional is taken as an experimental reference. 9) Ref. [  However, a larger relaxation occurs for the charged states with GGA þ U(d,p) calculation. We find that the charged states show the asymmetric relaxation. Taking V o 2þ as an example, the six Hf─Hf bonds are 3.85 Å (two bonds), 3.84 Å, 3.81 Å (two bonds), and 3.79 Å, respectively. This asymmetric relaxation by GGA þ U(d,p) is similar to that seen experimentally, [29] although the degree of relaxation is less obvious. The asymmetric relaxation is significant for the correct electron localization, with the electron mainly localized onto two of the adjacent Hf atoms rather than evenly on four of them. Although it is not ideal compared with the hybrid calculation, that is only localized on one Hf, [25] it is much better than the DFT calculation.
We study the HfO 2 interfaces with metals to test the calculation convergence and found the SBH using the core-level method, as shown in Figure 4. As in previous work, [30,31] we can generate metal-oxide interfaces with both oxygen-terminated polar (100) and nonpolar (110) surfaces, with a range of metals with different work functions [32] and using various surface unit cells to obtain lattice matching, as shown in Figure 4a,b. A pinning factor S ¼ ∂ϕ n / ∂Φ M with S % 0.9 is found for the polar and nonpolar faces, as shown in Figure 4c. This value of S is close to the value found in earlier calculations, [31] but we can now cover a wider range of metal lattice constants. The S is larger than earlier MIGs models. [24,33] This observation derived by GGA þ U(d,p) is also consistent with our previous work which used the GGA functional. [30] Note that the MIGS model is an empirical method and may be less accurate for some cases. The SBH supercell model previously used only small (1 Â 1) supercells, whereas now larger lateral supercells can be used. The point of these calculations is that there is no convergence problem for this GGA þ U(d,p) calculation. Furthermore, we see here there is an offset of 0.4 eV between polar and nonpolar fits to the SBH, which is smaller compared with the study by Tse et al. [30] Also, the SBH value is higher here, caused by the gap opening by GGA þ U(d,p).
This U(d,p) approach also allows us to study the effect of interfaces, grain and domain boundaries, and phase energies on FE and antiferroelectricity (AFE) in HfO 2 and its alloys, [34] without restriction due to the artificial lowering of bandgap in the GGA approach. The kinetics of phase transitions during annealing of HfO 2 are often assumed to follow Ostwald's rule of phases, [35] that the system steps sequentially from the least stable tetragonal phase to the most stable monoclinic phase via the orthorhombic FE phase. However, it was recently noted that the DFT energy for the electrode-dielectric system actually follows total energies once the strain and latticematching conditions are included for the FE oxide and TiN electrode, by the attraction of Ti sites to O sites in HfO 2 . [34] We built the ferroelectric-phase HfO 2 interface with compound metal TiN as the contact metal, as TiN is used in HfO 2 ferroelectric field-effect transistors. [36] To minimize the mismatch, the    [25] The calculated A 1 orbital of the O vacancy in þ2 state by GGA þ U(d,p) is inserted, with its localization mainly on two of the four adjacent Hf atoms, rather than uniformly.
www.advancedsciencenews.com www.pss-rapid.com 450 atoms. A second 30 Å-thick supercell with a vacuum slab on top is also used for comparison. To mimic the experiments, such large ferroelectric interfaces go through the ab initio molecular dynamics (AIMD) simulation for 1 ps at 600 K, followed by later relaxation and PDOS analysis. Only Γ point is used to reduce the computation cost. The GGA þ U(d,p) method and OTF potentials are then tested in an AIMD calculation of ferroelectric (100) HfO 2 between two  www.advancedsciencenews.com www.pss-rapid.com TiN layers or contact with one (111) TiN base layer. For the latter case, the free surface must be terminated with a nonpolar face or passivation layer. Figure 5 shows the final structures and corresponding PDOS of TiN/HfO 2 interface. The main HfO 2 layers are still the orthorhombic phase with the ferroelectric character after the MD process. The layer near TiN is a little distorted with defects but recovers after the relaxation, as shown in Figure 5a. These produce defect gap states spread down from the CBM but rapidly vanish when the layers are away from the interface, as shown in Figure 5b. Note that either the model with two TiN contacts or one TiN base layer gives the similar n-type SBH of %1.5 eV.
Also note that such large interface models can only be studied using the DFT-level calculation such as the accurate GGA þ U(d,p) method here, while impossible for any hybrid functional calculations. We built various TiN/HfO 2 FE and AFE interfaces, with the detailed atomic structures shown in Figure 6. For the FE interfaces, there is a large polarization-induced built-in electric field across the HfO 2 samples and grain boundaries, leading to the shifted band-edge positions in different HfO 2 layers, as shown in Figure 5b. With the AFE interfaces, the polarization field can be eliminated, thus increasing the stabilities. Further investigation on the ferroelectric interface as well as the grain boundaries or domains will be discussed in further detail elsewhere.
This GGA þ U(d,p) method is not only limited to the high-k oxides discussed here but is also applicable for a wide range of solids, whose gap is severally underestimated by DFT calculation. We stress that bandgap is not the only descriptor to evaluate the electronic properties of solids, and attention should also be paid to other important ones, such as some featuring band energies and the correct phase order. As this is still the DFT-level computation task, the accurate GGA þ U(d,p) lessens the burden of massive high-throughput or database-driven material discovery calculation by avoiding using costly hybrid functionals.

Conclusion
We have investigated the electronic structures, defect levels, and interface contact properties of various HfO 2 phases using the GGA þ U(d,p) scheme. It is found that the GGA þ U(d,p) method is more effective to obtain the reasonable calculated results than applying the U term only on cation-d orbitals. This observation proves that we can achieve fast and accurate calculation of electronic structures of solids with a computation time comparable with the DFT level and merits the possibility for high-throughput calculations with higher accuracy.