Phase Field Modeling of Domain Dynamics and Polarization Accumulation in Ferroelectric HZO

TL;DR

This paper presents a phase-field simulation framework to study polarization switching in ferroelectric HZO films, revealing how domain-wall instability influences accumulative polarization behavior under electric pulses.

Contribution

The work introduces a dynamic phase-field model incorporating domain-wall instability to analyze polarization accumulation in ferroelectric HZO, a novel approach for understanding switching dynamics.

Findings

Domain-wall instability affects polarization switching behavior.

Electric field and pulse timing influence polarization accumulation.

Domain-wall motion is key to understanding P-switching dynamics.

Abstract

In this work, we investigate the accumulative polarization (P) switching characteristics of ferroelectric (FE) thin films under the influence of sequential electric-field pulses. By developing a dynamic phase-field simulation framework based on time-dependent Landau-Ginzburg model, we analyze P excitation and relaxation characteristics in FE. In particular, we show that the domain-wall instability can cause different spontaneous P-excitation/relaxation behaviors that, in turn, can influence P switching dynamics for different pulse sequences. By assuming a local and global distribution of coercive field among the grains of an FE sample, we model the P-accumulation process in Hf0.4Zr0.6O2 (HZO) and its dependency on applied electric field and excitation/relaxation time. According to our analysis, domain-wall motion along with its instability under certain conditions plays a pivotal role in accumulative P-switching and the corresponding excitation and relaxation characteristics.

Figures

Equations

$$-{\rho }_n\frac{\partial P_n}{\partial t}=-K_P^n{\nabla }^2{P}_n-E_n^{app}+\hat{\alpha }{P}_n+\hat{\beta }{P}_n^3+\hat{\gamma }{P}_n^5$$

Full text

††preprint: AIP/123-QED

Ferroelectric (FE) materials, particularly Zr doped HfO2 (Hf1-xZrxO2:HZOMuller et al. (2012)) have drawn significant research interest in recent times due to CMOS process compatibilitySharma et al. (2017), thickness scalabilityWei et al. (2014); Khan et al. (2018) as well as many promising attributes of ferroelectric field effect transistorsSalahuddin and Datta (2008); Sharma et al. (2017) (FEFETs) for low-power logicSalahuddin and Datta (2008); Khan et al. (2014, 2011); Saha et al. (2017) and non-volatile memoriesTrentzsch et al. (2016); Chatterjee et al. (2017) applications. In addition, FEFETs can provide multiple non-volatile resistive states that harness the multi-domain FE characteristics, leading to the possibilities for multi-bit synapsesMulaosmanovic et al. (2017); Jerry et al. (2017) in a neuromorphic hardwareIndiveri and Liu (2015). Further, newly reported accumulative polarization ($P$)-switching processMulaosmanovic et al. (2018a) in ultra-thin FE leads to many appealing opportunities for novel applications like correlation detectionNi et al. (2018) and other non-Boolean computing paradigmsMulaosmanovic et al. (2018b). For such emerging applications of FEFETs, the $P$-switching dynamics in response to sub/super-coercive voltage pulse trains play an important role and are, therefore, critical to understand.

To that effect, this letter analyzes spatially local $P$-switching dynamics and its participation in globally observable $P$-accumulation characteristics in response to a pulse train. Our analysis is based on a dynamic phase field modelLi et al. (2001); Rabe, Ahn, and Triscone (2007) coupled with measured accumulation characteristics of HZO. By providing the spatial distribution of $P$ ($P$-map) in different electric field (E-field) excitation and relaxation steps, we discuss different types of $P$ excitation and relaxation processes and their corresponding dependency on E-field ($E$) amplitude ($E_{max}^{app}$), ON time (or excitation time $T_{on}$) and OFF time of the pulse (or relaxation time $T_{off}$). Finally, considering a coercive-field distribution among different FE grains, we analyze the experimental $P$-accumulation characteristics by utilizing our simulation framework.

Let us start by describing the experimentally observed $P$-switching characteristics in HZO. Fig. 1(a) shows the measured charge vs. E-field ($Q$-$E$) characteristics of a 10nm HZO film (x=0.6, grown by ALD with TiN capping layer as top and bottom contact). Here, $Q$=$P$+$\epsilon_{0}E$, where $\epsilon_{0}$ is vacuum permittivity. We observe excellent accumulative $P$-switching in HZO as the response of successive E-field stimulation (Fig. 1(b)), where the $P$-accumulation ($P^{acc}$) characteristics exhibit a strong dependence on the E-pulse properties. For example, we observe faster $P^{acc}$ with the increase in $E_{max}^{app}$ (Fig. 1(c)), increase in $T_{on}$ (Fig. 1(d)) and/or decrease in $T_{off}$ (Fig. 1(e)). Also, $P^{acc}$ saturates after a certain number of pulses. Such saturation occurs at higher $P$ with the increase in $E_{max}^{app}$, increase in $T_{on}$ and decrease in $T_{off}$. It is noteworthy that such accumulated-$P$ observed in the experiments is the average of locally accumulated-$P$ in different grains. Hence, to explain the experimental results described above, it is critical to understand the spatially local $P$-switching dynamics in an individual grain. We analyze such processes in detail based on our phase-field model, calibrated to the experiment.

Considering a thin ($\sim$10nm) FE film, we assume that the $P$ direction is only along the thickness ($z$-axis). Hence, $P_{x}$=0, $P_{y}$=0, $P_{z}\neq$0 and $P_{z}$ can have a spatial distribution in $x$-$y$ plane. However, we assume uniform $P_{z}$ along the $z$-axis ($dP_{z}/dz$=0) owing to the ultra-thin nature of the film. The temporal and spatial evolution of $P$ can be described by the time ($t$)-dependent Landau-Ginzburg (TDLG) equationRabe, Ahn, and Triscone (2007); Nambu and Sagala (1994): ${\delta F}/{\delta P_{z}}=-\rho({\partial P_{z}}/{\partial t})$. Here, $\rho$ is the kinetic coefficient and $F$ is the total energyNambu and Sagala (1994) of the system composed of free energy, domain-wall energy, electrostatic and electroelastic energy. Considering up to the 6th order terms in Landau’s free energy expansionLynden-Bell (1995), the normalized representation of TDLG equation within the FE is given by the following equation.

[TABLE]

Here, $P_{n}(=P_{z}/P_{C}$) and $E_{n}^{app}(=E^{app}/E_{C}$) are polarization and E-field normalized with respect to $E_{c}$ (coercive field) and $P_{C}$ ($P$ at $E$=$E_{C}$), respectively. $\hat{\alpha}$, $\hat{\beta}$ and $\hat{\gamma}$ are the normalized Landau coefficients, $\rho_{n}$ is the normalized kinetic coefficient and $K_{P}^{n}$ is the normalized domain-interaction parameter. Eqn. 1, is similar to the Euler-Lagrange equation of motionElder, Gould, and Tobochnik (1993), incorporates interactions among neighboring domains (via $K_{P}^{n}\nabla^{2}P_{n}$) giving rise to Klein-Gordon type field equationBandyopadhyay, Ray, and Gopalan (2006). In our simulations, we self-consistently solve eqn. 2 in a real space grid by considering Neumann boundary at the edgesCano and Jiménez (2010). We include a comprehensive discussion on parameter extraction, normalization and simulation methodology in the supplementary section. It is noteworthy that $P_{n}$ denotes normalized microscopic $P$ in each grid point, while the analogous quantity of experimentally measured $P$ is spatial average of $P_{n}$, denoted as $\bar{P}_{n}$. Also, $K_{P}^{n}\nabla^{2}P_{n}$ can be thought of as the local effective interaction E-field, $E_{n}^{int}$. Therefore, the $P$ switching depends on $E_{n}^{app}$+$E_{n}^{int}$. For instance, $P$-switching will occur for $E_{n}^{app}$+$E_{n}^{int}>$1 (since the normalized coercive field=1).

In general, $P$-switching can take place in two different ways, namely (i) direct nucleation and (ii) Domain-wall (DW) assisted nucleation. To understand such $P$-switching processes, let us start by considering the FE sample in Fig. 2(a), where region R1 exhibits $P_{n}=+|P_{n,r}|$ and R2 exhibits $P_{n}=-|P_{n,r}|$ and they are separated by a DW within which $P_{n}$ varies gradually along the $x$-axis (Fig. 2(b)). Here, $P_{n,r}$ is the remnant polarization. In this case, the domain structure is effectively 1D as $d^{2}P_{n}/dy^{2}=0$. Direct nucleation occurs for super-coercive applied fields ($|E_{n}^{app}|>1$), wherein region R2 will switch to $+P$ at once if E-field is applied for a sufficient time. On the other hand, DW assisted nucleation (which is the main focus of this work) is observed for sub-coercive applied fields ($|E_{n}^{app}|<1$), in which $E_{n}^{int}$ plays a key role. To explain this, let us consider $E_{n}^{app}=0$ and static condition ($dP_{n}/dt$=0). Hence, eqn. (1) can be written as, $E_{n}^{int}$=$K_{P}^{n}\nabla^{2}P_{n}$=$\hat{\alpha}_{n}P_{n}$+$\hat{\beta}_{n}P_{n}^{3}$+$\hat{\gamma}_{n}P_{n}^{5}$. Note, $E_{n}^{int}$ is localized and non-zero only within DW (Fig. 2(c)). Fig. 2(d) shows the relation between $E_{n}^{int}$ and $P_{n}$, signifying that the symmetric spatial distribution of $P_{n}$ provides a symmetric $E_{n}^{int}$ in the 1D case. Here, the symmetric $E_{n}^{int}$ plays a critical role by balancing the forces due to $E_{n}^{int}$ on the two sides of the DW, yielding stable and static DW for $E_{n}^{app}$=0. This can also be understood by noting that $E_{n}^{int}\leq$1, which leads to stable DW due to no $P$-switching. However, by applying a sub-coercive E-field, $0<|E_{n}^{app}|<1$, we can get a total local E-field, $|E_{n}^{app}$+$E_{n}^{int}|>$1 and that can eventually initiate $P$-switching. Such $P$-switching is a spatially local and temporally gradual process, which is referred to DW motion or DW assisted nucleation. (For more details on 1D $P$-switching, see supplementary section).

With the understanding of the 1D $P$-switching, let us now consider 2D FE dynamics for which, we analyze two cases by considering homogeneous $E_{C}$ (case-1) and a distribution of $E_{C}$ (case-2) in an FE grain. Let us start with case-1 and consider a square FE sample (Fig. 2(e)) where a circular region R1 exhibits $P_{n}=+|P_{n,r}|$, which is surrounded by R2 with $P_{n}=-|P_{n,r}|$. Here, the DW is 2D and the corresponding $\nabla^{2}P_{n}$ in the polar coordinate can be written as $[(\partial^{2}P_{n}/\partial r^{2})+(1/r)(\partial P_{n}/\partial r)]$. Note that, the $(1/r)(\partial P/\partial r)$ term exhibits a radial dependency and hence, $\nabla^{2}P_{n}$ becomes radially asymmetric, even for a symmetric radial distribution of $P_{n}$ (symmetric $\partial^{2}P_{n}/\partial r^{2}$ and $\partial P_{n}/\partial r$). Therefore, $E_{n}^{int}$ becomes spatially asymmetric (Fig. 3(f)), where $|E_{n}^{int}|>1$ at the inner interface and $|E_{n}^{int}|<1$ at the outer interface of DW. Such asymmetry in the $E_{n}^{int}$ causes the DW to undergo an effective inward force. Consequently, the DW becomes unstable and R1 region shrinks spontaneously with time (as $|E_{n}^{int}|>1$ at the interior end of the DW). Such spontaneous phenomena play an important role in the $P$-switching dynamics that we discuss subsequently. (Note, such a DW instability is in contrast with the 1D case that we discussed above where symmetric $E_{n}^{int}$ leads to stable DW).

Let us now consider a sequence of sub-coercive E-field pulses ($E_{n,max}^{app}$=0.8) applied to this sample of FE. Simulated transient $\bar{P}_{n}$ is shown in Fig. 3(a) and the initial $P$-map at $t$=0ns is shown in Fig. 3(b)-i, where the initially switched region (red) can be assumed as a pinned domain. After the arrival of first E-field pulse, R1 domain grows circularly, nucleating new lattices sequentially at the outer edge of the DW (Fig. 3(b)-ii). That implies an increase in R1 area and decrease in R2 area by an amount $\Delta A_{j}^{+}$ ($j$= E-field pulse number). The corresponding $P$-excitation characteristics (Fig. 3(a):(i-ii)) comprise of three different components (Fig. 3(c)), i.e. type-1: $-|P_{n,r}|$ to $-|P_{n,e1}|$ (in R2), type-2: $+|P_{n,r}|$ to $+|P_{n,e2}|$ (in R1) and type-3: $-|P_{n,r}|$ to $+|P_{n,e2}|$ (in $\Delta A_{j}^{+}$).

After the end of first E-field pulse, the DW propagation stops and the $P$ changes due to type 1 and 2 components get immediately relaxed to $-|P_{n,r}|$ and $+|P_{n,r}|$, respectively. We call these type-1 and 2 relaxations, respectively (Fig. 3(c)). Similarly, the newly nucleated part ($\Delta A_{j}^{+}$) also rapidly get relaxed to $+|P_{n,r}|$ by following type-2 relaxation. Corresponding transient relaxation in $\bar{P}_{n}$ can be seen in Fig. 3(a) (from point ii to ii(a)). Interestingly, followed by such rapid relaxation, there is another relaxation component that gradually reduces $\bar{P}_{n}$ until the arrival of next E-field pulse (Fig. 3(a): ii(a)-iii). Such spontaneous $P$-relaxation is the outcome of DW instability (due to $E_{n}^{int}$ asymmetry – discussed above) that causes spontaneous shrinking of R1 domain (Fig. 3(b): ii-iii). Let us define the spontaneous decrease in R1 area in the absence of E-field as $\Delta A_{j}^{-}$.

Now, due to sequential E-field pulses, rather than completely collapsing, DW moves further towards the grain boundary by following $P$-excitation and relaxation sequences. Once DW reaches sufficiently close to the grain boundary, spontaneous relaxation is not observed in the absence of E-field (during $T_{off}$). Instead, spontaneous $P$-excitation (Fig. 3(a):v(a)-vi) takes place. This is because, as the DW reaches near the edges, R2 domain becomes very narrow and hence, $\nabla^{2}P_{n}$ increases. Therefore, at the outer interface of the DW (alongside R2), $|E_{n}^{int}|>1$ and that causes an effective outward force in R2. As a result, R2 domain spontaneously switches to $+P$ (Fig. 3(b): iv-vi). After all the lattices switch to $+|P_{n,r}|$, transient $\bar{P}_{n}$ exhibits only type-2 excitation and relaxation.

Note that the $P$-maps at the beginning of different E-field pulses are different and play an important role in each $P$-excitation/relaxation behavior. Let us define the increase in $\bar{P}_{n}$ during each excitation period ($|\Delta^{+}\bar{P}_{n}|_{j}$). Recall that each $|\Delta^{+}\bar{P}_{n}|_{j}$ consists of three excitation components. It can be shown mathematically from eqn. 1 that $P$ change due to type-1 excitation is higher in magnitude than type-2 (Fig. 3(c)). Now, with the increase in pulse number ($j$), R1 area increases and R2 area decreases. Therefore, the contribution from type-2 excitation (in R1) increases and type-1 excitation (in R2) decreases. Hence, we expect an overall decrease in total excitation (type-1+2) with respect to (w.r.t.) $j$. Now, type-3 excitation (corresponding to R1 area increase during $j$-th excitation) can be written as, $\Delta A_{j}^{+}$=$\pi[(r_{j}+dr_{j})^{2}-r_{j}^{2}]$=$\pi[dr_{j}^{2}+2r_{j}dr_{j}]$. Here $r_{j}$ is the domain radius before $j$-th excitation and $dr_{j}$ is the increase in radius during $j$-th excitation. Note that a linear increase in $dr_{j}$ w.r.t. time gives rise to quadratic increase in R1 area. That implies that if we keep the E-pulse ON for a long time, the $\bar{P}_{n}$ dynamics will be quadratic w.r.t. time (gray lines in Fig. 3(a)). Similarly, assuming $dr_{j}$ as constant irrespective of the value of $j$, we can see that $\pi\times 2r_{j}dr_{j}$ increases with $j$ as $r_{j}$ increases. That implies an incremental change in $\Delta A_{j}^{+}$ and hence, type-3 excitation component increases with $j$. Note that the type-3 contribution is dominant over type-1+2 and therefore, $|\Delta^{+}\bar{P}_{n}|_{j}$ increase with the increase with $j$ up to $j$=6 (Fig. 3(d)). For, $j>$6, R1 domain reaches the grain boundary and the quadratic growth of R1 no longer holds true. Hence, type-3 contribution decreases significantly, leading to the domination of type-1+2 excitation and decrease in $|\Delta^{+}\bar{P}_{n}|_{j}$ with $j$. After R2 domain vanishes (or switched to $+P$ at $j$=9), $|\Delta^{+}\bar{P}_{n}|_{j}$ only exhibits type-2 excitation.

Similarly, the decrease in $\bar{P}_{n}$ during each relaxation period ($|\Delta^{-}\bar{P}_{n}|_{j}$) consists of type-1+2 relaxation and a spontaneous part. Like type-1+2 excitation, type-1+2 relaxation decreases as $j$ increases. However, the spontaneous component ($\Delta A_{j}^{-}$) behaves non-monotonically w.r.t. $j$. $\Delta A_{j}^{-}$ (decrease in R1 area) changes sign from positive ($+$) to negative ($-$) at $j$=6 as the spontaneous component changes from relaxation to excitation characteristics. Therefore, till the spontaneous contribution is relaxation (j$$\leq$$6), $|\Delta^{-}\bar{P}_{n}|_{j}$ decreases with the increase in $j$. Once the spontaneous contribution leads to excitation (j$$>$$7), $|\Delta^{-}\bar{P}_{n}|_{j}$ increases with $j$ and becomes constant at $j$=9 (with only type-2 relaxation in R1).

$\bar{P}_{n}$ at the end of each excitation-relaxation sequence, called accumulated polarization ($P_{n}^{acc}$) is shown in Fig. 3(a). Note that the change in $P_{n}^{acc}$ at each pulse is basically proportional to $\Delta A_{j}^{+}$-$\Delta A_{j}^{-}$. We discussed earlier that $\Delta A_{j}^{+}$ shows incremental increase with the increase in $j$, whereas, $\Delta A_{j}^{-}$ exhibits non-monotonic change along with a sign change from ‘+’ to ‘-’. Therefore, $P_{n}^{acc}$ initially increases slowly when $\Delta A_{j}^{-}$ is ‘+’ and once $\Delta A_{j}^{-}$ becomes ‘-’, then increases rapidly. On the other hand, the flat region (Fig. 3(a), $j\geq$9) in $P_{n}^{acc}$ signifies an absence in $P$-accumulation once the whole sample (or grain) completely switches to $+P$.

The trends in $P$ accumulation w.r.t the pulse attributes are illustrated in Fig. 3(e-g). With the increase in pulse amplitude ($E_{n,max}^{app}$), R1 domain grows more rapidly (increase in $dA^{+}/dt$) leading to faster accumulation (Fig. 3(e)). Similarly, with the increase in $T_{on}$, $\Delta A_{j}^{+}$ increases during each pulse and therefore, $P_{n}^{acc}$ saturates at a lower $j$ (Fig. 3(f)). Also, an increase in $T_{off}$ leads to an increase in spontaneous relaxation (increase in $\Delta A_{j}^{-}$). Consequently, larger number of pulses is required for $P_{n}^{acc}$ to get saturated (Fig. 3(g)). Note that, if $E_{n,max}^{app}$ or/and $T_{on}$ is/are very low, so that $(\Delta A_{j}^{+}$-$\Delta A_{j}^{-})<0$, then, rather than accumulation, R1 can get completely relaxed to $-P$. The same is true for a high $T_{off}$. Such scenarios can be seen in Fig. 3(e-g) (gray lines). Note that here, we assume the pinned-type (or initially nucleated) domain at the center of the grain. However, such pinned domain can be at random positions within the grain. In such cases too, the trends in $P$-switching dynamics remain the same as discussed (see supplementary section).

Now, we consider case-2, where we assume a Gaussian distribution of $E_{C}$ in an FE grain (Fig.4(a)) by considering a spatial distribution of $\hat{\alpha}$, $\hat{\beta}$ and $\hat{\gamma}$ (see supplementary section). Note that $E_{C}$ is assumed to be less near the grain boundary (Fig. 4(b)), which can be understood as the cause of strain relaxation near the edgesEmelyanov et al. (2002); Lee et al. (1999). Like the previous discussion, considering a sequence of E-field pulses ($E_{n,max}^{app}$=0.92), simulated $\bar{P}_{n}$ is shown in Fig. 4(c). Note that the FE grain was initially switched to $-|P_{n,r}|$ (Fig. 4(d): i). Therefore, $P$-switching occurs as a two-step process: (1) E-field induced nucleation and (2) E-field assisted domain growth. Once the first E-field pulse arrives, direct nucleation starts from the grain edges (with lower $E_{C}$) and propagates inward (Fig. 4(d): ii). After the end of E-field pulse, further nucleation stops and type-1-2 relaxation takes place (Fig. 4(c):ii-ii(a)) followed by a spontaneous relaxation (Fig. 4(c):ii(a)-iii and Fig. 4(d): iii) due to DW instability. However, after the 3rd pulse ($j>$3), spontaneous excitation occurs in the absence of E-field (Fig. 4(c): v(a)-vi), rather than spontaneous relaxation. Note that, the origin of spontaneous excitation in this case is not the instability of DW near the grain boundary. In contrast, when two DWs are sufficiently close, then the intermediate domain experiences a non-zero $E_{n}^{int}$ governed by both the DWs. Therefore, total $E_{n}^{int}>1$ at DW interfaces alongside the intermediate domain. Consequently, the intermediate domain becomes unstable and spontaneously switches to $+P$ (Fig. 4(d): v-vi). Such spontaneous excitation continues up to $j$=4, till all the lattices have switched to $+P$. Corresponding $|\Delta^{+}\bar{P}_{n}|_{j}$, $|\Delta^{-}\bar{P}_{n}|_{j}$ and $P_{n}^{acc}$ are shown in Fig. 4(e-h) that present similar trends like case-1. However, an important difference between these two cases is stronger spontaneous excitation and relaxation in case-2 compared to case-1 (see supplementary section for details), which yields relatively abrupt $P$-switching in case-2.

With the understating of $P$-excitation/relaxation processes in an FE grain, we now analyze the $P$-accumulation in HZO by considering an ensemble of grains. The global $E_{C}$ distribution for HZO (80$\mu$m$\times$80$\mu$m) is shown in Fig. 5(a), which we extract from the measured $P$-$E$ curves (discussed in supplementary section). Then, we use each sampled $E_{C}$ as the mean value of a local Gaussian distribution of $E_{C}$ in a grain (like case-2). Considering a large number of grains and multiplying each local $E_{C}$ distribution with the corresponding area fraction, the resultant global distribution of $E_{C}$ is shown in Fig. 5(a).

By considering a sequence of E-field pulses (for $E_{n,max}^{app}$=0.8,1.0,1.2), simulated $\bar{P}_{n}$ and corresponding $P_{n}^{acc}$ are shown in Fig. 5(b-c). While, the signatures of the dynamics of single grain (discussed above) are manifested in HZO (ensemble of grains), two important differences can be observed in HZO: (1) saturation of accumulated $P$ occurs at an intermediate value which increases for higher $E_{n,max}^{app}$, higher $T_{on}$ and or lower $T_{off}$ (Fig. 5(c-e)) and (2) for a long relaxation time, the overall $P$ does not relax completely (Fig.5 (b): gray dashed lines). The former observation is attributed to two processes. First, grains with low mean $E_{C}$ switch completely after sufficient number of pulses and therefore, do not contribute to $P$ accumulation further, leading to intermediate saturation. Second, grains with sufficiently high mean $E_{C}$ exhibit low initial nucleation for a given $E_{n,max}^{app}$ and $T_{on}$. Therefore, given a relaxation time, the grains with higher $E_{C}$ are more likely to relax completely and hence, do not participate in $P$-accumulation. Now, with the increase in $E_{n,max}^{app}$ and $T_{on}$, initial nucleation is enhanced, reducing the probability of complete relaxation in high $E_{C}$ grains . A decrease in $T_{off}$ also has a similar effect on relaxation. This results in the contribution of larger number of grains to $P$-accumulation leading to $P$ saturation at a higher value (Fig. 5(c-e)). For the second observation (incomplete relaxation for large $T_{off}$), the reason is attributed to low mean $E_{C}$ grains that completely switch during the excitation and hence, do not participate in spontaneous relaxation. Note that the large distribution of $E_{C}$ corresponds to the large area of our fabricated HZO sample. However, by scaling the area of HZO, less number of grains along with a compact global distribution of $E_{C}$ can be achieved. Therefore, the $P$-accumulation of a scaled HZO should exhibit less number of saturation levels as well as more abrupt $P$-switching.

In summary, we experimentally demonstrated the accumulative $P$-switching in HZO. Then, developing a phase-field model, we discuss the $P$-switching dynamics by analyzing different stimulated and spontaneous $P$-excitation/relaxation mechanisms governed by domain-domain interaction and DW instability. We attribute the strength and directional change in DW instability as one of the key factors for accumulative $P$-switching. Finally, considering an inter/intra-grain coercive-field distribution in our simulations, we describe the experimentally observed accumulative $P$-switching in HZO and their dependence on E-field pulse attributes.