The Role of Zr on Monoclinic and Orthorhombic HfxZryO2 Systems: A First-Principles Study

HfO2 shows different polymorphs, including monoclinic and orthorhombic ones, that exhibit singular properties. Moreover, the character of HfO2 is also influenced by the Zr atoms as a doping agent. Here, an extensive study of the monoclinic P21/c and the orthorhombic Pca21 polymorphs of HfO2, Hf0.75Zr0.25O2, and Hf0.5Zr0.5O2 is reported. For all six systems, density functional theory (DFT) methods based on generalized gradient approximations (GGAs) were first used; then the GGA + U method was settled and calibrated to describe the electrical and optical properties of polymorphs and the responses to the oxygen vacancies. Zr had different effects in relation to the polymorph; moreover, the amount of Zr led to important differences in the optical properties of the Pca21 polymorph. Finally, oxygen vacancies were investigated, showing an important modulation of the properties of HfxZryO2 nanostructures. The combined GGA and GGA + U methods adopted in this work generate a reasonable prediction of the physicochemical properties of o- and m-HfxZryO2, identifying the effects of doping phenomena.


Introduction
The discovery of ferroelectricity in hafnium oxides (HfO 2 ) thin films [1] has been of crucial importance, thanks to the possibility of overcoming many of the problems related to perovskites-based field-effects transistor (FET) technology. Since the introduction in a manufacturing process by Intel in 2007, HfO 2 has been vastly used as a high-k material, displaying a full CMOS (complementary metal-oxide-semiconductor) compatibility. The ferroelectric behavior originates from the development of nanoscale thin films [2]. The nature of the ferroelectricity is most likely related to the formation of a non-centrosymmetric polar crystal phase that is stable under specific conditions. As already mentioned, the ferroelectricity behavior of HfO 2 is similar to that of perovskites materials but with a few peculiarities: HfO 2 has a relatively wide bandgap (5.3 eV), large band offset with Si (less parasitic leakage), large coercive field (1 MVcm −1 ), and low permittivity. Importantly, the absence of an interfacial dead layer in HfO 2 makes this material a promising candidate in thin-film technology, in contrast to perovskites-based material. Moreover, the mid-range dielectric constant allows the switching at reasonable voltage, even though the field strength required to obtain the reversal polarization is higher if compared to perovskites [3][4][5].
When an external electric field is applied to a ferroelectric material, a permanent polarization and hysteresis are induced; such an effect was found for Si-doped [1], Aldoped [6], and Y-doped HfO 2 films [7]. However, HfO 2 has been continuously studied to improve ferroelectricity and to extend its applicability in the semiconductor industry. In this context, the doping of HfO 2 with Zr atoms was effective in obtaining ferroelectricity, although neither the pure HfO 2 nor the pure ZrO 2 films exhibit the same singular properties. It is well-known that the monoclinic structure (space group P2 1 /c) of HfO 2 is the most stable; however, the presence of Zr determines the formation of an orthorhombic phase (Pca2 1 non-centrosymmetric space groups) which is closely related to ferroelectricity [8].
In the present work, we report an extensive theoretical study regarding the two mentioned polymorphs by doping different doses of Zr. Thus, six different systems were studied, the monoclinic and orthorhombic polymorphs of HfO 2 , Hf 0.75 Zr 0.25 O 2 , and Hf 0.5 Zr 0.5 O 2 ( Figure 1).  [6], and Y-doped HfO2 films [7]. However, HfO2 has been continuously studied to improve ferroelectricity and to extend its applicability in the semiconductor industry. In this context, the doping of HfO2 with Zr atoms was effective in obtaining ferroelectricity, although neither the pure HfO2 nor the pure ZrO2 films exhibit the same singular properties. It is well-known that the monoclinic structure (space group P21/c) of HfO2 is the most stable; however, the presence of Zr determines the formation of an orthorhombic phase (Pca21 non-centrosymmetric space groups) which is closely related to ferroelectricity [8].
In the present work, we report an extensive theoretical study regarding the two mentioned polymorphs by doping different doses of Zr. Thus, six different systems were studied, the monoclinic and orthorhombic polymorphs of HfO2, Hf0.75Zr0.25O2, and Hf0.5Zr0.5O2 (Figure 1). A wide comparison of different properties was performed by using the density functional theory (DFT) approach based on generalized gradient approximation (GGA) exchange-correlation functional [9]. Even if the DFT approach is known to be reliable for ground-state properties' predictions, different studies indicate that conventional exchange-correlation (xc) functionals often underestimate the bandgap of semiconductors concerning the experimental studies [10]. To overcome this problem, a computational strategy that is able to include the correlation effect with reasonable computational effort is the GGA + U approach. Using this method, the correlation effect is incorporated through the on-site Coulomb interaction Hubbard (U) value that needs to be calibrated. The inclusion of corrective terms confers a qualitative improvement in comparison to the GGA based DFT approach for excited-state properties (such as DOS and optical properties [11]).
The on-site U for Hf 5 d , Zr 4 d , and O 2 p electrons was considered even if the U correction is typically only applied on the d and f orbitals of metal entities. The inclusion of 2p electrons of O entities in the U treatment allows us to better describe these strong Coulomb interactions, making the U p values of O to be crucial in metal oxides [12,13]. Hf U d , Zr U d , and O U p values were settled by using on-site repulsion of 6, 5.8, and 4 eV, respectively, fitting these values to the DFT approach adopted. A wide comparison of different properties was performed by using the density functional theory (DFT) approach based on generalized gradient approximation (GGA) exchange-correlation functional [9]. Even if the DFT approach is known to be reliable for ground-state properties' predictions, different studies indicate that conventional exchangecorrelation (xc) functionals often underestimate the bandgap of semiconductors concerning the experimental studies [10]. To overcome this problem, a computational strategy that is able to include the correlation effect with reasonable computational effort is the GGA + U approach. Using this method, the correlation effect is incorporated through the on-site Coulomb interaction Hubbard (U) value that needs to be calibrated. The inclusion of corrective terms confers a qualitative improvement in comparison to the GGA based DFT approach for excited-state properties (such as DOS and optical properties [11]).
The on-site U for Hf 5 d , Zr 4 d , and O 2 p electrons was considered even if the U correction is typically only applied on the d and f orbitals of metal entities. The inclusion of 2p electrons of O entities in the U treatment allows us to better describe these strong Coulomb interactions, making the U p values of O to be crucial in metal oxides [12,13]. Hf U d , Zr U d , and O U p values were settled by using on-site repulsion of 6, 5.8, and 4 eV, respectively, fitting these values to the DFT approach adopted.
For all of these reasons, the combined employment of DFT and DFT + U functional was adopted to give the most reasonable energy stability and structural parameters of the systems. In detail, the effects of Zr doping on lattice parameters were elucidated by using the DFT-GGA approach, and then the electrical and the optical properties of all six models were calculated by adopting the DFT-GGA + U method.
It is well-established that oxide defects adversely affect the functionality and reliability of a wide range of semiconductor devices strongly affecting the physical-chemical properties of the material and thus decisively impact the performances of the respective devices. With this aim, the defects' formation was also considered with the DFT + U approach to identify the effects of oxygen vacancies on the structural and electrical properties of HfO 2 and Hf x Zr y O 2 polymorphs. Since the PCa2 1 phase and the Zr amount are correlated to ferroelectricity, differences in the defects' propensities are expected between polymorphs.

Materials and Methods
A quantum ATK atomic-scale modeling platform was used to model all polymorphs and to perform all calculations [14]. Six structures, namely three monoclinic (P2 1 /c) and three orthorhombic (Pca2 1 ) polymorphs, were modeled: (i) HfO 2 with 0% of Zr substitution, (ii) Hf 0.75 Zr 0.25 O 2 with 25% of Zr substitution, and (iii) Hf 0.5 Zr 0.5 O 2 with 50% of Zr substitution with respect to the total Hf amount. To discover the effects of different amounts of Zr in each polymorph, the Hf ions were randomly substituted, and all six models were simulated. The electron basis was expanded in linear combination, using the atomic orbital (LCAO) method for Hf, Zr, and O entities, resembling the SIESTA formalism [15]. All simulations were carried out by using the Perdew-Burke-Ernzerhof (PBE) GGA density functional for the electron xc energy [16]. For each atom, the ionic cores are represented by norm-conserving (NC) PseudoDojo (PDj) pseudopotentials [17]. The 5d 2 and 6s 2 electrons of Hf and 4d 2 and 5s 2 electrons for Zr are explicitly treated as valence, while the 2s 2 and 2p 4 electrons are the obvious valence electrons for O atoms. Correct on-site repulsion values were calibrated, and then the Hubbard (U) values of 6, 5.8, and 4 eV for Hf d , Zr d , and O p valence electrons were chosen, respectively, to reproduce the bandgap values of HfZrO 2 . Since the description of charge density in the core region could overlap valence density when U values are included, the choice of the correct NC PDj is greatly important because the charge density of the electron core was reproduced within a short cutoff radius, thus avoiding an underestimation of the xc energy. To model the polymorphs, the periodic boundary conditions (PBCs) were used along all axes; in this way, it is possible to avoid problems with boundary effects caused by the finite size and to reduce the calculation time while maintaining high accuracy. The energy cutoff was fixed at 1200 eV, and the Brillouin-zone integration was performed over a 15 × 15 × 15 k-points grid for the modeled P21/c and Pca21 polymorphs. These parameters assure the total energy convergence of 5.0 × 10 −6 eV/atom, the maximum stress of 2.0 × 10 −2 GPa, and the maximum displacement of 5.0 × 10 −4 Å. The polarization of o-Pca2 1 polymorph was obtained by the modern theory of polarization [18] and the Berry phase operator method. The total polarization is assumed to be the sum of the electronic and ionic contributions. The first one (P i ) is calculated by using a simple classical electrostatic sum of point charges, as shown in Equation (1): where Ω is the unit cell volume, Z v ion is the valence charge, and r ν is the position vector of the ν atom.
The second contribution (P e ) to the polarization is obtained as Equation (2): where the sum runs over occupied bands, and k is parallel to the direction of polarization. The G term is a reciprocal lattice vector in the same direction. The states U k,n > are the cell-periodic parts of the Bloch functions, ψ k,n (r) = u k,n (r) e ikr . The last integral is known as the Berry phase [19]. Oxygen vacancies were explicitly treated for each unit cell by putting two O atoms in a ghost configuration. This approach allows for the maintenance of the basic functions associated with the atoms selected when removing them to create the vacancies. The energy formation of a defect X in charge state q is defined as Equation (3): where E tot [X q ] is the total energy derived from a supercell calculation containing the defect X, and E tot [HfZrO 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 were added to (n i > 0) or removed from (n i < 0) the supercell to form the defect, and the µ i corresponds to the 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 elastic and/or electrostatic interactions between supercells.

Results and Discussions
As previously reported, the GGA/PBE and GGA/PBE + U methods were used in combination to calculate the different properties of six systems. The results are organized as follows: The geometry optimizations, lattice parameters, ground state cohesive energy formation, and tensile stress values were calculated by using the DFT approach and are reported in the first part. Then the second section is dedicated to the fitting of the U values to further use. Next, the electrical and optical properties of the six optimized systems are investigated by using the GGA + U d + U p approach, and the respective results are reported in the third section. In the fourth section, the same approach is used to model oxygen vacancies, calculate the defect formation energy, and compute bandgap structure variations for each system.

Effects of Zr on Lattice Parameters
HfO 2 and Hf x Zr y O 2 are characterized by different polymorphs, the monoclinic (m−) and orthorhombic (o−) ones, with space groups of P2 1 /c and Pca2 1 , respectively ( Figure 1). Each of the two polymorphs was used to evaluate the effect of the presence of Zr. Table 1 reports the values of the calculated lattice parameters for the monoclinic and orthorhombic configurations of (i) HfO 2 with 0% of Zr substitution, (ii) Hf 0.75 Zr 0.25 O 2 , and (iii) Hf 0.5 Zr 0.5 O 2 with 25% and 50% of Zr substitution to the total amount of Hf ions, respectively. Lattice energy minimized for m-HfO 2 was obtained by the optimization of the atomic positions and systematically altering the size and angle of the unit cell. After optimization of the lattices, unit cell dimensions of 5.068, 5.135, and 5.292 Å were found for a, b, and c, respectively. The lattice vectors slightly decreased by imposing Zr in the system.
The Pca2 1 polymorph is directly correlated to the ferroelectricity [22], which is due to the formation of the non-centrosymmetric polar orthorhombic phase. The displacement of four polarization-related oxygen atoms per unit cell leads to the polarization of the phase. The optimized lattices for o-HfO 2 were 5.231, 5.008, and 5.052 Å for a, b, and c vectors, respectively; this is in line with other results [23].
The effect of Zr on the Pca2 1 unit cell is generally more evident than that observed for the P2 1 /c polymorph, and, in particular, a sensitive decrease in the c vector was detected in o-Hf 0.5 Zr 0.5 O 2 . Since the c vector follows the direction of spontaneous polarization of the phase, these data suggest that Zr has a significant role to modulate the ferroelectric property of this polymorph.
In Figure 2, the m-phases of HfO 2 , Hf 0.75 Zr 0.25 O 2 , and Hf 0.5 Zr 0.5 O 2 systems are reported, whose bond angles undergo only moderate changes when Zr replaces Hf atoms. Inside the unit cell, the angle Hf-O-Hf (106.71 • ) slightly increased in Zr-O-Hf (108.20 • and 108.34 • for 25% and 50% of Zr substitution, respectively). The same trend, but less pronounced, is observed for the O-Hf-O angle (73.11 • ), which increased in O-Zr-O (73.43 • and 74.24 • for 25% and 50% of Zr substitution, respectively). This assumption is in line with the negligible differences found for the cell vectors as the amount of Zr increases. More substantial effects can be observed for the o-polymorphs, where the Zr intercalation has a high impact on the bond angles with the asymmetric O entities ( Figure 2). This influence is directly linked to the Zr amount in the Hf 0.5 Zr 0.5 O 2 , and the results show an important rearrangement of atom coordinates along the c vector. Based on this evidence, the replacement of 50% of Hf with Zr atoms has the peculiar role of stabilizing the o-polymorph while preserving the necessary asymmetry for the ferroelectric behavior.
Lattice energy minimized for m-HfO2 was obtained by the optimization of the atomic positions and systematically altering the size and angle of the unit cell. After optimization of the lattices, unit cell dimensions of 5.068, 5.135, and 5.292 Å were found for a, b, and c, respectively. The lattice vectors slightly decreased by imposing Zr in the system.
The Pca21 polymorph is directly correlated to the ferroelectricity [22], which is due to the formation of the non-centrosymmetric polar orthorhombic phase. The displacement of four polarization-related oxygen atoms per unit cell leads to the polarization of the phase. The optimized lattices for o-HfO2 were 5.231, 5.008, and 5.052 Å for a, b, and c vectors, respectively; this is in line with other results [23].
The effect of Zr on the Pca21 unit cell is generally more evident than that observed for the P21/c polymorph, and, in particular, a sensitive decrease in the c vector was detected in o-Hf0.5Zr0.5O2. Since the c vector follows the direction of spontaneous polarization of the phase, these data suggest that Zr has a significant role to modulate the ferroelectric property of this polymorph.
In Figure 2, the m-phases of HfO2, Hf0.75Zr0.25O2, and Hf0.5Zr0.5O2 systems are reported, whose bond angles undergo only moderate changes when Zr replaces Hf atoms. Inside the unit cell, the angle Hf-O-Hf (106.71°) slightly increased in Zr-O-Hf (108.20° and 108.34° for 25% and 50% of Zr substitution, respectively). The same trend, but less pronounced, is observed for the O-Hf-O angle (73.11°), which increased in O-Zr-O (73.43° and 74.24° for 25% and 50% of Zr substitution, respectively). This assumption is in line with the negligible differences found for the cell vectors as the amount of Zr increases. More substantial effects can be observed for the o-polymorphs, where the Zr intercalation has a high impact on the bond angles with the asymmetric O entities ( Figure 2). This influence is directly linked to the Zr amount in the Hf0.5Zr0.5O2, and the results show an important rearrangement of atom coordinates along the c vector. Based on this evidence, the replacement of 50% of Hf with Zr atoms has the peculiar role of stabilizing the o-polymorph while preserving the necessary asymmetry for the ferroelectric behavior.  In order to gain deeper insight into the behavior of the phases and to underline the role of Zr, the PBE-GGA approach was again used to compute the ground-state cohesive energy of the polymorphs ( Figure 3A). The direct comparison between m-and o-HfO 2 phases outlined higher stability for the m-polymorph, as already reported in Reference [24]; negligible energetic variations were found for P2 1 /c phases when Zr replaces with Hf atoms. Therefore, the zirconium ion exhibits an excellent doping agent for the hafnium-based oxides [9]. An opposite trend was observed for the orthorhombic polymorphs, in which the Zr doping drastically decreased the formation energy. This effect appeared to be Zrcontent dependent; that is, the higher the Zr percentage, the lower the formation energy, confirming the high impact of zirconium on the phase properties and the stabilization of the ferroelectric crystal.
In order to gain deeper insight into the behavior of the phases and to underline the role of Zr, the PBE-GGA approach was again used to compute the ground-state cohesive energy of the polymorphs ( Figure 3A). The direct comparison between m-and o-HfO2 phases outlined higher stability for the m-polymorph, as already reported in Reference [24]; negligible energetic variations were found for P21/c phases when Zr replaces with Hf atoms. Therefore, the zirconium ion exhibits an excellent doping agent for the hafniumbased oxides [9]. An opposite trend was observed for the orthorhombic polymorphs, in which the Zr doping drastically decreased the formation energy. This effect appeared to be Zr-content dependent; that is, the higher the Zr percentage, the lower the formation energy, confirming the high impact of zirconium on the phase properties and the stabilization of the ferroelectric crystal. The stabilization effect of Zr atoms to the orthorhombic structure is also confirmed by internal stress values ( Figure 3B), which quantify the internal constriction factors in the unit cell. The m-polymorphs do not show a variation of the internal stress after the Zr inclusion, as is predictable by the minor changes in the lattice parameters when moving from undoped HfO2 to Hf0.75Zr0.25O2 and Hf0.75Zr0.25O2. On the contrary, DFT calculations detected internal stress of unit cells for the o-HfO2 (0% of Zr); the stress value drastically decreases when zirconium is included, reaching 7.6 × 10 −4 eV/Å at 50% of Hf replaced by Zr.
Therefore, the increment on Zr concentration seems to induce a structural transition and an internal stabilization for the orthorhombic polymorph. These results are consistent with experimental data indicating tensile stress on the ferroelectric films [25] and showing that a composition close to 50% Hf and 50% Zr tends to favor ferroelectricity [26].

DFT + U Calibration
As previously described, the DFT approach is used to compute lattice parameters, cohesive energy formation, and internal stress values. Since the GGA-PBE method is not reliable to compute excited state properties as dielectric functions, the inclusion of corrective terms to the canonical DFT approach is necessary to obtain the partial density of states (PDOS) and bandgap variations induced by oxygen vacancies. Through the inclusion of on-site Coulomb interaction Hubbard (U) terms, it is possible to better describe the chemical bonds involving d and p electrons of Hf, Zr, and O atoms and to improve qualitatively the electronic descriptions. Zr U d , Hf U d , and O U p values were settled to 5.8, 6, and 4 eV, respectively. To avoid non-physical states, the U values were calibrated by computing the already known bandgap values of m-and o-HfZrO2 polymorphs. Since the U values fitted Therefore, the increment on Zr concentration seems to induce a structural transition and an internal stabilization for the orthorhombic polymorph. These results are consistent with experimental data indicating tensile stress on the ferroelectric films [25] and showing that a composition close to 50% Hf and 50% Zr tends to favor ferroelectricity [26].

DFT + U Calibration
As previously described, the DFT approach is used to compute lattice parameters, cohesive energy formation, and internal stress values. Since the GGA-PBE method is not reliable to compute excited state properties as dielectric functions, the inclusion of corrective terms to the canonical DFT approach is necessary to obtain the partial density of states (PDOS) and bandgap variations induced by oxygen vacancies. Through the inclusion of on-site Coulomb interaction Hubbard (U) terms, it is possible to better describe the chemical bonds involving d and p electrons of Hf, Zr, and O atoms and to improve qualitatively the electronic descriptions. Zr U d , Hf U d , and O U p values were settled to 5.8, 6, and 4 eV, respectively. To avoid non-physical states, the U values were calibrated by computing the already known bandgap values of m-and o-HfZrO 2 polymorphs. Since the U values fitted under one set of technical assumptions could not be used for different DFT approaches, identifying the optimal values to use together with specific pseudopotentials represents a crucial step. It is useful to emphasize that there are no data reported in the literature about the use of the PseudoDojo library for the proposed HfZrO 2 polymorph. Even if the coordinates and the structural parameter effects are different in the HfZrO 2 polymorphs under study, the Hf 4+ , Zr 4+ , and O 2entities do not change. This means that, since U depends on the projectors used, in our case, the same U values can be applied for more than one polymorph [27]. To be sure about the use of univocal U values for both polymorphs, the bandgap values calculated, including/excluding U values, were reported and compared to one another, and then the optimal values were obtained (Figure 4). Following the calibration of the U values, the effective Hubbard in GGA + U calculation tended to upper shift the minimum conduction band, providing a larger gap between the lowest conduction band and the highest valence band; this behavior was observed for both HfZrO 2 polymorphs. With the previously reported U, the bandgap values obtained were 5.68 and 5.76 eV for the m-and o-phases, respectively; these results are in line with the experimental results of Reference [28]. These calculations confirmed the reliability of the fitted U values to compute the excited states' properties for both of the polymorphs.
identifying the optimal values to use together with specific pseudopotentials represents a crucial step. It is useful to emphasize that there are no data reported in the literature about the use of the PseudoDojo library for the proposed HfZrO2 polymorph. Even if the coordinates and the structural parameter effects are different in the HfZrO2 polymorphs under study, the Hf 4+ , Zr 4+ , and O 2-entities do not change. This means that, since U depends on the projectors used, in our case, the same U values can be applied for more than one polymorph [27]. To be sure about the use of univocal U values for both polymorphs, the bandgap values calculated, including/excluding U values, were reported and compared to one another, and then the optimal values were obtained (Figure 4). Following the calibration of the U values, the effective Hubbard in GGA + U calculation tended to upper shift the minimum conduction band, providing a larger gap between the lowest conduction band and the highest valence band; this behavior was observed for both HfZrO2 polymorphs. With the previously reported U, the bandgap values obtained were 5.68 and 5.76 eV for the m-and o-phases, respectively; these results are in line with the experimental results of Reference [28]. These calculations confirmed the reliability of the fitted U values to compute the excited states' properties for both of the polymorphs.

Electrical and Optical Properties Calculation
In the next step, DFT + U approximation was used to investigate the electrical and optical properties of m-P21/c and o-Pca21 HfxZryO2 polymorphs.
To determine the relative values of charges from wave functions, focusing on a qualitative point of view, the Mulliken population analysis [29] was performed ( Table 2). For the P21/c HfO2, the transferred charge from the Hf to O atoms is 0.840 atomic unit (a.u.). Following the doping, a proportional decrease in charge transferred into O is evident with 0.826 and 0.812 for 25% and 50% doped m-polymorphs, respectively. A similar trend was observed for Pca21 polymorphs when Zr was absent. Instead, when Zr was included in the polymorph, the charge transfer rose remarkably in O atoms, with 0.882 and 0.981 for 25% and 50% doped o-polymorphs, respectively. These results suggest that Zr incorporation into crystals leads to an increase in charge transfer from Hf and Zr into O atoms, and this effect is evident in o-phases, in particular, on Hf0.5Zr0.5O2.

Electrical and Optical Properties Calculation
In the next step, DFT + U approximation was used to investigate the electrical and optical properties of m-P2 1 /c and o-Pca2 1 Hf x Zr y O 2 polymorphs.
To determine the relative values of charges from wave functions, focusing on a qualitative point of view, the Mulliken population analysis [29] was performed ( Table 2). For the P2 1 /c HfO 2 , the transferred charge from the Hf to O atoms is 0.840 atomic unit (a.u.). Following the doping, a proportional decrease in charge transferred into O is evident with 0.826 and 0.812 for 25% and 50% doped m-polymorphs, respectively. A similar trend was observed for Pca2 1 polymorphs when Zr was absent. Instead, when Zr was included in the polymorph, the charge transfer rose remarkably in O atoms, with 0.882 and 0.981 for 25% and 50% doped o-polymorphs, respectively. These results suggest that Zr incorporation into crystals leads to an increase in charge transfer from Hf and Zr into O atoms, and this effect is evident in o-phases, in particular, on Hf 0.5 Zr 0.5 O 2 .
The GGA/PBE + U approach was also used to identify differences in bandgap values by computing the partial density of states (PDOS) and to compare the behavior of the different phases ( Figure 5). For all systems, only one feature in the valence band was evidenced. Between an energy range of −7.5 and −2 eV, the main contribution of valence band maximum (VBM) originated from O 2 p orbitals. On the other hand, the energy ranges from 2.2 eV upward represent the conduction band maximum (CBM), which was mainly composed of Hf 5 d orbitals. Upon doping into the two pure systems by introducing dopant Zr atoms, the VBM is still mainly attributed to O 2 p , while the CBM depends on Hf 5 d and Zr 4 d orbitals. More, with the inclusion of Zr, impure states are formed between the Fermi level and CBM ( Figure 5E,F).  In general, the intensity of the PDOS peaks is sensitive to the geometrical structure around the Hf atoms. Obviously, as the Zr concentration increases, the CMB is divided into two parts from the energy range of around 5.00 eV. Our inspection of Figure 5C,D reveals that Hf 3d and Zr 3d orbitals have more contributions in PDOS plots, while O 2p possesses less content. Furthermore, the Zr defect enhanced the peak intensity for Pca21 more than the P21/c HfO2 polymorph. As Zr doping increased from 25% to 50% ( Figure 5E,F), similar results appeared in the CMB state of both polymorphs by dropping peak intensities; however, the intensity of the Zr 3d peak in P21/c HfO2 became weaker than one calculated for the Pca21 HfO2 structure, meaning that the former compound is more affected by In general, the intensity of the PDOS peaks is sensitive to the geometrical structure around the Hf atoms. Obviously, as the Zr concentration increases, the CMB is divided into two parts from the energy range of around 5.00 eV. Our inspection of Figure 5C,D reveals that Hf 3d and Zr 3d orbitals have more contributions in PDOS plots, while O 2p possesses less content. Furthermore, the Zr defect enhanced the peak intensity for Pca2 1 more than the P2 1 /c HfO 2 polymorph. As Zr doping increased from 25% to 50% ( Figure 5E,F), similar results appeared in the CMB state of both polymorphs by dropping peak intensities; however, the intensity of the Zr 3d peak in P2 1 /c HfO 2 became weaker than one calculated for the Pca2 1 HfO 2 structure, meaning that the former compound is more affected by the presence of Zr, while the latter structure undergoes slight changes.
One of the most interesting behaviors from PDOS is the change in the Zr contribution, moving from Hf 0.75 Zr 0.25 O 2 to Hf 0.5 Zr 0.5 O 2 in o-polymorph. The VBM and CBM energies were also chosen to estimate the bandgap values. The monoclinic structures are characterized by bandgap values of 5.24, 5.62, and 5.68 eV for 0%, 25%, and 50% of Zr doses. In the o-HfO 2 , a similar trend was observed that is relevant to the Hf and O contributions. The extrapolated bandgap is 5.63, 5.81, and 5.76, moving from 0%, 25%, and 50% of Zr contributions, respectively.
The effects of Zr on the real part (ε 1 ) and imaginary part (ε 2 ) of dielectric function for m-and o-HfZrO 2 were calculated on a wide energy range ( Figure 6). As reported, the ε 1 value measured for m-polymorphs ( Figure 6A) at 0 eV was 5.73, 5.02, and 4.33 for HfO 2 , Hf 0.75 Zr 0.25 O 2 , and Hf 0.5 Zr 0.5 O 2 , respectively. It is worth mentioning that the maximum exhibited constants are 8.2, 8.9, and 9.4, respectively, reaching the maximum static dielectric constant at a lower energy range. This indicates that these systems could be easily polarized. The peaks below zero reveal the metallic behavior of the targeted systems at an energy range of 8.4 to 15 eV, which is ascribed to the reflection of the electromagnetic radiation due to the anomalous dispersion approach. On the other hand, the ε 2 was in the order of 10 −4 for all m-phases at 0 eV, showing differences between them up to 4 eV ( Figure 6B); this is in agreement with other studies [23]. ized. The peaks below zero reveal the metallic behavior of the targeted systems at an energy range of 8.4 to 15 eV, which is ascribed to the reflection of the electromagnetic radiation due to the anomalous dispersion approach. On the other hand, the ε2 was in the order of 10 −4 for all m-phases at 0 eV, showing differences between them up to 4 eV ( Figure 6B); this is in agreement with other studies [23]. The o-phase of HfxZryO2 shows a different behavior due to its ferroelectric character that influences the response of the material to the applied field. In detail, the calculated ε1 has different values depending on the Zr amount, moving from 5.88 (without Zr) to 36.06 (in Hf0.5Zr0.5O2) at 0 eV ( Figure 6C).
These values are higher compared to those of the m-structures and are in line with the experimental studies that remark the high values of dielectric constant for the ferroelectric phase [30]. The calculated ε2 remained in the order of 10 −3 up to 3.36 eV ( Figure 6D). The optical transition of the undoped o-polymorph situated in 6.1 eV can be attributed to the intrinsic transition between O 2 p states (at valence band) and Hf 5 d states (at conduction band). More probably, this change can be ascribed the incorporation of Zr and the shift in the transition to a higher energy range, from 6.1 to 8.2 eV, revealing the bandgap increase in doped polymorphs. To conclude, Zr does not have marked effects on the optical properties of m-phases, while it appeared mandatory to modulate the dielectric properties of the ferroelectric o-polymorphs. The o-phase of Hf x Zr y O 2 shows a different behavior due to its ferroelectric character that influences the response of the material to the applied field. In detail, the calculated ε 1 has different values depending on the Zr amount, moving from 5.88 (without Zr) to 36.06 (in Hf 0.5 Zr 0.5 O 2 ) at 0 eV ( Figure 6C).

Oxygen Vacancies Effects
These values are higher compared to those of the m-structures and are in line with the experimental studies that remark the high values of dielectric constant for the ferroelectric phase [30]. The calculated ε 2 remained in the order of 10 −3 up to 3.36 eV ( Figure 6D). The optical transition of the undoped o-polymorph situated in 6.1 eV can be attributed to the intrinsic transition between O 2 p states (at valence band) and Hf 5 d states (at conduction band). More probably, this change can be ascribed the incorporation of Zr and the shift in the transition to a higher energy range, from 6.1 to 8.2 eV, revealing the bandgap increase in doped polymorphs. To conclude, Zr does not have marked effects on the optical properties of m-phases, while it appeared mandatory to modulate the dielectric properties of the ferroelectric o-polymorphs.

Oxygen Vacancies Effects
To better describe the role of Zr on properties of different Hf x Zr y O 2 phases, oxygen vacancies were included in the calculation, putting two O entities in ghost configuration. This approach consists of leaving the basis functions associated with an atom when removing it to create the vacancy. Since each unit cell is constructed of twelve atoms, and eight of them are oxygen entities, the explicit assumption of two considered vacancies leads to an important defect of materials. To quantify the effect of vacancies, the defect formation energy was computed for each polymorph. The calculation of defect formation energy includes (i) the vibrational (phonon) contributions because the creation of a defect modifies the chemical bonds and, thus, the bond strength; (ii) the electronic contributions, which are commonly small for semiconductors but can be sizable for metals; and (iii) the magnetic excitations. Defect-formation energies were identified (Table 3) by computing the total energy formation of each system with and without defects. The values obtained were always positive, and this is not a surprise because, in the case of negative values of the defect formation energy, the not-defected crystals would be unstable. Moreover, polymorphs exhibited different aptitudes to suffer defects, and Zr had an important role in modifying the tendency of different systems to undergo oxygen vacancies. The monoclinic phase always shows a major propensity to suffer oxygen vacancies, and the inclusion of Zr seemed to slightly decrease the formation energy of oxygen vacancies. The reason could be related to the different strength of chemical bonds of Zr with respect to Hf, since the first has a less energetic electronic level than the second. On the other hand, the orthorhombic phase was found to be less likely to undergo oxygen vacancies, and the reason could be ascribed to the ferroelectric behavior of this polymorph. The ferroelectricity, which was correctly reproduced by the simulations, could be at the basis of this energetic difference with the monoclinic phase since the oxygen entities that were placed in the ghost configuration were directly involved in the spontaneous polarization of the phase. Moreover, the inclusion of Zr strongly increased the energy of the oxygen vacancies' formation, and this is probably due to the internal stress modulated by Zr, which finds the best values with 50% of metal. This result is extremely important to identify the relative stability of the ferroelectric Hf x Zr y O 2 polymorphs and to predict the best conditions in the vision of devices' development. The effect of oxygen vacancies on the electrical properties of the polymorphs was also investigated. With the aim to remark upon the differences, the band structures with and without defects were calculated by using the previously described GGA/PBE + U approach. Since no differences were found when different amounts of Zr were included, only the results for Hf 0.5 Zr 0.5 O 2 polymorphs (with and without vacancies) are reported (Figure 7). An increase in the numbers of valence and conduction bands was found in the defective systems of both polymorphs, and the calculated bandgap was always lower than that computed for the pure HfxZryO2 phases. Another evident effect of oxygen vacancies was that the bandgap underwent a shrinkage, and this is directly related to the delocalized states in the valence band, in which the main contribution originated from O 2p orbitals.

Conclusions
HfO2 was found with different polymorphs that confer unique properties to the materials; in the present study, the monoclinic and orthorhombic phases were investigated. In addition, Zr can be used as a dopant agent in different percentages influencing the properties of the systems; thus, it can tune and modulate the devices performances. Although different theoretical and experimental techniques have been already implemented to study the chemical-physical features of these two polymorphs, still a systematic computations investigation of the efficacy of Zr on HfO2 material is missing. In this work, the capabilities of the atomistic simulation methods were discussed, shedding light on the possibility to extrapolate different information. With this aim, six different models were generated: for each P21/c and Pca21 polymorph, HfO2, Hf0.75Zr0.25O2, and Hf0.5Zr0.5O2 were considered. An extensive comparison of different properties was performed by using the its properties, highlighting the negligible effect of the Zr defect. Moreover, no internal stress was observed when Zr entities replaced Hf atoms. A different trend was detected for the Pca2 1 polymorph, in which an important decrease along the z-axis was identified in Hf 0.5 Zr 0.5 O 2 , with important changes in the bond angles. Such a variation stabilized the phase with a significant decrease in terms of internal stress values.
To better describe the electrical and optical properties and the responses to the oxygen vacancies, the GGA + U method was settled and calibrated by using U Hf d = 6 eV, U Zr d = 5.8 eV, and U O p = 4 eV. The PDOS calculation showed that impurity states were located between the CBM and the Fermi level, and this behavior was more evident for o-polymorph. Moreover, following the Mulliken population analysis, the charge-transfer calculations to O atoms indicated completely different results depending on the phases. For the m-phases, Zr doping induced a gradual decrease of charge transfer, while the opposite trend was observed for o-systems. Other important differences were detected by computing the dielectric functions. Minor variations were observed between HfO 2 , Hf 0.75 Zr 0.25 O 2 , and Hf 0.5 Zr 0.5 O 2 m-polymorphs; on the contrary, distinctive ε 1 values without applied field were detected for all the o-polymorphs under study, and variations in features for ε 1 and ε 2 were identified in the high energy range. Finally, the role of oxygen vacancies on materials' electrical properties was evaluated, showing a narrowing of the bandgap regardless of the polymorph type. The m-phases showed a higher propensity to suffer oxygen vacancies, and this behavior gradually decreases when the Zr amount has increased. On the other hand, the o-polymorph resulted in being less prone to these defects, and Zr acts as a phase-stabilizing agent and, at the same time, increases the propensity to suffer vacancies.
All of these results provide important information about the physical-chemical properties of Hf x Zr y O 2 materials. The combined use of GGA/PBE and GGA/PBE + U methods can provide important guidelines in predicting material properties and guiding the development of HfZrO 2 -based devices.