Phase Properties of Different HfO 2 Polymorphs: A DFT-Based Study

: Background: Hafnium Dioxide (HfO 2 ) represents a hopeful material for gate dielectric thin ﬁlms in the ﬁeld of semiconductor integrated circuits. For HfO 2, several crystal structures are possible, with different properties which can be difﬁcult to describe in detail from an experimental point of view. In this study, a detailed computational approach has been shown to present a complete analysis of four HfO 2 polymorphs, outlining the intrinsic properties of each phase on the basis of atomistic displacements. Methods: Density functional theory (DFT) based methods have been used to accurately describe the chemical physical properties of the polymorphs. Corrective Hubbard (U) semi-empirical terms have been added to exchange correlation energy in order to better reproduce the excited-state properties of HfO 2 polymorphs. Results: the monoclinic phase resulted in the lowest cohesive energy, while the orthorhombic showed peculiar properties due to its intrinsic ferroelectric behavior. DFT + U methods showed the different responses of the four polymorphs to an applied ﬁeld, and the orthorhombic phase was the least likely to undergo point defects as oxygen vacancies. Conclusions: The obtained results give a deeper insight into the differences in excited states phenomena in relation to each speciﬁc HfO 2 polymorph.


Introduction
Silicon Dioxide (SiO 2 ) is well known as a desirable gate insulator, but it shows important limits for maintaining reliability when the thickness of the gate insulator becomes a few nanometers. Current leakage in the Complementary Metal Oxide-Semiconductor (CMOS) is significantly large for very thin SiO 2 films [1,2]. This leakage represents a serious problem from the viewpoint of the reliability of advanced electronic devices and the loss of electric power of materials. In this optic, the development of next-generation CMOS devices represents a very challenging point in the scaling down of their size. Thanks to its relatively high dielectric constant, large band gap, and good thermal and chemical stabilities, Hafnium Dioxide (HfO 2 ) represents one of the most promising candidates among many high-k materials [3,4].
It is known that HfO 2 can adopt different crystal structures in the solid state, along with the chemical physical properties depending on them. The monoclinic structure (space group P2 1 /c) ( Figure 1A) is known to be the most stable, and its dielectric constant is reported to be about 22 [5]. In the monoclinic phase, oxygen atoms can be either threeor fourfold-coordinated, while the Hf atoms can be seven-or eightfold-coordinated. The orthorhombic phase with the space group of Pca2 1 ( Figure 1B) represents another important metastable HfO 2 polymorph; this phase is widely recognized as the most responsible for ferroelectricity behavior [6,7]. On the other hand, higher dielectric constants, reaching about 30 [8,9] have been attributed to other metastable structures, such as the cubic (Fm3m) cubic (Fm3m) ( Figure 1C) and the tetragonal (P42/nmc) ( Figure 1D) examples. In these two phases, oxygen is fourfold coordinated, and the Hf atom is eightfold coordinated. The correlation between the atomic displacements and the phase properties is of extreme significance to better understand the application of HfO2 in the manufacture of advanced electronic devices. In this work, a massive and detailed investigation of the four different HfO2 polymorphs has been conducted. Different aspects of different crystals, i.e., lattice parameters, cohesive energy, density of states (DOS), optical properties, and the effects of oxygen vacancies have been elucidated for each HfO2 polymorph using different DFT methods based on generalized gradient approximation (GGA) for exchange correlation (xc) energy description [10]. Since the DFT approach is reliable for ground state properties predictions, different studies indicate that the resulting HfO2 band gap obtained by the GGA based DFT approach is about 3.65 eV, smaller than the experimental value of 5.7 eV [11]. This means that correlation effects in HfO2 cannot be ignored; for this reason, a computational strategy able to include the correlation effect with reasonable computational effort is represented by the GGA + U approach, in which the correlation effect is incorporated through on-site Coulomb interaction "U". This method is based on the inclusion of some corrective terms that confer a qualitative improvement compared with the GGA based DFT approach for excited state properties such as DOS and optical properties [12]. The GGA + U approach allows reproduction of experimental data in a more accurate way with respect to those obtained with the hybrid functional approach, requiring much lower computational costs [13,14].
For the mentioned reasons, the GGA + U method has been employed in this study in combination with the DFT GGA method. With the aim to better describe these strong Coulomb interactions, not only the on-site U for Hf 5d electrons, but also that for O 2p electrons have been considered [15,16]. Typically, the U correction is only applied on the d and f orbitals of metal entities. However, in certain cases the band gap is still underestimated if compared with the experimental result even with high U values. To obtain the ideal band structure, in addition to the d orbitals in transition metals (U d ), a few theoretical The correlation between the atomic displacements and the phase properties is of extreme significance to better understand the application of HfO 2 in the manufacture of advanced electronic devices. In this work, a massive and detailed investigation of the four different HfO 2 polymorphs has been conducted. Different aspects of different crystals, i.e., lattice parameters, cohesive energy, density of states (DOS), optical properties, and the effects of oxygen vacancies have been elucidated for each HfO 2 polymorph using different DFT methods based on generalized gradient approximation (GGA) for exchange correlation (xc) energy description [10]. Since the DFT approach is reliable for ground state properties predictions, different studies indicate that the resulting HfO 2 band gap obtained by the GGA based DFT approach is about 3.65 eV, smaller than the experimental value of 5.7 eV [11]. This means that correlation effects in HfO 2 cannot be ignored; for this reason, a computational strategy able to include the correlation effect with reasonable computational effort is represented by the GGA + U approach, in which the correlation effect is incorporated through on-site Coulomb interaction "U". This method is based on the inclusion of some corrective terms that confer a qualitative improvement compared with the GGA based DFT approach for excited state properties such as DOS and optical properties [12]. The GGA + U approach allows reproduction of experimental data in a more accurate way with respect to those obtained with the hybrid functional approach, requiring much lower computational costs [13,14].
For the mentioned reasons, the GGA + U method has been employed in this study in combination with the DFT GGA method. With the aim to better describe these strong Coulomb interactions, not only the on-site U for Hf 5d electrons, but also that for O 2p electrons have been considered [15,16]. Typically, the U correction is only applied on the d and f orbitals of metal entities. However, in certain cases the band gap is still underestimated if compared with the experimental result even with high U values. To obtain the ideal band structure, in addition to the d orbitals in transition metals (U d ), a few theoretical studies of the U effect on p orbitals (U p ) in oxygen have been discussed and has been shown to be crucial in oxides [17,18]. The choice of U values is not unambiguous. Though this could be extracted from first-principles calculations, it is nontrivial to determine its value a priori. Hence, U is often determined by fitting a specific set of experimental data, such as band gaps and structural properties and, in general, more than one U value can be used [19]. It is known that the hafnium d orbital is too high in energy. On the other hand, the bottom of the conduction band is mainly composed of the oxygen 2p orbitals. It is very important to note that, since the GGA + U p + U d method aims to correct a localized electron in the d state in a correlated system, adding the parameter Hubbard U in the p-state is also a crucial step. For this reason, in this work Hf U d and O U p values have been respectively settled using on-site repulsion of 6 eV and 4 eV, respectively, fitting these values to the DFT approach adopted. In this way, the GGA + U formalism leads to an improved description of critical aspects such as electronic properties, as already reported in other studies of metal dioxides [20].
Another important aspect to consider is the generation of point defects and impurities, that strongly affect the physical chemical properties of polymorphs and have a decisive impact on their performance in applications as semiconductor devices. For this reason, defects formation has been considered in this work, and in particular, the DFT + U approach has been used again to identify the effects of oxygen vacancies on the structural and electrical properties of HfO 2 polymorphs. The DFT and DFT + U methods have been used to make an exhaustive comparison between the four polymorphs, underlining all possible differences on the basis of the properties of phases. The results obtained have shown a direct correlation between the chemical-physical properties of HfO 2 systems and the atomistic positions, suggesting a way to study in detail the properties of different polymorphs to be used for engineering HfO 2 based devices. Moreover, a direct comparison between the results obtained with and without the U values has been reported to underline the capabilities and the limits of both methods and the important and powerful combination between them.

Materials and Methods
The modelling of all HfO 2 polymorphs and all DFT and DFT + U calculations has been performed using the QuantumATK atomic scale modelling platform [21]. The electron basis has been expanded in linear combination using the atomic orbital (LCAO) method for both Hf and O entities, resembling the SIESTA formalism [22]. The Perdew-Burke-Ernzerhof (PBE) GGA density functional has been used for all simulations to estimate the electron exchange-correlation (xc) energy [23]. The ionic cores have been represented by normconserving (NC) PseudoDojo pseudopotentials for both Hf and O atoms [24]. For Hf, the 5d 2 and 6s 2 electrons have been explicitly treated as valence electrons, while the 2s 2 and 2p 4 electrons were the obvious valence electrons for O atoms. Pseudo-Dojo pseudopotentials are a modern NC type with multipole projectors for each angular momentum, with the aim to ensure high accuracy. Moreover, using this type of NC pseudopotential, the charge density of the electron core has been reproduced within a short-cut-off radius; this aspect is very important, because the xc energy is not linear, so the charge density in the core region could overlap with valence density when U values are included, leading to an underestimation of the xc energy. Correct on-site repulsion values have been calibrated, and the values of 6 eV and 4 eV to Hf d and O p valence electrons, respectively, have been chosen, allowing us to reproduce the bandgap values of HfO 2 ( Figure 2). All polymorphs have been modelled and periodic boundary condition (PBC) have been used along all axes to avoid problems with boundary effects caused by finite size and to reduce the calculation time, maintaining a high accuracy degree. The energy cut-off has been fixed at 1200 eV and the Brillouin-zone integration has been performed over a 15 × 15 × 15 k-points grid for the m-HfO 2 , o-HfO 2 , c-HfO 2 , and t-HfO 2 polymorphs modelled. This set of parameters assure the total energy convergence of 5.0 × 10 −6 eV/atom, the maximum force of 1 × 10 −2 eV/Å, the maximum stress of 2.0 × 10 −2 GPa, and the maximum displacement of 5.0 × 10 −4 Å. The polarization of o-HfO 2 has been obtained by the modern theory of polarization [25] and the Berry phase operator method [26].
Crystals 2022, 12, x FOR PEER REVIEW 4 of 1 maximum displacement of 5.0 × 10 −4 Å. The polarization of o-HfO2 has been obtained b the modern theory of polarization [25] and the Berry phase operator method [26]. The DFT-GGA method has been used as first approach to obtain the optimized ge ometries and to compute the cohesive energy of HfO2 polymorphs. In a second step, th calibrated U values are not indicated for geometry optimization and for cohesive energ calculations because of the explicit inclusion of terms for d electrons for the xc energy which can underestimate this aspect for ground state phenomena such as those mentione above. On the other hand, the inclusion of U values in calculation is suitable for the repro duction of excited states phenomena. For this reason, in the first step the canonical DF approach has been used to compute ground state system properties such as lattice param eters and cohesive energy, and then U d and U p values have been included [27]. The result ing DOS and optical properties of m, c, o, and t-HfO2 obtained by means of GGA + U d U p have been reported. Oxygen vacancies have been explicitly treated for all polymorph putting two O atoms in the unit cell in a ghost configuration. This approach allows mainte nance of the basic functions associated with the atoms selected when removing them t create the vacancies. The formation energy of a defect X in charge state q has been define as: Etot [X q ] is the total energy derived from a supercell calculation containing the defec X, and Etot [HfO2] is the total energy for the perfect crystal using an equivalent supercel The integer ni indicates the number of atoms of type i (host atoms or impurity atoms) tha have been added to (ni > 0) or removed from (ni < 0) the supercell to form the defect, an the μi are the corresponding chemical potentials of these species. Chemical potentials rep resent the energy of the reservoirs with which atoms are being exchanged. The analog o the chemical potential for "charge" is given by the chemical potential of the electrons, i.e the Fermi energy E F . Finally, Ecorr is a correction term that accounts for finite k-point sam pling in the case of shallow impurities, or for elastic and/or electrostatic interactions be tween supercells. All simulations have been carried out using the QuantumATK atomi scale modelling platform. The DFT-GGA method has been used as first approach to obtain the optimized geometries and to compute the cohesive energy of HfO 2 polymorphs. In a second step, the calibrated U values are not indicated for geometry optimization and for cohesive energy calculations because of the explicit inclusion of terms for d electrons for the xc energy, which can underestimate this aspect for ground state phenomena such as those mentioned above. On the other hand, the inclusion of U values in calculation is suitable for the reproduction of excited states phenomena. For this reason, in the first step the canonical DFT approach has been used to compute ground state system properties such as lattice parameters and cohesive energy, and then U d and U p values have been included [27]. The resulting DOS and optical properties of m, c, o, and t-HfO 2 obtained by means of GGA + U d + U p have been reported. Oxygen vacancies have been explicitly treated for all polymorphs putting two O atoms in the unit cell in a ghost configuration. This approach allows maintenance of the basic functions associated with the atoms selected when removing them to create the vacancies. The formation energy of a defect X in charge state q has been defined as: is the total energy derived from a supercell calculation containing the defect X, and E tot [HfO 2 ] is the total energy for the perfect crystal using an equivalent supercell. The integer n i indicates the number of atoms of type i (host atoms or impurity atoms) that have been added to (n i > 0) or removed from (n i < 0) the supercell to form the defect, and the µ i are the corresponding chemical potentials of these species. Chemical potentials represent the energy of the reservoirs with which atoms are being exchanged. The analog of the chemical potential for "charge" is given by the chemical potential of the electrons, i.e., the Fermi energy E F . Finally, E corr is a correction term that accounts for finite k-point sampling in the case of shallow impurities, or for elastic and/or electrostatic interactions between supercells. All simulations have been carried out using the QuantumATK atomic scale modelling platform.

Results
As previously reported, a DFT-GGA approach has been used to calculate the optimal lattice parameters after geometry optimizations and to compute the ground state cohesive energy. The DOS and the electrical properties of the four optimized systems have been investigated using the GGA + U d + U p approach. The same approach has been used again to model oxygen vacancies, to calculate the defect formation energy, and to compute bandgap structure variations for each HfO 2 polymorph.

Structure Analysis
The space group of m-HfO 2 is P2 1 /c, and the optimized values obtained after DFT optimization were 5.068 Å, 5.135 Å, and 5.292 Å for, respectively, a, b, and c vectors of the unit cell ( Table 1). The detected distances between Hf and O were 2.27 Å, 2.15 Å, and 2.04 Å. The optimized values o-HfO 2 (space group Pca2 1 ) of a, b, and c vectors were 5.231 Å, 5.008 Å, and 5.052 Å, respectively, while the detected distances between Hf and O were 2.14 Å and 2.05 Å. The detected value for cubic phase (space group FM3m) was 5.062 Å, and the obtained distance between Hf and O atoms in the cubic phase was 2.2 Å. Finally, for the t-polymorph, (space group P4 2 /nmc) with unit cell consisting of a regular prism with squared basis, the optimized values were in line with the previously described cubic phase, because a = b = 5.062 Å, while c = 5.22 Å (Table 1), and the obtained distances between Hf and O ions were 2.2 Å, as well as in the c-HfO 2 phase.  To be sure about the reliability of the DFT-GGA computational method adopted, a direct comparison with experimental data has been performed, and the effects of U values on the optimization of different HfO 2 polymorphs have been reported. In all cases, an increase of lattice size of 0.2 Å has been observed, leading to results that deviate from the experimental trend.
The PBE-GGA approach has been used again to compute the cohesive energy of all the HfO 2 polymorphs investigated. From the computational evidence the m-HfO 2 was the most stable polymorph (−46.98 eV/unit cell) (Figure 2), and o-, and t-HfO 2 phases have been found to be low energy systems, with cohesive energy values of −46.64 and −45.76 eV/unit cell, respectively. The c-HfO 2 polymorph showed low cohesive energy, with values of −45.51 eV/unit cell.

Electrical Analysis
The optimized lattice parameters obtained from the DFT-GGA calculation have been used as a starting point to detect the effects of atomistic displacements on the electrical properties of HfO 2 polymorphs and then to investigate the effects of oxygen vacancies.
For this purpose, U corrective terms have been included, and a preliminary calibration of the optimal values has been performed. As previously mentioned, the chosen values of 6 eV and 4 eV for U d and U p terms respectively have been applied, and correct bandgap values [32] have been obtained for all the four polymorphs ( Figure 3). values of 6 eV and 4 eV for U d and U p terms respectively have been applied, and cor bandgap values [32] have been obtained for all the four polymorphs ( Figure 3). After the calibration, the PBE-GGA + U approach has been used to identify di ences in the bandgap values, that have been obtained from partial and total DOS com ing the behavior of the HfO2 investigated phases (Figure 4). For all polymorphs, two tures in the valence band have been evidenced. In the m-HfO2 the calculated band value was 5.246 eV. In the o-HfO2 a similar trend with respect to the previous phase been observed, also considering the Hf and O contributions. The calculated bandgap 5.627 eV, and a different trend has been detected for the c-HfO2; two parts in the con tion band have been found, the first from 2.783 to 4.189 eV, the second starting from 5 eV. The resulted band gap value was 4.753 eV. A similar trend in Hf and O contribut has been observed for t-HfO2, but very different values in bandgap have been detected. lower part of the conduction band was between 2.827 eV to 4.216 eV, while the upper started from 5.462 eV, with a resulting bandgap of 6.046 eV.  After the calibration, the PBE-GGA + U approach has been used to identify differences in the bandgap values, that have been obtained from partial and total DOS comparing the behavior of the HfO 2 investigated phases ( Figure 4). For all polymorphs, two features in the valence band have been evidenced. In the m-HfO 2 the calculated bandgap value was 5.246 eV. In the o-HfO 2 a similar trend with respect to the previous phase has been observed, also considering the Hf and O contributions. The calculated bandgap was 5.627 eV, and a different trend has been detected for the c-HfO 2 ; two parts in the conduction band have been found, the first from 2.783 to 4.189 eV, the second starting from 5.538 eV. The resulted band gap value was 4.753 eV. A similar trend in Hf and O contributions has been observed for t-HfO 2 , but very different values in bandgap have been detected. The lower part of the conduction band was between 2.827 eV to 4.216 eV, while the upper one started from 5.462 eV, with a resulting bandgap of 6.046 eV. values of 6 eV and 4 eV for U d and U p terms respectively have been applied, and correct bandgap values [32] have been obtained for all the four polymorphs ( Figure 3). After the calibration, the PBE-GGA + U approach has been used to identify differences in the bandgap values, that have been obtained from partial and total DOS comparing the behavior of the HfO2 investigated phases (Figure 4). For all polymorphs, two features in the valence band have been evidenced. In the m-HfO2 the calculated bandgap value was 5.246 eV. In the o-HfO2 a similar trend with respect to the previous phase has been observed, also considering the Hf and O contributions. The calculated bandgap was 5.627 eV, and a different trend has been detected for the c-HfO2; two parts in the conduction band have been found, the first from 2.783 to 4.189 eV, the second starting from 5.538 eV. The resulted band gap value was 4.753 eV. A similar trend in Hf and O contributions has been observed for t-HfO2, but very different values in bandgap have been detected. The lower part of the conduction band was between 2.827 eV to 4.216 eV, while the upper one started from 5.462 eV, with a resulting bandgap of 6.046 eV. To outline the differences in the phase behavior in response to an applied field, dielectric properties as a function of frequencies have been calculated including inter-band To outline the differences in the phase behavior in response to an applied field, dielectric properties as a function of frequencies have been calculated including inter-band and intra-band transitions of band dispersions. In this case also, both the GGA-PBE and GGA-PBE + U approaches have been used to highlight the importance of the inclusion of U terms to simulate exited state properties of systems. In all cases the electronic static dielectric constant ε 0 underwent variations depending on the HfO 2 polymorph, and a typical behavior of the real part ε 1 (Re (ε)) and the imaginary part ε 2 (Imm (ε)) in function of frequency has been observed.
Based on the complex dielectric functions obtained with GGA-PBE values ( Figure 5), important discrepancies from experimental evidence have been found. In detail, the magnitudes of maximum values of the ε 2 were always higher (size 8, 9, and 12 for t-, o-, and c-HfO 2 polymorphs, respectively) than the experimental values of HfO 2 . Regarding the ε 1 behavior, in this case also the magnitudes obtained along the applied field were always higher than the experimental evidence. and intra-band transitions of band dispersions. In this case also, both the GGA-PBE and GGA-PBE + U approaches have been used to highlight the importance of the inclusion o U terms to simulate exited state properties of systems. In all cases the electronic stati dielectric constant ε0 underwent variations depending on the HfO2 polymorph, and a typ ical behavior of the real part ε1 (Re (ε)) and the imaginary part ε2 (Imm (ε)) in function o frequency has been observed. Based on the complex dielectric functions obtained with GGA-PBE values ( Figure 5) important discrepancies from experimental evidence have been found. In detail, the mag nitudes of maximum values of the ε2 were always higher (size 8, 9, and 12 for t-, o-, and c HfO2 polymorphs, respectively) than the experimental values of HfO2. Regarding the ε behavior, in this case also the magnitudes obtained along the applied field were alway higher than the experimental evidence. Using the GGA-PBE + U approach (Figure 6), the ε0 measured for m-HfO2 was 5.73 and the ε1 increased with frequency obtaining a maximum when the photon energy wa 4.96 eV, corresponding to 1200 THz. The ε2 showed vary low values along the whole rang reported, showing a sensitive increase after 1100 THz. The ε0 for o-HfO2 has been found to be 5.67, while the ε1 exhibited a maximum at 920 THz (3.8 eV), then dropped sharply with the rise of photon energy, showing another peak at 1152 THz (4.76 eV). The ε showed two peaks at 963 and 1164 THz (corresponding to 3.98 and 4.81 eV, respectively with increasing values along with increasing frequencies. Interestingly, the c-HfO showed an ε1 trend that seemed to be similar to that observed for o-HfO2 ( Figure 6A). In particular, the ε0 value (5.89) was in line with the those previously measured, while a peak at 934 THz (3.86 eV) and a maximum value of ε1 in correspondence of 1045 THz (4.32 eV have been observed; then this decreased quickly along with rise in frequency showing another peak at 1127 THz (4.66 eV). The similar behavior between o-and c-HfO2 phase is evident also when considering ε2 trend, which showed three peaks at 952, 1083, and 1149 THz, corresponding to 3.94, 4.48, and 4.75 eV, respectively ( Figure 6B). The tetragona phase showed further different behavior with respect to the previous three HfO2 poly morphs. These discrepancies have to be interpreted in terms of presence of peaks and magnitudes. The measured ε0 was 3.66 and ε1 showed lower values with respect to thos observed for the previous systems, and only one peak at 1103 THz (4.56 eV) has been Using the GGA-PBE + U approach (Figure 6), the ε 0 measured for m-HfO 2 was 5.73, and the ε 1 increased with frequency obtaining a maximum when the photon energy was 4.96 eV, corresponding to 1200 THz. The ε 2 showed vary low values along the whole range reported, showing a sensitive increase after 1100 THz. The ε 0 for o-HfO 2 has been found to be 5.67, while the ε 1 exhibited a maximum at 920 THz (3.8 eV), then dropped sharply with the rise of photon energy, showing another peak at 1152 THz (4.76 eV). The ε 2 showed two peaks at 963 and 1164 THz (corresponding to 3.98 and 4.81 eV, respectively) with increasing values along with increasing frequencies. Interestingly, the c-HfO 2 showed an ε 1 trend that seemed to be similar to that observed for o-HfO 2 ( Figure 6A). In particular, the ε 0 value (5.89) was in line with the those previously measured, while a peak at 934 THz (3.86 eV) and a maximum value of ε 1 in correspondence of 1045 THz (4.32 eV) have been observed; then this decreased quickly along with rise in frequency showing another peak at 1127 THz (4.66 eV). The similar behavior between o-and c-HfO 2 phases is evident also when considering ε 2 trend, which showed three peaks at 952, 1083, and 1149 THz, corresponding to 3.94, 4.48, and 4.75 eV, respectively ( Figure 6B). The tetragonal phase showed further different behavior with respect to the previous three HfO 2 polymorphs. These discrepancies have to be interpreted in terms of presence of peaks and magnitudes.
The measured ε 0 was 3.66 and ε 1 showed lower values with respect to those observed for the previous systems, and only one peak at 1103 THz (4.56 eV) has been observed. The ε 2 denoted an intermediate behavior between those detected for m-and o-HfO 2 , showing a gradual increase from 850 THz (3.51 eV) upwards. Since the GGA + U approach could fail in the description of dielectric functions, to sure about the values obtained including U values, the refractive index of monoclinic HfO has been calculated ( Figure 7). This property has already been simulated [33] and is a ready known from an experimental point of view [34]. This allowed us to better note t reliability of the method used and the consistency of the fitted U values, since the resu obtained from our simulation were in perfect agreement with other calculated and expe imental data, confirming the accuracy of this method.  Since the GGA + U approach could fail in the description of dielectric functions, to be sure about the values obtained including U values, the refractive index of monoclinic HfO 2 has been calculated ( Figure 7). This property has already been simulated [33] and is already known from an experimental point of view [34]. This allowed us to better note the reliability of the method used and the consistency of the fitted U values, since the results obtained from our simulation were in perfect agreement with other calculated and experimental data, confirming the accuracy of this method. Since the GGA + U approach could fail in the description of dielectric functions, to be sure about the values obtained including U values, the refractive index of monoclinic HfO2 has been calculated ( Figure 7). This property has already been simulated [33] and is already known from an experimental point of view [34]. This allowed us to better note the reliability of the method used and the consistency of the fitted U values, since the results obtained from our simulation were in perfect agreement with other calculated and experimental data, confirming the accuracy of this method.

Oxygen Vacancies Effects
With the aim to better describe the properties of different HfO2 phases in the presence of point defects, oxygen vacancies have been modelled and the defect formation energy has been computed for each polymorph. Defect formation energies have been identified (Table 2) computing the total energy formation of a system with and without defects fol-

Oxygen Vacancies Effects
With the aim to better describe the properties of different HfO 2 phases in the presence of point defects, oxygen vacancies have been modelled and the defect formation energy has been computed for each polymorph. Defect formation energies have been identified (Table 2) computing the total energy formation of a system with and without defects following the previously described approach used to compute the ground state cohesive energy of HfO 2 polymorphs. As expected, the values obtained were always positiveotherwise the host crystals would be unstable. From the computational results, the phase that seemed to be most prone to vacant oxygens was the cubic phase, with 0.149 eV, while the phase less likely to suffer this type of defect was the orthorhombic phase, with 0.315 eV.
Finally, the band structures with and without oxygen vacancies have been computed for all HfO 2 polymorphs (Figure 8). DFT + U results showed that the numbers of valence bands and conduction bands significantly increased in all defected systems, and the band gap was lower than that of the pure HfO 2 systems.  From the computational results, the phase that seemed to be most prone to vacant oxygens was the cubic phase, with 0.149 eV, while the phase less likely to suffer this type of defect was the orthorhombic phase, with 0.315 eV.
Finally, the band structures with and without oxygen vacancies have been computed for all HfO2 polymorphs (Figure 8). DFT + U results showed that the numbers of valence bands and conduction bands significantly increased in all defected systems, and the band gap was lower than that of the pure HfO2 systems.

Discussion
The HfO2 can exhibit different polymorphs, and the chemical physical properties of material are directly linked to the atomistic disposition. The polymorphs are defined by crystallographic space group symmetries, leading to the existence of many crystalline forms. Each of these could present peculiar properties that can influence the use of this material in small device fabrication. For this reason, the DFT-GGA approach has been used first for all HfO2 phases to identify the geometry optimization and to compute the ground state cohesive energy. The first polymorph studied was m-P21/c phase, and the obtained lattice parameters outlined an optimal accommodation of the atoms in the unit cell [35,36]. The second model investigated was the o-Pca21 phase. Since more than one orthorhombic phase is known for HfO2, the choice of the specific o-phase is due to the

Discussion
The HfO 2 can exhibit different polymorphs, and the chemical physical properties of material are directly linked to the atomistic disposition. The polymorphs are defined by crystallographic space group symmetries, leading to the existence of many crystalline forms. Each of these could present peculiar properties that can influence the use of this material in small device fabrication. For this reason, the DFT-GGA approach has been used first for all HfO 2 phases to identify the geometry optimization and to compute the ground state cohesive energy. The first polymorph studied was m-P2 1 /c phase, and the obtained lattice parameters outlined an optimal accommodation of the atoms in the unit cell [35,36]. The second model investigated was the o-Pca2 1 phase. Since more than one orthorhombic phase is known for HfO 2 , the choice of the specific o-phase is due to the controversial origin of ferroelectricity for the HfO 2 -based film stems, which is known to be related to the competition between the possible ferroelectric and non-ferroelectric HfO 2 phases [37]. However, though many polar phases, such as the orthorhombic with space groups of Pca2 1 and Pmn21 [38] and the rhombohedral phase with the R3m space group, [39] have been supposed as the possible origins of the ferroelectricity of HfO 2 -based thin films, the orthorhombic phase with the space group Pca2 1 is most widely recognized to be responsible for its ferroelectricity [40,41]. The origin of the ferroelectricity in HfO 2 is related to the formation of the non-centrosymmetric polar orthorhombic phase, as well as the Pca2 1 o-polymorph. The displacement of four polarization-related oxygen atoms per unit cell (Figure 9) leads to the polarization of o-HfO 2 . DFT results showed optimal lattice parameters, detecting the asymmetry and the spontaneous polarization of the optimized model.
Crystals 2022, 12, x FOR PEER REVIEW [39] have been supposed as the possible origins of the ferroelectricity of HfO2-ba films, the orthorhombic phase with the space group Pca21 is most widely recogniz responsible for its ferroelectricity [40,41]. The origin of the ferroelectricity in Hf lated to the formation of the non-centrosymmetric polar orthorhombic phase, as the Pca21 o-polymorph. The displacement of four polarization-related oxygen at unit cell ( Figure 9) leads to the polarization of o-HfO2. DFT results showed optim parameters, detecting the asymmetry and the spontaneous polarization of the op model. The cubic HfO2 phase (space group Fm3m) is fully characterized by a singl constant a. The a value for the perfect bulk c-HfO2 has been obtained after geome mization, and the lattice parameters are in good agreement with other theoretica [42][43][44]. In particular, the DFT approach showed that only the Hf-O distance c along with the structure being preserved. These Hf-O bond lengths are consist those obtained by means of experimental and computational works [45]. The tetra HfO2 has been also considered, since a series of phase changes from the t-phas phase to the o-phase are possible [46][47][48]. The DFT result suggested an optimizatio bond angles leading to a more suitable system configuration within the availab [49]. With the aim to perform a rigorous comparison between DFT and DFT+U res lattice constants of systems have been obtained using both approaches (Table 1). T parison is very important to better describe the effects of the two treatments on simulated data The obtained results suggested that the addition of U terms to th ergy led to an overestimation of the nanoparticle lengths of about 0.2 Å for all poly with an important discrepancy to the experimental data. This is due to the ov wavefunction of the valence electron described with xc energy, including U term eV and U p = 4 eV), with the all-core electron density described with NC PseudoDo dopotentials. This overlap causes an underestimation of xc energy, leading to wr tice parameters. Even if the use of lower U values could give more reliable resu The cubic HfO 2 phase (space group Fm3m) is fully characterized by a single lattice constant a. The a value for the perfect bulk c-HfO 2 has been obtained after geometry optimization, and the lattice parameters are in good agreement with other theoretical values [42][43][44]. In particular, the DFT approach showed that only the Hf-O distance changed along with the structure being preserved. These Hf-O bond lengths are consistent with those obtained by means of experimental and computational works [45]. The tetragonal t-HfO 2 has been also considered, since a series of phase changes from the t-phase or m-phase to the o-phase are possible [46][47][48]. The DFT result suggested an optimization in the bond angles leading to a more suitable system configuration within the available space [49]. With the aim to perform a rigorous comparison between DFT and DFT + U results, the lattice constants of systems have been obtained using both approaches (Table 1). This comparison is very important to better describe the effects of the two treatments on the final simulated data The obtained results suggested that the addition of U terms to the xc energy led to an overestimation of the nanoparticle lengths of about 0.2 Å for all polymorphs, with an important discrepancy to the experimental data. This is due to the overlap of wavefunction of the valence electron described with xc energy, including U terms (U d = 6 eV and U p = 4 eV), with the all-core electron density described with NC PseudoDojo pseudopotentials. This overlap causes an underestimation of xc energy, leading to wrong lattice parameters. Even if the use of lower U values could give more reliable results, this comparison confirmed the reliability of the DFT-GGA approach for the investigation of ground state phenomena as lattice constant values, without the necessity to include corrective terms.
The DFT-GGA approach has been used to compute the cohesive energy of polymorphs. Each HfO 2 phase has a thermodynamic "internal energy", which directly influences the behavior of different HfO 2 systems. The most stable phase detected was the m-polymorph, and this is directly linked to the angles and the bond distances optimized in the previous geometry computation. The ferroelectric o-phase resulted in a metastable phase in bulk crystals since the cohesive energy was close to the monoclinic P2 1 /c polymorph. The revealed distances between atoms and the asymmetric unit cell seems to stabilize the orthorhombic structure. For these reasons, the m-and the o-phases were shown to be the most stable polymorphs [50]. Following a gradual increase in the cohesive energy, m-phase and o-phase would transform to the t-phase and then to the c-phase, which showed the highest cohesive energy, probably because the atoms are placed in a more constrained way due to non-optimized angle bonds. To conclude, from GGA calculation, HfO 2 shows different phase stability on the basis of the atom positions in the space, and this aspect is really important to predict the explicit properties of each polymorph. Even if some energy information is already known, this kind of investigation allowed us to have a starting point to investigate the energy effects of oxygen vacancies on polymorphs, as discussed below.
To better describe the excited state properties of different HfO 2 phases modelled, total and partial DOS have been computed including the U values in the DFT approach. Since the Hubbard U values computed under one set of technical assumptions cannot be used for different DFT calculations, the finding of the optimal values to use with the PseudoDojo library is of useful importance, since no data are present in the literature. In all four HfO 2 polymorphs studied, the Hf 4+ and O 2− entities do not change, meaning that, since U depends on the projectors used as well as other implementations, the same U values can be useful for different polymorphs [50]. To be sure about this, bandgap values with and without U values have been reported, and optimal values have been obtained. Following the calibration of the U values to use, the effective Hubbard U in GGA + U calculation tended to upper shift the minimum conduction band to provide a larger gap between the lowest conduction band and the highest valence band, and this behavior has been observed for all four polymorphs. The chosen values of 6 eV and 4 eV for U d and U p , respectively, allowed us to obtain the correct shift for each polymorph, obtaining the optimal bandgap values found in the literature.
Focusing on the DOS simulated, the band structure tended to have a different appearance for all four systems. This suggests that the properties of HfO 2 are defined by the short range rather than the long-range atomic ordering, and the absolute bandgap values calculated were in very good agreement with experimental evidence. The m-and o-polymorphs showed similar DOS profiles, which means that these first two phases showed some similarities, in line with an intrinsic stability detected in the cohesive energy computation. On the other hand, a lower bandgap value has been found for the c-phase, and this result could be in line with the lower stability of the phase, denoted with a higher cohesive energy then that observed for the previous systems.
Optical properties have been investigated comparing again the DFT and DFT + U approaches. Results showed that the peaks detected along the frequency range were the same using both methods, with an important difference in magnitude, because the optical properties detected with the DFT approach were always higher than the literature data, while the inclusion of U terms allowed reliable values. This result outlines the importance of including a U corrective term at this point of the investigation, confirming that canonical DFT calculations fail in accurate description of excited states phenomena such as electrical properties. Analyzing the differences between the polymorphs via the dielectric function investigation, an interesting similarity has been found between the e 1 and e 2 trend of o-and c-HfO 2 . In both cases, ε 2 had more pronounced peaks around 90, 1050, and 1150 THz. This could be ascribed to the ferroelectric behavior of the o-polymorph, which should influence the response of the phase to the applied field. In fact, the spontaneous polarization of this polymorph is directly reflected to the behavior of the material in response to the simulated applied field. The t-HfO 2 showed lower values of ε 1 , in line with a higher bandgap value with respect to the those observed for the other phases. Results suggested that the ε 1 and ε 2 values depend directly on the polymorph, where ε 1 always assumed more intense values in o-and c-HfO 2 phases, representing the greater polarization capacity of the material in these two phases. To verify the consistency of the dielectric functions simulated, the refractive index of m-HfO 2 has been simulated and compared with other calculated data based on the use of the GGA + U method, and then showing a perfect agreement along the energy range of the applied filed. The magnitudes detected have been compared with experimental data, showing the reliability of the method and its capability to detect different responses of different polymorphs to an applied field.
Finally, the effects of the oxygen vacancies have been investigated following the DFT + U approach already used. The reason for this kind of approach is to identify if different polymorphs may have different aptitudes to suffer defects, because this directly influences the properties of polymorphs, leading to a loss of performance of semiconductor devices. The calculation of defect formation energy includes (i) vibrational (phonon) contributions, because the creation of a defect modifies the chemical bonds and thus the bond strength in its vicinity, (ii) electronic contributions, which are commonly small for semiconductors but can be sizable for metals, and (iii) magnetic excitations. As expected, different responses have been detected on the basis of different atoms' displacement, and in particular the c-phase was the most prone to undergo vacant oxygens. The lower stability of the polymorph with respect to the other models could lead an intrinsic instability that is directly linked to a major propensity of the oxygen vacancies. On the other hand, the o-phase was the least subject to oxygen vacancies, possibly ascribed to the ferroelectric behavior of the o-polymorph, since the oxygen entities are directly involved in this behavior and the spontaneous polarization is directly linked to a strong overlap between the wavefunctions. The bandgaps structures showed important changes when the oxygen vacancies occurred, and the band gap value is narrowed along with increase in oxygen vacancy concentrations due to the delocalized states in the valence band. In more detail, the overlap of the nonlocalized oxygen vacancy states with valence band raises the position of the valence band, resulting in a band gap narrowing. These results clearly show that the band gap of all HfO 2 phases decreased due to oxygen vacancies, and that all HfO 2 polymorphs investigated showed a similar response to the oxygen vacancy onset.

Conclusions
HfO 2 shows many polymorphs with different material properties. This large variation promises many prominent applications of HfO 2 based devices. This work shows the potentialities of the atomistic simulation method exposed and the possibility of obtaining more and different information. In particular, a complete comparison of four different polymorphs has been reported using GGA and GGA + U methods in combination. Since some information is already known, a complete comparison of all the properties investigated between the four most important HfO 2 most polymorphs is new. Moreover, the capabilities and the limits of the GGA and GGA + U methods used have been underlined, laying the foundations for their correct use and providing the elements to understand when to use different approaches, showing the correct way to combine them and their impact on the reliability of the final results. The optimal balance between the accuracy degree and reasonable simulation time required makes this combined method a valid alternative to more expensive simulation methods such as those based on semi-hybrid functionals, or purely ab-initio. The differences between monoclinic (P2 1 /c), orthorhombic (PCa2 1 ), cubic (Fm3m), and tetragonal (P4 2 /nmc) HfO 2 polymorphs have been exhaustively discussed on the basis of transition phenomena between them. The DFT-GGA approach has been used to obtain lattice parameters of geometry optimized structures, and the results coincided well with the experimental value and other theoretical results. In addition, the computation of the cohesive energy, which depends on the specific polymorph, revealed that m-and o-HfO 2 phases are observed as the most stable.
To better describe electrical properties that are correlated to the excited states of different polymorphs, the on-site Coulomb interactions have been included for 5d orbital of Hf atom and 2p orbital of O atom with the aim to solve the discrepancies between experimental and predicted excited states properties, improving the description of the electronic properties. Optimal bandgap values for all the four polymorphs have been obtained with U d = 6 eV and U p = 4 eV, suggesting that these values are reliable to use in combination with the settled approach. The electron structure and optical properties of HfO 2 phases have been calculated by the calibrated GGA + U method. Since the atomistic composition is always the same, the results obtained suggested different densities of states and bandgap values in relation to the polymorph, confirming the importance of investigating the atomistic effects on the materials properties. Differences between polymorphs have been also observed in the measuring of static dielectric constant value (ε 0 ) and in the function of an applied field, since ε 1 and ε 2 had different trends in relation to the phase. Different ε 0 values have been detected and features in ε 1 and ε 2 values have been identified in the high frequencies range for o-and c-HfO 2 polymorphs. To better identify the role of U terms, lattice parameters and optical properties have been calculated using GGA and GGA + U, and the results have been compared and discussed.
Finally, the role of oxygen vacancies as electron acceptors has been investigated. The results suggested that oxygen vacancies play important roles in modulating the electrical activity of the HfO 2 nanostructures, prompting a decrease in bandgap of different polymorphs. The defect formation calculation showed a peculiar propensity of the c-HfO 2 phase for oxygen vacancies, while the o-HfO 2 resulted was the least prone to these kinds of defect. The changing in the bandgap structures did not outline different responses of different phases, showing an increase in the number of the valence and conduction bands and a narrowing in the bandgap values for all systems. All these findings provide further insight into the influences of oxygen defects observed in HfO 2 structures outlining the importance of the first principal calculations; in particular, the combined approach of GGA and GGA + U methods can be used as a guideline for engineering HfO 2 based devices, accurately predicting their properties before production.