Phase Competition in HfO 2 with Applied Electric Field from First Principles
PPhase Competition in HfO with Applied Electric Field from First Principles Yubo Qi, and Karin M. Rabe
Department of Physics & Astronomy, Rutgers University,Piscataway, New Jersey 08854, United States (Dated: September 29, 2020)In this work, the results of first-principles density–functional–theory calculations are used toconstruct the energy landscapes of HfO and its Y and Zr substituted derivatives as a functionof symmetry–adapted lattice-mode amplitudes. These complex energy landscapes possess multiplelocal minima, corresponding to the tetragonal, oIII ( P ca ), and oIV ( P mn ) phases. We findthat the energy barrier between the non–polar tetragonal phase and the ferroelectric oIII phasecan be lowered by Y and Zr substitution. In Hf . Zr . O with an ordered cation arrangement, Zrsubstitution makes the oIV phase unstable, and it become an intermediate state in the tetragonalto oIII phase transition. Using these energy landscapes, we interpret the structural transformationsand hysteresis loops computed for electric-field cycles with various choices of field direction. Theimplications of these results for interpreting experimental observations, such as the wake–up andsplit–up effects, are also discussed. These results and analysis deepen our understanding of theorigin of ferroelectricity and field cycling behaviors in HfO –based films, and allow us to proposestrategies for improving their functional properties. I. INTRODUCTION
Hafnia (HfO ) has long been recognized as a high– κ material, valuable in complementary metal–oxide–semiconductor (CMOS) applications [1, 2]. In 2011, Si–doped HfO thin films under mechanically encapsulatedcrystallization were observed to be ferroelectric [3], whichled to a renewal of scientific interest and much recent the-oretical and experimental research [3–9].HfO is a binary oxide with a number of polymorphs.At high temperatures ( T > has ahigh–symmetry F m m cubic fluorite structure, in whicheach Hf atom is located at the center of an oxygencube [Fig.1 (a)]. As the temperature decreases (2773K > T > X − mode pattern [Fig. S1 (a)] to give a tetrag-onal P /nmc structure [Fig.1 (b)] in which four Hf–Obonds become slightly shorter while the other four be-come slightly longer [11]. Below 2073 K, bulk HfO isfound in a P /c monoclinic phase, in which additionaldistortions reduce the coordination number of Hf from8 to 7. However, all these known bulk phases are cen-trosymmetric and thus nonpolar. The ferroelectricity ob-served in HfO thin films is attributed to an orthorhom-bic P ca (denoted oIII) phase with polarization along[001], as has been confirmed by both experimental andtheoretical studies [11–15]. A competing [011]–polarizedorthorhombic P nm (denoted oIV) phase [Fig.1 (d)]has been proposed on the basis of first–principles cal-culations [15], but has not yet been experimentally re-ported. These facts indicate that HfO has a complexenergy landscape with multiple local minima. Investigat-ing the state switching between different local minima,especially switching driven by an electric field, is of greatimportance in understanding the origin of ferroelectricityand the behavior in applied electric fields.In this paper, we report the results of first-principles FIG. 1. (a) Cubic, (b) tetragonal, (c) oIII and (d) oIV phasesof HfO . Along the viewing direction, nearer and farther Hfatoms are represented by dark blue and light blue spheres,respectively. Similarly, nearer and farther O atoms are repre-sented by red and pink spheres, respectively. calculations for HfO and substitutional derivatives (Ydoped HfO and Hf . Zr . O ), with the constructionof their energy landscapes as a function of selectedsymmetry–adapted lattice modes. Specifically, we con-sider uniform electric fields applied in various directionswith respect to the crystal axes and map out the struc-tural evolution and hysteresis loops for selected electric a r X i v : . [ c ond - m a t . m t r l - s c i ] S e p field profiles. We find that Y and Zr substitution mod-ifies the energy landscape of pure HfO . In particular,the substitutions change the minimum barrier path fromthe nonpolar tetragonal to the ferroelectric oIII phaseand between up and down ferroelectric variants. Theseresults lead to a better understanding of wake–up andswitching behavior in these systems, enabling electric–field control of competing phases in HfO and its deriva-tives. II. METHODS
We carry out Density functional theory (DFT) calcu-lations with the
Quantum–espresso [16] package forstructural relaxations under finite electric fields and the
ABINIT [17] package for structure optimizations withfixed lattice modes. Optimized norm–conserving localdensity approximation (LDA) pseudopotentials were gen-erated using the
Opium package [18, 19]. A 4 × × k –point mesh was used to sample theBrillouin zone for the conventional 12–atom cell with cor-responding meshes for the supercells. The plane–wavecutoff energy was 50 Ry [20]. Structural relaxations wereperformed with a force threshold of 5.0 × − Hartreeper Bohr. Relaxed structural parameters for the variousbulk phases of pure HfO are reported in supplementarymaterials (SM) section I. The computed lattice constantsare identical to those in our previous work [11] and arein good agreement with the results in other studies.Here, we consider pure HfO and two substitutionalderivatives: Y doped HfO (Y–HfO ) and Hf . Zr . O ,as shown in Fig. 2. The conventional unit cell of pureHfO contains 12 atoms or 4 formula units [Fig.2 (a)].For Y–HfO , we double this unit cell along the y axisand substitute one of the Hf atoms by an Y atom, corre-sponding to a 12.5 % doping concentration [Fig. 2 (b)].For Hf . Zr . O , we consider a structure with 1:1 layeredcation ordering, as shown in Fig. 2 (c). Previous exper-imental work has shown that this superlattice structurewith alternating Zr and Hf cation layers exhibits elec-trical behavior quite similar to that of random–cationHf . Zr . O [21].We describe the low–symmetry structures in each sys-tem by homogeneously straining them back to the lat-tice of the ideal cubic fluorite reference structure anddecomposing the atomic displacements from the refer-ence structure into the symmetry–adapted lattice modesof the reference structure [11, 12]. The atomic displace-ment patterns of the lattice modes and the method forcalculating the mode amplitudes are described in SM sec-tion II. In TABLE S6, the lattice–mode patterns and am-plitudes relevant to the optimized tetragonal, oIII andoIV phases are listed. The oxygen atom positions in thetetragonal structure are obtained from an anti–polar X − mode, with displacements of oxygen chains along the 4–fold direction, here taken as the x direction, alternat-ing in a 2D checkerboard pattern [Fig. S1(a)]. Starting FIG. 2. High–symmetry cubic structures of the (a) bulk HfO ,(b) Y doped HfO , and (c) Hf . Zr . O unit cells. Hf, Zr, Oand Y atoms are represented by blue, green, red and yellowspheres respectively. As in Fig. 1, the darkness of the colorindicates the spatial distance along the view direction; closeratoms are represented by dark colors, which farther ones arerepresented by light colors. from the tetragonal structure, the polar phases are gen-erated by the 3–fold–degenerate zone–center polar Γ − mode, corresponding to uniform displacement of the Hfatoms relative to the simple cubic oxygen network (Fig.S1(c)). We refer to these modes as Γ x , Γ y and Γ z , withamplitudes u x , u y and u z , respectively. The oIII struc-ture, whose polarization is along the [001] direction, hasonly one non–zero polar mode amplitude u z (cid:54) = 0. In theoIV structure, u y = u z (cid:54) = 0, indicating a [011]–directedpolarization.The polarization can be estimated from the linear de-pendence of polarization P on the amplitude of the polarmode in the high–symmetry reference structure as [22] P i = e Ω Z ∗ u i , (1)where i = x, y, z , e is the electronic charge, Ω is thevolume of the supercell, and Z ∗ is the mode Born effectivecharge (SM section III).There are two relevant components of the anti–polar X +5 mode. The X +5 ,y ( q (cid:107) ˆ x ) mode, with amplitude A y ,describes the displacement of oxygen layers along the y axis, alternating in the x direction, as shown in Fig. S1(d). This amplitude is nonzero both in the oIII and oIVstructures. The X +5 ,z ( q (cid:107) ˆ y ) mode, with amplitude B z ,describes the displacement of oxygen layers along the z axis, alternating in the y direction, as shown in Fig. S1(c). This amplitude is nonzero only in the oIII structure. B z and u z combine to give zero and non–zero oxygendisplacements in alternating layers in the oIII structure.For constructing the energy landscape as a functionof the lattice–mode amplitudes, we generate structureswith selected amplitudes of the lattice modes, and then,by using the ABINIT package, relax the structures withthe constraint that the selected mode amplitudes stayfixed while the amplitudes of other modes compatiblewith the resulting symmetry are optimized (SM sectionIII). In certain cases to be discussed below, when the lat-tice mode amplitude is zero, the structure is relaxed withthe lower symmetry corresponding to a small nonzerovalue of the lattice mode amplitude. This ensures thatthe energy plotted as a function of lattice amplitude iscontinuous through zero amplitude, which would not ingeneral be the case if the zero-amplitude energy werecomputed at a saddle point protected by high symme-try. Then, we fit the calculated energies with Landaupolynomials with the lattice–mode amplitudes as orderparameters, and use the polynomials to generate the en-ergy landscapes.We determine the effects of uniform applied electricfield by relaxing the structure with added forces F i onthe ions: F i = Z ∗ i e · E , where Z ∗ i is the effective chargetensor of ion i and E is the applied electric field (SM sec-tion III). To construct electric–field hysteresis loops, anelectric field cycle is applied to each of the systems consid-ered with the tetragonal structure as the initial state. Ineach cycle, we fix a direction for the electric field and thenincrease the magnitude of the electric field in steps from 0to + E max in steps, then decrease it to − E max , and thenincrease it again to + E max . Here, E max = 13 . crystal cansustain in our DFT calculations. Beyond E max = 13 . E = 1 . E is reduced to 0.5 MV/cm. For each value of E , the optimized structure from the previous step is usedas the starting structure for relaxation. Hysteresis loopsare plotted showing the polar–mode amplitudes of theoptimized structures as a function of electric field.The lattice–mode–amplitude energy landscape can beused to understand the electric–field cycling behaviorsto first order in the field. To see this, we start with agiven nonzero electric field and relax the structure as de-scribed above. Next, the resulting structure is relaxed atzero field, constraining the values of u y and u z , to de-termine one point in the lattice–mode–amplitude energylandscape. We find that the amplitudes of the nonpolarmodes in the structure relaxed under a nonzero electricfield, and in the structure computed with a zero electricfield and constrained values of u y and u z , are the same toa good approximation. Therefore, the effect of nonzero field can be described by adding a term − P · E , where P is the spontaneous polarization [23], to the computedzero-field energy landscape. III. RESULTS
We begin by considering the energy landscape of pureHfO in zero applied field, identifying the local minimathat correspond to competing phases and analyzing thebarriers that separate them. We then present resultsfor electric–field hysteresis loops of pure HfO , focusingon the structural changes and switching behavior as themagnitude of the field changes. Finally, we extend ourdiscussion to the Y–HfO and Hf . Zr . O systems, andshow how Y and Zr substitutions change the energy land-scape and electric–field cycling behaviors. A. Pure HfO The energy landscape of pure HfO as a function of thetwo polar mode amplitudes u y and u z is shown in Fig. 3(a). There are three types of local minima, correspondingto the tetragonal, oIII, and oIV phases. The optimalpath from the tetragonal phase to the oIII phase runsalong the horizontal axis, with the optimal value of u y at each u z being zero. Similarly, the optimal path fromthe tetragonal phase to the oIV phase runs along the line u y = u z .The nature of the barriers between the local minimacan be understood by plotting the energy and the relaxedamplitudes A y and B z of the X +5 modes (Fig. S1) for twostraight–line paths in ( u z , u y ) space, as shown in Fig. 3(b) and Fig. 3 (c). Along the straight line path ( u ,0 + )[Fig. 3 (b)], a cusp at u = 0 .
12 ˚A separates the twolocal minima corresponding to the tetragonal and oIIIphases, with a change in space group symmetry from
Aba
2, for the tetragonal phase with a constrained polardistortion, to
P ca , for the oIII phase. In the Aba A y grows linearly with u ,while B z remains zero [Fig. 3 (b) and Fig. S1 (b)]. Inthe P ca phase, found for u > .
12 ˚A, A y ≈ B z (cid:54) = 0,with jumps in both A y and B z across the cusp. The jumpin B z can be understood by plotting the energy landscapeas a function of u z and B z , with other symmetry–allowedmodes, including A y , allowed to relax [Fig. 4 (a)]. For u z < .
12 ˚A, the energy as a function of B z has a singleminimum at B z =0, while for u z > .
12 ˚A, it is a doublewell [Fig. 4 (b)].To understand the nature of the barrier between thetetragonal phase and the oIV phase, we show the energyand the mode amplitudes A y and B z as a function of u along the straight–line path ( u, u ) connecting the tetrag-onal phase to the oIV phase [Fig. 3 (b)]. Nonzero u lowers the symmetry to P nm . With increasing u , A y grows linearly with u z and B z remains zero. The P nm oIV phase is metastable, with an energy of 31 meV/f.u. FIG. 3. (a) Energy landscape of HfO as a function of u y and u z ; We focus on the two paths indicated by the arrows. (b)Energy, A y , and B z as a function of u along the path b fromthe tetragonal to the oIII phase. Dots are calculated pointsand the solid line is an interpolation. (c) Energy, A y , and B z as a function of u along the path c from the tetragonal to theoIV phase. FIG. 4. (a) Energy landscape of HfO as a function of u z and B z (b) Energy as a function of B z , for u z smaller and largerthan the critical value 0.12 ˚A respectively. above the tetragonal phase. The energy barrier from theoIV phase is only 6 meV/f.u. and is smooth rather thancusp–like.To investigate the electrical field cycling behavior, weapplied two electric field cycles to the tetragonal phase,one with electric field along the [001] direction and theother along [011]. The polar–mode–amplitude hysteresisloops are shown in Fig. 5 (a) and the evolution of thefree energy landscape with electric field is shown in Figs.S4 and S5. For E // [001], a tetragonal to polar oIII phasetransition occurs at E T → P = 11 MV/cm. After the elec-tric field is returned to zero, the system is trapped at thelocal minimum corresponding to the oIII structure. Asthe electric field becomes more negative, the polarizationswitches at E P →− P = −
13 MV/cm. Note that after theinitial transition to the oIII phase, field cycling switchesthe system back and forth between up and down vari-ants as in a conventional ferroelectric. For E // [011], thecritical field for inducing a tetragonal to polar oIV phasetransition is E T → P = 7 MV/cm. After the electric fieldis returned to zero, the system is trapped at the localminimum corresponding to the oIV structure. As theelectric field becomes more negative, this local minimumdisappears and at E P → T = − . E P →− P = − B. Y–HfO The Y–doped HfO × × , with nonzero values for the X − mode.Here, we would like to note that for the phases we study,we find that heterovalent Y doping adds a hole at the topof he valence band without introducing any states in thegap, which is consistent with Ref. [24].The energy of Y–HfO as a function of u z is shownin Fig. 6 (a), where we see three phases separated bycusps, with one local minimum at u z = 0 ˚A and a secondat u z = 0 .
17 ˚A. Relative to pure HfO , both the value of u z at the first cusp and the height of the energy barrierseparating the two local minima are lower in Y–HfO .For values of u z below the first cusp value u z < .
06 ˚A,the energy as a function of u z is very close to that of pureHfO . In both systems, with increasing u z , A y coupleswith u z linearly and B z remains zero. At the first cusp u z = 0 .
06 ˚A, both A y and B z exhibit boosts, indicatinga structural transformation with a change in symmetry.A representative structure with u z = 0 .
08 ˚A is shown inFig. 6 (b). The supercell can be divided into two partswith distinct behavior. The unsubstituted part of thesupercell has the atomic displacements characteristic ofthe polar oIII structure, which is responsible for the fer-roelectricity. Here, the adjacent oxygen layers (indicatedby the arrows) have different bonding environments. Asa result, they are prone to have different displacements,leading to a non–zero B z , which is the character of theoIII structure. The part of the supercell containing the Yatom is less distorted. At the second cusp u z = 0 .
10 ˚A,there is a downward jump in B z , indicating another phasetransition. Fig. 6 (b) shows the structure of the polarphase at u = u min = 0 .
18 ˚A. The structure in the unsub-stituted part of the supercell closely resembles the polaroIII structure, while the part of the supercell containingthe Y atom is close to the oIV phase. The decrease in B z is attributed to the formation of the oIV phase in theY-substituted part, with zero B z as indicated in Fig. 3(c).Relative to pure HfO , both the value of u z at thecusp and the height of the energy barrier in Y–HfO arelower. indicating a smaller critical electric field E T → P for triggering the tetragonal to polar phase transition. The behavior in applied electric field exhibits correspond-ing differences. The detailed evolution of the free energylandscape with electric field is shown in Fig. S6. Here,because of the added hole, the density of states at theFermi level is nonzero and we investigate the change ofstructure under an electric field with an accompanyingnonzero field–induced current. As shown in Fig. 5 (b),for E // [001], the critical field E T → P = 7 . E T → P = 11 MV/cm]in pure HfO , as expected from the lower barrier. Futher-more, the saturation polarization is smaller than in pureHfO .Fig. 7 shows the energy profile of Y–HfO as a func-tion of u = u y = u z . For small u , A y couples with u z linearly while B z remains approximately zero, very closeto the behavior in pure HfO . Similar to the energy pro-file of Y–HfO along the [001] direction, Y doping alsointroduces intermediate states, here two rather than one,and the phase changes along the transformation path areindicated by the singular points of the B z profile (seeSM Fig. S3 for the representative structure of each in-termediate state). Even though the relative energy ofthe oIV structure is lowered in Y–HfO , we see that thecurvatures of the local minima, which correlate to themaximum slopes in the energy profile, are approximatelythe same as those in pure HfO , indicating that the crit-ical fields for inducing the phase transitions are similar[Fig. 5 (a) and (b)]. C. Hf . Zr . O The energy landscape of Hf . Zr . O as a function ofthe two polar mode amplitudes u y and u z is shown inFig. 8 (a). As in pure HfO , there are local minima cor-responding to the tetragonal and oIII phases. However,the cation layering normal to the y direction breaks the y – z symmetry of pure HfO , so the two local minima foroIII, with polarization in the y direction and polariza-tion in the z direction respectively, are not equivalent.In addition, the polarization of z–polarized oIII polar-ized in the z direction has a small nonzero y component,which resulting in two symmetry–related variants, onewith positive y component and one with negative y com-ponent. We also see that in Hf . Zr . O , the oIV phaseis unstable. Structural optimization starting from an oIVstructure leads to the oIII minimum, as shown in Fig. 8(a).In Fig. 8 (b), we plot the energy of Hf . Zr . O asa function of u z , with the amplitudes of other modes attheir optimized values, shown in Fig. 8 (a) and (b). Asin pure HfO , this energy profile also exhibits two localminima corresponding to the tetragonal and oIII phases,very little changed from the pure case. However, thecusp seen in pure HfO is cut off, substantially loweringthe energy barrier between the two phases. The lower-ing of the energy barrier results from the change in the FIG. 5. Hysteresis loops of (a) pure HfO , (b) Y–HfO and (c) Hf . Zr . O , with electric field along [001] and [011] directionsrespectively. The evolution of the free energy landscape with electric field is shown in SM section V. In pure HfO , the criticalfield for polarization flipping ( E P →− P = −
13 MV/cm) is very close to the maximum electric field which the crystal can sustain( E max = 13 . // [011], u represents the magnitude of the polar modes projected on the directionof the electric field. The color of the segment (red, blue and green) indicates the state is in the local minimum correspondingto the tetragonal, oIII, oIV structure respectively. For Y–HfO , the blue segment represents that the structure containing anoIII part. path through the energy landscape, which goes throughthe oIV phase region, and around the peak between thetetragonal and oIII phase [see Fig. 8 (a)]. This is in con-trast to pure HfO , in which the optimal path runs alongthe horizontal axis.These features significantly influence the electric–fieldcycling behavior of Hf . Zr . O . As shown in Fig. 5 (c),both [001] and [011] directed electric fields induce a tran-sition from the tetragonal phase to the z –polarized oIIIphase, with the critical field E T → P being much smaller if E is applied along the [011] direction (3.5 MV/cm vs. E// [011]exhibits a complex and unconventional shape. Once an[011]–directed electric field has induced the tetragonalto z –polarized oIII phase transition, the system staystrapped in the local oIII minimum as the field is increasedto E = 13 . E = − . − z –direction. In Fig. S8, we show this in more de-tail, connecting the complexity in the hysteresis loop ofHf . Zr . O to the specific features of its energy land-scape.The surprising fact that the critical field for switch-ing to the oIII phase is lower along [011] than along thepolarization direction [001] is related to the complexityof the energy landscape. Indeed, an electric field along [01¯1] is also more effective in driving the state out of theoIII minimum, compared with a [00¯1] or [0¯1¯1] orientedone. In Fig 8(c), we show the polar mode amplitude hys-teresis loop for an electric field cycle in which the field isfirst increased from zero along [011], decreased to zero,and then increased along [01¯1]. Because the energy land-scape is the same under u z → − u z , the hysteresis loop issymmetrical. The magnitude of the coercive field for flip-ping the polarization is 6.5 MV/cm, smaller than those (7and 12.5 MV/cm) with electric field along [001] and [011]directions. These results open possibilities for improvingthe electrical performance of HfO by modulating the di-rection of applied electric field.We note that the results we have presented are basedon calculations with the LDA functional. Using the gen-eralized gradient approximation (GGA) changes the rela-tive energy of each phase (see SM section I, TABLE S4),but does not destabilize any of the phases we considered.Therefore, the energy landscapes generated by GGA pos-sess the same multiple local minima, with an electric fieldsimilarly inducing phase transitions between them. IV. DISCUSSION
We have demonstrated that the tetragonal phase ofHfO can transform to the ferroelectric oIII or oIV phaseunder an electric field, providing insights into the emer-gence of ferroelectricity and electric-field-cycling behav-iors. In the following subsections, we will relate our com-putational results to experimental observations. FIG. 6. (a) Energy of Y–HfO as a function of u z with otherdegrees of freedom relaxed. The same profile for pure HfO isalso plotted (dashed line) for comparison. (b) Crystal struc-ture of Y–HfO for u z =0.08 and 0.18 ˚A. The structures of theoIII and oIV phase in pure HfO are also shown for compari-son.FIG. 7. Energy, A y , and B z as a function of u z along thepath of the tetragonal to oIV phase transition. Along thistransformation path, u y is fixed to be equal to u z . The rep-resentative structures of each intermediate states are shownin Fig. S3. A. Wake–up and ferroelectricity
FIG. 8. (a) Energy landscape of Hf . Zr . O with respectto the amplitudes of the u y and u z polar modes; (b) Energyprofiles of Hf . Zr . O as a function of u z with energy min-imized with respect to u y . The values of u y at each u z areshown by the black curve with a white edge in (a). The valuesof A y and B z at each u z are shown in the lower panel. (c)Hysteresis loop in which the forward electric is along the [011]direction, and the backward one is along [01-1]. The wake–up effect, which refers to the fact that an as–grown HfO –based thin film needs several electric fieldcycles to establish its full value of switching polarization,has been widely observed in experiments [25–28]. Thiswake–up effect has been attributed to various factors,including structural change at the interface [26, 27, 29],and the internal bias fields caused by oxygen vacancies[30–32].In addition, an electric–field–induced nonpolar to po-lar phase transition is expected to play a significantrole [26, 33–35]. In experiments, doped HfO films areusually deposited at high temperature [3–5, 9]. High tem-perature, surface energy and dopants all promote the for-mation of the tetragonal phase [24, 36]. Upon cooling,the thin films can be kinetically trapped in the tetrago-nal phase, since it is metastable, with no unstable phononmodes [11]. However, the oIII structure has a lower en-ergy than the nonpolar tetragonal phase, so that once themore stable oIII structure is induced by an electric field, aelectric field opposite to the polarization will tend to flipthe polarization to the corresponding oIII phase ratherthan drive the system back to the tetragonal phase.In our calculations, for pure HfO the critical elec-tric field E T → P for triggering this phase transition is 11MV/cm. Our calculations also show that Y and Zr dop-ing can cut off the cusp in the energy profiles by introduc-ing intermediate states, and thus lower the energy barrierand critical fields ( E T → P = 7 . and4 MV/cm in Hf . Zr . O ). Moreover, we note that sincethe energy barriers are not very high, thermal vibrationcan facilitate the phase transition and further lower thecritical fields. With the simple assumption that a phasetransition will be driven by thermal vibration if (1) com-pared with the current state, the energy of another stateis lower by more than k B T / k B T / ) and 2.2MV/cm (in Hf . Zr . O ), which are the same order ofmagnitude as the experimental values [3, 4, 9]. B. Multiple jumps in the P – E loop In a conventional ferroelectric hysteresis loop, the po-larization jumps at the coercive field. However, in HfO based thin films, jumps in polarization can occur at twoor more values of the electric field, leading to an irregularshape of the P – E loop (as in Fig. 5). This phenomenonhas been discussed in Ref. [31], and referred as the ‘split– up’ effect. This ‘split–up’ effect has been attributed tothe different activation energies in grains with differentoxygen vacancy concentrations [31]. Here, we would liketo point out that the competing phases in the energylandscape can also contribute to this effect. As shown inFig. 5, state hopping between different (tetragonal, oIIIand oIV) local minima can lead to multiple boosts. More-over, our results demonstrate that electric field cyclingbehaviors are also field–direction dependent. In HfO -based thin films composed of multiple grains [27, 28, 37],grains with different lattice orientations may follow dif-ferent transformation paths with different critical fields. V. CONCLUSION
In this study, we investigated the energy landscapes ofpure HfO , Y–doped HfO , and Hf . Zr . O as a func-tion of the polar modes u y and u z . The complex energylandscapes of these systems are found to possess multiplelocal minima, corresponding to the tetragonal, oIII, andoIV phases respectively. Moreover, we find that Y and Zrdoping can lower the energy barriers between the non–polar tetragonal phase and the ferroelectric oIII phase, byintroducing intermediate states. In Hf . Zr . O with anordered cation arrangement, the oIV phase becomes un-stable, and a new transformation path for the tetragonalto oIII phase transition is opened, with a lower energybarrier. We connect our results to experimental obser-vations, such as the wake–up effect and P – E loop withirregular shape, and also propose dependence of behav-ior on the direction of applied electric field. This workprovides useful insights about the origin of the ferroelec-tric phase in HfO –based films, and suggests strategiesfor improving their ferroelectric performance. ACKNOWLEDGMENTS
We thank Sebastian E. Reyes–Lillo, Fei–Ting Huang,Xianghan Xu, and Sang–Wook Cheong for valuable dis-cussions. This work was supported by Office of NavalResearch Grant N00014–17–1–2770. Computations wereperformed by using the resources provided by the High–Performance Computing Modernization Office of the De-partment of Defense and the Rutgers University ParallelComputing (RUPC) clusters. [1] M. T. Bohr, R. S. Chau, T. Ghani, and K. Mistry, IEEEspectrum , 29 (2007).[2] M. Gutowski, J. E. Jaffe, C.-L. Liu, M. Stoker, R. I.Hegde, R. S. Rai, and P. J. Tobin, Appl. Phys. Lett. ,1897 (2002).[3] T. B¨oscke, J. M¨uller, D. Br¨auhaus, U. Schr¨oder, andU. B¨ottger, Appl. Phys. Lett. , 102903 (2011). [4] T. S. B¨oscke, S. Teichert, D. Br¨auhaus, J. M¨uller,U. Schr¨oder, U. B¨ottger, and T. Mikolajick, Appl. Phys.Lett. , 112904 (2011).[5] S. Mueller, C. Adelmann, A. Singh, S. Van Elshocht,U. Schroeder, and T. Mikolajick, ECS J. Solid State Sci.Technol. , N123 (2012). [6] S. Mueller, J. Mueller, A. Singh, S. Riedel, J. Sundqvist,U. Schroeder, and T. Mikolajick, Adv. Funct. Mater. ,2412 (2012).[7] J. M¨uller, T. B¨oscke, D. Br¨auhaus, U. Schr¨oder,U. B¨ottger, J. Sundqvist, P. K¨ucher, T. Mikolajick, andL. Frey, Appl. Phys. Lett. , 112901 (2011).[8] J. M¨uller, U. Schr¨oder, T. B¨oscke, I. M¨uller, U. B¨ottger,L. Wilde, J. Sundqvist, M. Lemberger, P. K¨ucher,T. Mikolajick, et al. , J. Appl. Phys. , 114113 (2011).[9] J. Muller, T. S. Boscke, U. Schroder, S. Mueller,D. Brauhaus, U. Bottger, L. Frey, and T. Mikolajick,Nano Lett. , 4318 (2012).[10] R. Ruh, H. J. Garrett, R. F. Domagala, and N. M. Tal-lan, J. Am. Ceram. Soc. , 23 (1968).[11] S. E. Reyes-Lillo, K. F. Garrity, and K. M. Rabe, Phys.Rev. B , 140103 (2014).[12] H.-J. Lee, M. Lee, K. Lee, J. Jo, H. Yang, Y. Kim, S. C.Chae, U. Waghmare, and J. H. Lee, Science (2020).[13] X. Sang, E. D. Grimley, T. Schenk, U. Schroeder, andJ. M. LeBeau, Appl. Phys. Lett. , 162905 (2015).[14] M. H. Park, Y. H. Lee, H. J. Kim, Y. J. Kim, T. Moon,K. D. Kim, J. M¨uller, A. Kersch, U. Schroeder, T. Miko-lajick, et al. , Adv. Mater. , 1811 (2015).[15] T. D. Huan, V. Sharma, G. A. Rossetti Jr, and R. Ram-prasad, Phys. Rev. B , 064111 (2014).[16] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car,C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococ-cioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris,G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis,A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari,F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello,L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P.Seitsonen, A. Smogunov, P. Umari, and R. M. Wentz-covitch, J. Phys.: Condens. Matter , 395502 (2009).[17] X. Gonze, J.-M. Beuken, R. Caracas, F. Detraux,M. Fuchs, G.-M. Rignanese, L. Sindic, M. Verstraete,G. Zerah, F. Jollet, M. Torrent, A. Roy, M. Mikami,P. Ghosez., J.-Y. Raty, and D. Allan, Comp. Mater.Sci. , 478 (2002).[18] http://opium.sourceforge.net.[19] J. W. Bennett, Phys. Procedia , 14 (2012).[20] H. J. Monkhorst and J. D. Pack, Phys. Rev. B , 5188(1976).[21] S. L. Weeks, A. Pal, V. K. Narasimhan, K. A. Littau,and T. Chiang, ACS Appl. Mater. Interfaces , 13440(2017).[22] K. M. Rabe, C. H. Ahn, and J.-M. Triscone, Physics ofFerroelectrics: A Modern Perspective (Springer–Verlag, 2007).[23] N. Sai, K. M. Rabe, and D. Vanderbilt, Phys. Rev. B , 104108 (2002).[24] C.-K. Lee, E. Cho, H.-S. Lee, C. S. Hwang, and S. Han,Phys. Rev. B , 012102 (2008).[25] D. Zhou, J. Xu, Q. Li, Y. Guan, F. Cao, X. Dong,J. M¨uller, T. Schenk, and U. Schr¨oder, Appl. Phys. Lett. , 192904 (2013).[26] M. Peˇsi´c, F. P. G. Fengler, L. Larcher, A. Padovani,T. Schenk, E. D. Grimley, X. Sang, J. M. LeBeau, S. Sle-sazeck, U. Schroeder, et al. , Adv. Funct. Mater. , 4601(2016).[27] E. D. Grimley, T. Schenk, X. Sang, M. Peˇsi´c,U. Schroeder, T. Mikolajick, and J. M. LeBeau, Adv.Electron. Mater. , 1600173 (2016).[28] D. Martin, J. M¨uller, T. Schenk, T. M. Arruda, A. Ku-mar, E. Strelcov, E. Yurchuk, S. M¨uller, D. Pohl,U. Schr¨oder, et al. , Adv. Mater. , 8198 (2014).[29] M. Hoffmann, U. Schroeder, T. Schenk, T. Shimizu,H. Funakubo, O. Sakata, D. Pohl, M. Drescher, C. Adel-mann, R. Materlik, et al. , J. Appl. Phys. , 072006(2015).[30] T. Schenk, M. Hoffmann, J. Ocker, M. Pesic, T. Mikola-jick, and U. Schroeder, ACS Appl. Mater. Interfaces ,20224 (2015).[31] T. Schenk, U. Schroeder, M. Pesic, M. Popovici, Y. V.Pershin, and T. Mikolajick, ACS Appl. Mater. Interfaces , 19744 (2014).[32] C. Mart, K. Kuhnel, T. Kampfe, M. Czernohorsky,M. Wiatr, S. Kolodinski, and W. Weinreich, ACS Appl.Mater. Interfaces , 2612 (2019).[33] M. Hoffmann, U. Schroeder, C. K¨unneth, A. Kersch,S. Starschich, U. B¨ottger, and T. Mikolajick, Nano En-ergy , 154 (2015).[34] M. H. Park, H. J. Kim, Y. J. Kim, Y. H. Lee, T. Moon,K. D. Kim, S. D. Hyun, and C. S. Hwang, Appl. Phys.Lett. , 192907 (2015).[35] R. Batra, T. D. Huan, J. L. Jones, G. Rossetti Jr, andR. Ramprasad, J. Phys. Chem. C , 4139 (2017).[36] R. Batra, H. D. Tran, and R. Ramprasad, Appl. Phys.Lett. , 172902 (2016).[37] M. Hoffmann, M. Peˇsi´c, K. Chatterjee, A. I. Khan,S. Salahuddin, S. Slesazeck, U. Schroeder, and T. Miko-lajick, Adv. Funct. Mater.26