Origin of Interstitial Doping Induced Coercive Field Reduction in Ferroelectric Hafnia

Ferroelectric hafnia (HfO₂) has emerged as a promising candidate for integrating ferroelectric functionalities into integrated circuits, enabled by its robust nanoscale ferroelectricity and exceptional silicon compatibility [1, 2, 3]. However, the high coercive field ($\mathcal{E}_{\rm c}$) required to switch the polarization in HfO₂-based thin films remains a critical issue, impeding the commercialization of this silicon-compatible ferroelectric [3, 4]. Typically, polycrystalline thin films of hafnia fabricated through atomic-layer deposition display an $\mathcal{E}_{\rm c}$ exceeding 1 MV/cm [5]. High-quality epitaxial thin films, obtained by pulsed laser deposition, can show even higher $\mathcal{E}_{\rm c}$ values of $\approx$ 2–5 MV/cm [6, 7, 8, 9]. The need for applying high electric fields, close to the material’s breakdown strength, for polarization reversal seriously limits its field cycling endurance. Reducing $\mathcal{E}_{\rm c}$ without sacrificing the nanoscale ferroelectricity of hafnia is a pressing challenge for the wider adoption of HfO₂-based ferroelectrics.

Various strategies such as doping [10, 11], strain [12, 13], and superlattices [14] have been theoretically explored to lower the barrier for polarization switching, aiming to reduce $\mathcal{E}_{\rm c}$. Investigations based on density functional theory (DFT) revealed that the polarization switching process in the ferroelectric othorhombic ($O$) $Pca2_{1}$ phase often involves an intermediate tetragonal ($T$) $P4_{2}/nmc$ phase, and the energy difference between these two phases could serve as a measure of the switching barrier [11, 15]. Substitutional Si has been suggested as an effective dopant for reducing $\mathcal{E}_{\rm c}$ due to its intrinsic $sp^{3}$ bonding with oxygen, which helps stabilize the intermediate $T$ phase [11, 16]. Conversely, most other substitutional dopants have shown limited effectiveness in influencing $\mathcal{E}_{\rm c}$ based on DFT calculations [11]. Recently, a substantial reduction in $\mathcal{E}_{\rm c}$, down to $\approx$0.65 MV/cm, has been achieved in thin films of Hf(Zr)_1+xO₂, a composition rich in hafnium-zirconium [Hf(Zr)] [17]. This reduction was attributed to the interstitial doping of Hf and Zr atoms into a polar rhombohedral ($R$) $R3m$ phase. DFT calculations show that a 7%-strained $R$-phase Hf_1.08O₂ exhibits a switching barrier of 7.6 meV/atom, significantly lower than that of 20 meV/atom in stoichiometric $R$-HfO₂.

However, assuming a polar $R$ phase in Hf_1.08O₂ films grown via magnetron sputtering leads to notable discrepancies between theory and experiment [18, 19]. The unstrained $R$ phase is actually a nonpolar cubic $P\bar{4}3m$ phase [3, 20]; the polar $R3m$ symmetry only emerges under an equibiaxial compressive strain within the (111) crystallographic plane. Although the relative stability between the $R$ and $O$ phases can be reversed with 8% interstitial doping [17], highly-doped $R$-Hf_1.08O₂ without strain remains nonpolar; a giant compressive strain of 7% is nevertheless required to induce a polarization of $\approx$15 $\mu$C/cm², still below the experimental value of 22 $\mu$C/cm² in polycrystalline films [17]. Moreover, X-ray diffraction (XRD) patterns of Hf(Zr)_1+xO₂ thin films show only a slight offset from both unstrained $R$ and $O$ phases [21], indicating the absence of high strains.

In this work, we combine DFT and deep potential molecular dynamics (DPMD) simulations [22] to investigate how interstitial Hf doping reduces the coercive field in ferroelectric HfO₂. We find that interstitial Hf dopants diminish the energy difference between the polar orthorhombic $Pca2_{1}$ phase and the intermediate tetragonal $P4_{2}/nmc$ phase, lowering the switching barrier and coercive field. Unlike the $R$ phase, which demands additional compressive strains for polarization, the $Pca2_{1}$ phase of Hf_1.03O₂, with a moderate doping concentration of 3%, already exhibits a low switching barrier and a high polarization. We propose that the $Pca2_{1}$ phase of HfO₂ with interstitial Hf defects better explains experimental observations. Large-scale MD simulations at finite temperatures with a DFT-derived machine learning force field demonstrate that 0.5% interstitial doping in $Pca2_{1}$ HfO₂ is enough to reduce the switching field significantly. Moreover, after a pre-poling treatment during which interstitial dopants induce the formation of mobile domain walls, the switching field can drop below 1 MV/cm, comparable to experimental results. Further DFT calculations covering various interstitial dopants reveal a negative correlation between the switching barrier and dopant size.

All DFT calculations are carried out using the Vienna ab initio simulation package (VASP) [23] with the projector augmented-wave (PAW) method [24, 25] and the Perdew-Burke-Ernzerhof (PBE) exchange correlation functional [26]. The cut-off energy for the plane-wave basis is set to 600 eV. The Brillouin zones of the 12-atom pseudocubic unit cell and the 2$\times$2$\times$2 supercell for doped systems are sampled by $\Gamma$-centered 4$\times$4$\times$4 and 2$\times$2$\times$2 Monkhorst-Pack [27] $k$-point meshes, respectively. All the structures are fully optimized until the atomic forces converge to 0.01 eV/Å. The polarization values are determined using the Berry phase method [28, 29]. The polarization switching pathway is determined using the nudged elastic band (NEB) method [30]. The bond order analysis is performed with the crystal orbital Hamilton population (COHP) method [31] implemented in LOBSTER [32]. To investigate the ferroelectric switching in Hf_1+xO₂ at finite temperatures, we perform isobaric-isothermal ensemble ($NPT$) MD simulations using a deep neural network-based force field [33, 34], with model accuracy verification and sample inputs available via an online notebook [35]. The electric fields are included in MD simulations using the “force method” [36, 37, 15], where an additional force $\mathcal{F}_{i}$ is added on the ion $i$ according to $\mathcal{F}_{i}=Z_{i}^{*}\cdot\mathcal{E}$, with $Z_{i}^{*}$ being the Born effective charge (BEC) tensor of ion $i$.

We start with an analysis of the phase stability of Hf_1+xO₂ with interstitial Hf dopants (Hf_i), considering the monoclinic ($M$) $P2_{1}/c$, polar orthorhombic ($O$) $Pca2_{1}$, tetragonal ($T$) $P4_{2}/nmc$, rhombohedral ($R$) $R3m$, and cubic ($C$) $Fm\bar{3}m$ phases. Two doping concentrations, 3.125% and 6.25% (abbreviated as 3% and 6%), are modeled by intercalating one and two excess Hf atom within the hollow sites in a 96-atom 2$\times$2$\times$2 supercell, respectively. As shown in Fig. 1(a), pristine HfO₂ polymorphs exhibit an inverse correlation between the thermodynamic stability and crystal symmetry.With increasing Hf_i concentration, the energies of these phases become closer, suggesting higher Hf_i formation energy in lower-symmetry $M$ and $O$ phases but lower in higher-symmetry $C$ and $R$ phases. This dependence is due to different effective void sizes in each phase [see the inset of Fig. 1(a)]. Specifically, the $M$ and $O$ phases have smaller voids due to alternating fourfold-coordinated (4C) and threefold-coordinated (3C) oxygen atoms, hindering Hf_i integration and resulting in higher formation energy.

Consistent with previous DFT results [17], our calculations show that, compared to the $T$ phase, Hf_i doping destabilizes the lower-symmetry $M$ and $O$ phases while imposing minimal impact on the relative stability of the $R$ phase. This suggests that the $R$ phase would become energetically more favorable than the $O$ phase above a critical Hf_i concentration, if it remains stable. However, the Hf_i-doped $R$ phase is nonpolar, requiring a (111) in-plane compressive strain to achieve a polarized $R$ phase. Using a 36-atom hexagonal supercell with a 7% in-plane compressive strain and an excess Hf atom, Wang et al. reported a polarized $R$-Hf_1.08O₂ with a low switching barrier of 7.6 meV/atom [17]. Notably, such a high strain state should correspond to a large out-of-plane interplanar spacing, $d_{111}=3.11$ Å. This does not align with the experimental XRD data of Hf_1+xO₂ films, which shows a (111) peak at $2\theta=30.1\degree$, indicating only a 0.1$\degree$ shift compared to the stoichiometric $O$ phase (with $d_{111}=$2.96 Å, see Fig. S1).

In contrast to the $R$ phase, the $O$ phase is inherently polar and does not depend on applied strain for polarization. The ferroelectric switching barrier in the $O$ phase can be estimated from the energy difference between the $T$ phase and $O$ phase ($\Delta E=E^{T}-E^{O}$) [11]. As illustrated in Fig. 1(a), interstitial doping elevates the energy of the $O$ phase relative to the $T$ phase, thereby reducing $\Delta E$. Figure 1(b) plots the NEB energy variation along the switching pathway in the $O$ phase for different doping concentrations. The switching barrier decreases drastically from 28 meV/atom in pristine HfO₂ to 12.3 meV/atom in Hf_1.03O₂, and further drops to 7.5 meV/atom in Hf_1.06O₂.

As summarized in Fig. 1(c), for the $O$ phase, both the switching barrier and polarization decrease with increasing Hf_i concentration. This reduction in polarization is attributed to the diminished local displacements of polar oxygen atoms surrounding Hf_i. Further increasing Hf_i concentration destabilizes the $O$ phase, causing it to transform into a $Pbcn$-like structure during DFT geometry relaxation (see Fig. S2). Importantly, $O$-Hf_1.06O₂, with a lower Hf_i concentration, exhibits a switching barrier (7.5 meV/atom) comparable to that of 7%-strained $R$-Hf_1.08O₂. It also features a higher spontaneous polarization along the [111] direction, matching closely to the experimental value. Although large local strains in polycrystalline thin films cannot be ruled out, the $Pca2_{1}$ phase with moderate interstitial Hf doping serves as a minimal model to explain experimental observations. We note that Hf_i exhibits a negative BEC of $-2e$, indicating that its slight movement induces a change in polarization equivalent to the displacement of two negative elementary charges. This counterintuitive result is likely due to the crowding effect of the confining cage around Hf_i. Nevertheless, since the Hf_i is strongly confined, it does not change the polarization switching mechanism identified with DFT-based NEB calculations.

Limited by computational cost, our DFT investigations focused on small supercells. This raises a legitimate question: how does the reduced energy difference between the $O$ and $T$ phases, seen in zero-Kelvin DFT calculations, manifest in larger supercells at elevated temperatures that better represent the experimental sample conditions? To address this, we perform large-scale MD simulations at finite temperatures for an $O$-Hf_1.005O₂ supercell of 20,772 atoms, considering different Hf_i dopant distributions. The force field is a deep neural network-based model potential [15], which accurately reproduces various properties including the energy variation along the switching pathway in Hf_1+xO₂ (see Fig. S3). All simulations are conducted at 400 K to facilitate the switching process on (sub)nanosecond time scales, thereby mitigating computational cost. We gauge the ease of switching by determining the lowest field strength ($\mathcal{E}_{\rm s}$) that triggers polarization switching within 200 ps in MD simulations. As shown in Fig. 2(a), the switching field is 5.3 MV/cm for pristine single-domain HfO₂. Introducing uniformly distributed 0.5% Hf_i dopants decreases $\mathcal{E}_{\rm s}$ to 4.4 MV/cm. The local enrichment of Hf_i dopants further reduce $\mathcal{E}_{\rm s}$: the cluster-like and two-dimensional (2D) distributions of Hf_i lead to $\mathcal{E}_{\rm s}$ of 2.8 and 3.1 MV/cm, respectively (see dopant distributions in Fig. S4).

We discover that a simple pre-poling treatment can bring $\mathcal{E}_{\rm s}$ to an even lower value of $\approx$0.8 MV/cm. The switching process extracted from MD simulations is presented in Fig. 2(c-f), with the color scheme explained in Fig. 2(b). The starting configuration has 3C oxygen atoms displaced upward (P^-) relative to 4C oxygen atom (NP). The 2D-confined Hf_i dopants are distributed randomly within the $xz$ plane between P^- and NP layers. A main impact of Hf_i is the slight downward displacement of 4C atoms [left inset of Fig. 2(c)]. Upon poling under an electric field of $\mathcal{E}=3.1$ MV/cm, the Hf_i enriched layer facilitates the nucleation of 180$\degree$ domain walls characterized by oppositely polarized 3C atoms, called type-P walls [see Fig. 2(d)]. Locally, the type-P wall resembles a $Pbcn$ unit cell, consistent with DFT calculations showing a spontaneous transformation of the $O$ phase to a $Pbcn$-like structure at high Hf_i concentration (see Fig. S2). We confirm that Hf_i dopants reduce the type-P domain-wall energy (see Fig. S5). Consequently, polarization switching is driven by the motion of a type-P wall, while the other type-P wall is pinned by Hf_i dopants. When the external field is turned off, the mobile type-P wall remains stable. After pre-poling, a much lower electric field of 0.8 MV/cm is sufficient to drive the domain-wall motion again.

Above results show that the value of $\Delta E$, easily accessible through DFT calculations at the unit cell level, can indeed serve as a useful descriptor for $\mathcal{E}_{\rm c}$. To identify more types of dopants for regulating $\mathcal{E}_{\rm c}$ in ferroelectric HfO₂, we perform high-throughput DFT calculations covering diverse dopant atoms, $X=$ Hf, Zr, Ti, La, Ta, Y, Nb, Sn, Ge, Si, Al, Ga, In, Sb, Mg, Ca, Cu, Zn, with a fixed concentration of 3%. As shown in Fig. 3(a), $\Delta E$ between the $T$ and $O$ phases with interstitial dopants ($X_{\rm i}$) is plotted against the dopant atomic radius, revealing a clear negative correlation. Interstitial dopants with large atomic radii, such as Hf, Zr, and Ti, tend to reduce the $\Delta E$, while smaller interstitial atoms such as Si, Ge, and Sn tend to increase $\Delta E$, potentially leading to higher $\mathcal{E}_{\rm c}$.

This negative correlation between $\Delta E$ and dopant size for interstitial doping contrasts sharply with substitutional doping, where a volcano-like dependence (mostly a positive correlation) of $\Delta E$ on the ionic radius of the substitutional dopant was observed [11]. We further examine substitutional doping with atom $X$ ($X_{\rm Hf}$) and plot $\Delta E$ as a function of the ionic radius of $X$ [see Fig.3(b)]. Different from interstitial doping, the $\Delta E$ of HfO₂ under substitutional doping increases with dopant size, consistent with previous DFT results [11]. Notably, only a few small substitutional dopants, like Si_Hf, effectively lower $\Delta E$. This has been attributed to the formation of stable $sp^{3}$ bonds between the substitutional Si dopant and neighboring oxygen atoms in the $T$ phase [11].

To comprehend the diametrically opposing effects of Hf_i and Si_i on $\Delta E$ and $\mathcal{E}_{\rm c}$, we analyze the bonding of interstitial atoms in the $T$ and $O$ phases of HfO₂. As depicted in Fig. 4, the intercalated atom is surrounded by eight oxygen atoms. We quantify the bonding strength of an interstitial dopant based on the bond-valence conservation principle [38, 39]. According to this principle, each atom $i$ has a preferred atomic valence $V_{0,i}$, which is often the nominal oxidation state of the atom (e. g., $+4$ for Hf). The actual atomic valence, $V_{i}$, is obtained by summing the individual bond valence ($V_{ij}$) for bonds between the atom and its neighbors ($j$), $V_{i}=\sum_{j}V_{ij}$. A smaller deviation of $V_{i}$ from $V_{0,i}$ reflects stronger overall bonding strength. The bond valence for each oxygen-dopant pair is calculated as $V_{ij}=\exp{[(R_{\rm o}-R_{ij})/B]}$, where $R_{ij}$ is the bond length, $R_{\rm o}$ and $B$ are Brown’s empirical parameters [40]. The calculated $V_{i}$ for Hf_i is 3.19 in the $T$ phase and 3.02 in the $O$ phase, indicating stronger bonding of Hf_i in the $T$ phase. This is further supported by the projected COHP analysis. As shown in the right panel of Fig. 4, the total integration of COHP (ICOHP) for Hf_i$-$O pairs is $-$20.0 and $-$19.0 eV in the $T$ and $O$ phases, respectively. The higher magnitude of ICOHP in the $T$ phase confirms stronger bonding between Hf_i and oxygen atoms.

The difference between the energy of the $O$ and $T$ phases containing Hf_i ($E^{O}[{\rm Hf_{i}}]$ and $E^{T}[{\rm Hf_{i}}]$) can be directly related to the differences in the formation energy of Hf_i ($E_{f}[{\rm Hf_{i}}]$) as follows:

where $\mu_{\rm Hf}$ is the chemical potential of Hf and $E_{0}$ is the energy of the pristine crystal. The relatively stronger Hf_i$-$O bonding in the $T$ phase suggests a lower formation energy of Hf_i than that in the $O$ phase ($E^{T}_{f}[{\rm Hf_{i}}]<E^{O}_{f}[{\rm Hf_{i}}]$) [41], which is responsible for the reduced $\Delta E$ compared to the undoped value ($\Delta E_{0}$). In the case of Si_i, the pCOHP curve reveals a stronger anti-bonding character, which is corroborated by the small values of $V_{i}$. Importantly, both the magnitudes of ICOHP and $V_{i}$ are larger in the $O$ phase than those in the $T$ phase. This points to a lower formation energy of Si_i in the low-energy $O$ phase, which will further increases the energy difference between $T$ and $O$ phases.

In summary, our investigations combining both zero-Kelvin DFT calculations and large-scale MD simulations at finite temperatures establish the link between interstitial doping and coercive field in ferroelectric $Pca2_{1}$ HfO₂. Unit-cell-level DFT calculations act as a mean-field-like analysis, demonstrating that the coercive field reduction from interstitial Hf doping is due to the lower defect formation energy in the intermediate $P4_{2}/nmc$ phase. MD simulations, better representing experimental conditions, suggest that interstitial Hf dopants likely promote switching by forming mobile $Pbcn$-like domain walls. With pre-poling, a switching field of $<$1.0 MV/cm can drive polarization reversal within subnanoseconds. Finally, high-throughput DFT calculations reveal a negative correlation between the switching barrier and the size of the interstitial dopant. The comprehensive understanding of interstitial doping effects in hafnia offers useful guidelines for optimizing the coercive field in this silicon-compatible ferroelectric oxide.