Multi-Domain Negative Capacitance Effects in Metal-Ferroelectric-Insulator-Semiconductor/Metal Stacks: A Phase-field Simulation Based Study

We analyze the ferroelectric domain-wall induced negative capacitance (NC) effect in Metal-FE-Insulator-Metal (MFIM) and Metal-FE-Insulator-Semiconductor (MFIS) stacks through phase-field simulations by self-consistently solving time-dependent Ginzburg Landau equation, Poisson’s equation and semiconductor charge equations. Considering Hf0.5Zr0.5O2 as the ferroelectric material, we study 180° ferroelectric domain formation in MFIM and MFIS stacks and their polarization switching characteristics. Our analysis signifies that the applied voltage-induced polarization switching via soft domain-wall displacement exhibits non-hysteretic characteristics. In addition, the change in domain-wall energy, due to domain-wall displacement, exhibits a long-range interaction and thus, leads to a non-homogeneous effective local negative permittivity in the ferroelectric. Such effects yield an average negative effective permittivity that further provides an enhanced charge response in the MFIM stack, compared to Metal-Insulator-Metal. Furthermore, we show that the domain-wall induced negative effective permittivity is not an intrinsic property of the ferroelectric material and therefore, is dependent on its thickness, the gradient energy coefficient and the in-plane permittivity of the underlying insulator. Similar to the MFIM stack, MFIS stack also exhibits an enhanced charge/capacitance response compared to Metal-Oxide-Semiconductor (MOS) capacitor. Simultaneously, the multi-domain state of the ferroelectric induces non-homogeneous potential in the underlying insulator and semiconductor layer. At low applied voltages, such non-homogeneity leads to the co-existence of electrons and holes in an undoped semiconductor. In addition, we show that with the ferroelectric layer being in the 180° multi-domain state, the minimum potential at the ferroelectric-dielectric interface and hence, the minimum surface potential in the semiconductor, does not exceed the applied voltage (in-spite of the local differential amplification and charge enhancement).

www.nature.com/scientificreports www.nature.com/scientificreports/ According to Landau-Khalatnikov equation 1 , FE polarization (P) versus electric-field (E) characteristics exhibit an unstable negative slope (negative dP dE / ) region as shown in Fig. 1(c). According to ref. 1 , such an unstable region in FE can be stabilized in a heterogeneous system (i.e. FE-dielectric (DE) stack) so that a homogeneously suppressed polarization (P = 0) can be obtained by suppressing the depolarization energy. However, under certain conditions, it may be more natural to form multiple domains with positive and negative P separated by domain-walls (DWs) to suppress the depolarization energy of the system 2 . Recently, DW motion-based P -switching in multi-domain (MD) FE has been identified as a possible mechanism for obtaining static effective negative capacitance in FE [2][3][4][5] . Such DW-induced NC effect has been theoretically predicted in refs. [6][7][8] showing that the soft-DW displacement can lead to an effective negative permittivity of FE in the presence of the interfacial dead layer. Further, a similar effect has been discussed and analyzed through phase-field simulations predicting a hysteresis-free NC path in FE by considering a moving DW in a metal-FE-metal capacitor 9 and DE-FE-DE superlattice 10 . Additionally, an analytical model for DW-induced NC has been proposed for DE-FE-DE superlattice in ref. 10 suggesting that the effective NC path is dependent on the DE thickness (T de ), which contrasts with the analytical model presented in ref. 11 . Our phase-field simulations show that the soft-DW motion-based effective NC path in FE is independent of T de , but it depends on the in-plane permittivity of the DE layer, which is in agreement with ref. 11 . To further identify such inter-dependencies of FE NC behavior on the properties of the constituent FE and DE layers in such heterostructures, we extensively analyze DW-induced NC effect in MFIM ( Fig. 1(a)) based on phase-field simulations (beyond what has been explored so far) and establish its dependence on FE thickness, gradient energy coefficient, and DE permittivity and thickness. Furthermore, we, for the first time, develop a self-consistent 2D phase-field simulation framework for Metal-Ferroelectric-Insulator-Semiconductor (MFIS) stack ( Fig. 1(b)). Utilizing our framework, we investigate DW induced NC effect in the MFIS stack, its effect on the semiconductor potential and its dependency on key material/device parameters.
In our phase-field simulation, we solve the 2D time (t)-dependent Ginzburg-Landau (TDGL) equation 10 , Poisson's equation and semiconductor charge equations, yielding self-consistent solutions for polarization (P x z t ( , , )), potential (Φ x z t ( , , )) and charge (ρ x z t ( , , )), where z and x are the axis parallel to the thickness and length of the stack, respectively. Even though the simulation framework can capture the transient dynamics of polarization switching, in this work, we focus on analyzing the steady-state solution of P x z ( , ) for an applied voltage. For the FE material, we consider Hf 0.5 Zr 0.5 O 2 (HZO) and the corresponding Landau's free energy coefficients (α, β and γ) are extracted from measured P-E characteristics 12 shown in Fig. 1(c). For the gradient energy coefficient (g ), a range of values are considered (see Methods section) as the actual value is still unknown for www.nature.com/scientificreports www.nature.com/scientificreports/ HZO. We assume that the spontaneous P direction in FE is along the thickness of the film (z-axis), which is parallel to the c-axis of the orthorhombic crystal phase of HZO 13,14 . For DE, we consider SiO 2 , Al 2 O 3 and HfO 2 , and for the semiconductor, we consider undoped silicon. The details of the simulation methodologies, parameters and calibration are discussed in the 'Methods' section. In addition, the symbolic representations and corresponding equations of different parameters are summarized in a table shown in Fig. 1(d). Utilizing this framework, we analyze the multi-domain formation and the applied voltage-dependent P-switching in MFIM and MFIS stacks (shown in Fig. 1(a,b)), which we discuss in the subsequent sections.

Multi-domain formation in MFIM stack
Let us start by considering an MFIM stack with an applied voltage (V app ) of 0 V. It is well known that in MFIM stack, spontaneous polarization (P) appears at the FE-DE interface, which leads to a voltage drop across the DE layer 10 . As a result, an E-field appears in FE opposite to the P direction, called depolarization field, which leads to an increase in the depolarization energy density. Let us define the E-field along the thickness of the FE film (z-axis) as E z fe , and the depolarization energy density as . However, f dep can be suppressed with the formation of periodic 180° domains of alternating P-directions (P↓ and P↑). Simulation result considering such a scenario is shown in Fig. 2 Here, ↑ (↓) sign denotes the +z (−z) direction. In this multi-domain (MD) state, the magnitude of the local E z fe , (at a particular point in the FE) is greatly reduced due to stray fields (in-plane E-field, E x fe , ) between P↑ and P↓ domains, as shown in Fig. 2

(b). While this decrease in local E z fe
, is larger near the domain walls (DWs) compared to inside of the domains, the suppression of average E z fe , is significant across the entire length of the stack (along the x-direction). The resultant decrease in f dep , however, comes at the cost of DW energy density (F dw ), which is comprised of (a) the electrostatic energy density ( f E ). Note that the magnitude of P inside of a domain also varies along the z-axis exhibiting a minima at the DE interface and gradually increases in the bulk FE (away from the DE interface). This induces a bound charge density dP dz / b ρ = − and further suppresses the E z fe , (and hence, f dep ) inside of the domain. However, this additional suppression of E z fe , occurs at the cost of an increase in the z-component of gradient energy density ( = × f g dP dz ( / ) z grad , 2 ). Our simulations show that f z grad , occurs in FE both in the MD (co-existing P↑ and P↓) and poled (either P↑ or P↓) states. In the MD state, achieved by suppressing f dep at the cost of f dw and f z grad , while minimizing the overall system energy, the intricate interactions of these energy components with each other as well as the free energy ( f free ) play a key role in determining the www.nature.com/scientificreports www.nature.com/scientificreports/ charge response of FE in the MFIM stack and its dependence on the device/material parameters, as discussed subsequently. Now, let us first consider the implication of FE thickness, T fe on the formation of MD state. The P configurations of FE in MFIM stacks shown in Fig. 2(c), suggest that the number of domains and DWs (within a certain length) increases with the decrease in T fe . As T fe decreases, f z grad , increases as a similar P variation along z-axis (i.e. similar P maxima in the bulk and minima in the interface) occurs within a lower T fe . One of the possible ways to reduce f z grad , could be decreasing P variation by increasing P magnitude in the interface, but this would increase f dep . As an alternative, an increase in the number of DWs leads to reduced domain width and thus, reduces the P magnitude in the bulk region. Therefore, an increase in f z grad , due to T fe scaling is mitigated by increasing the number of DWs. In this case, suppression of f dep becomes more significant inside of a domain (as P decreases in magnitude) and also on an average (as the number of DWs increases). At the same time, with decreasing T fe , as the number of DWs increases, the nature of DW changes from hard to soft type (Fig. 2(c)). Here, the term hard-DW implies that the spatial variation in P is localized and abrupt (dP dx / is high) only near the DWs. Thus, f x grad , is non-zero in the DW and zero (or very small) within the domains. In contrast, in a soft-DW, the P variation is more gradual (dP dx / is low) and the effects of DW diffuses along the length-scale of a domain. That implies f x grad , becomes non-zero inside of a domain. Similar to the effect of T fe , the gradient energy coefficient (g ) also determines the number of DWs in FE. As f x grad , is one of the components of f dw , a decrease in g leads to lower DW energy cost and, thus, leads to larger number of domains and DWs. Such an increase in the number of DWs density with the decrease in g is shown in Fig. 2(d). Further, the nature of DW changes from hard to soft type as g increases. This is because, for higher g , dP dx / decreases to compensate for the f x grad , and thus the P-variation becomes more gradual and diffuses within the domain. The nature of DW (i.e. soft or hard) in MFIM for different g and T fe is illustrated in Fig. 2(e) showing that the critical T fe for hard to soft DW transition decreases with a decrease in g. Further, if T fe is scaled below a critical value, a single domain (SD) state with homogenous P = 0 stabilizes ( Fig. 2(c): T fe = 1 nm), where the suppression of f dep occurs at the cost of f free rather than f dw . For suppressing f dep , if f dw is higher than f free then the SD state is preferred over the MD state. In addition, the critical T fe , for MD to SD transition, decreases with a decrease in g as shown in Fig. 2(e). As f dw decreases with g , a lower T fe is needed to go beyond f free . Note that if g is very small, the critical T fe can potentially become so small that the SD state or soft MD state may not be physically realizable. For example, if g m V C 0 1 10 / 9 3 < . × − , the critical T fe for SD state (0.25 nm) and soft MD state (0.5 nm) is lower than or comparable to a unit cell height of HZO.
In the above discussion, the stability of SD state with P = 0 may or may not imply that the ferroelectricity of the FE layer will be retained and the same Landau coefficients can be used to calculate the polarization switching characteristics. Rather, one needs to calculate the f free of the ferroelectric phase at P = 0 along with the paraelectric and anti-ferroelectric phases (i.e. orthorhombic, monoclinic and tetragonal phase for HZO). In such a case, the state with lowest f free will be stabilized and if the stable phase is either paraelectric or anit-ferroelectric, then the Landau coefficients for ferroelectric phase can not be used to simulate the polarization switching characteristics. However, such calculation requires first-principle simulation, which is out of the scope of our current study. Hence, we limit our analysis to the MD states.
Further, in the MD state, the nature of DW plays an important role in E-field driven DW motion. To displace the hard-DW, the applied E-field needs to be higher than a critical value called coercive field of DW motion, E c dw , . Therefore, the hard-DW motion based P-switching is hysteretic, because, the applied E-field needs to be higher than the positive E c dw , for forward DW motion and lower than the negative E c dw , for reverse DW motion 15 . In contrast, | | E c dw , is infinitesimally small (~0) for soft-DW 15 and hence, non-hysteretic DW motion is possible. As in this work, our focus is on analyzing the non-hysteretic NC effect, we restrict our discussion only for soft-DW motion based P-switching.

MD-NC Effect in MFIM Stack
Let us begin by discussing P-switching in the MFIM stack with soft-DW. The simulated average charge density (Q) vs applied voltage (V app ) characteristics is shown in Fig. 3 Here, Q is calculated as the average displacement field at the Metal-DE interface based on following equation.
where, E z de , is the z-component of E-field at Metal-DE interface and l is the length of the stack. For | | < V 2 V app , a continuous Q-V app path exists when the FE is in the MD state and P-switching takes place through DW motion (see Fig. 3(b)). If | | V app is increased above 2 V, MD state (P) switches to the poled state (either P↑ or P↓). Now, with decreasing V app | |, MD state forms from the poled state at a lower | | V app (0.9 V) and that induces a hysteresis in the Q-V app characteristics (discussed in the Supplementary Section). Therefore, for non-hysteretic operation, the MD state needs to be retained by limiting V app . Interestingly, in the MD state, Q is higher in the MFIM stack compared to the MIM (Metal-Insulator-Metal) at the same V app as shown in Fig. 3(a). This implies that the effective capacitance of the MFIM stack is higher than MIM. In a static scenario, such a phenomena is only possible if the FE www.nature.com/scientificreports www.nature.com/scientificreports/ layer acts as an effective negative capacitor. The Q V fe avg − characteristics are shown in Fig. 3(a). In the MD state, as the potential drop across the FE layer is non-homogeneous along the x-axis ( Fig. 2(b)), the average potential drop across the FE layer (V fe avg ) is calculated using the following equation.
is the FE-DE interface potential Q is calculated by taking the average displacement field at the FE-DE interface which provides exactly same Q as Eq. 1 for the same V app . Figure 3( Fig. 3(a) www.nature.com/scientificreports www.nature.com/scientificreports/ P↓ domains exhibit E↑ and P↑ domains exhibit E↓). Note that f x grad , is non-zero inside of the domain (due to DW diffusion in soft-DW) and that causes the P to decrease in magnitude (discussed earlier). Now, with the increase in V app , P↓ domains grow and P↑ domains shrink in size, due to positive stiffness of DW motion 8 . As the DW moves away from P↓ domain and towards the P↑ domain, f x grad , in P↓ domain decreases and in P↑ domain increases. Due to this effect as well as increase in V app , the magnitude of local P in P↓ domain increases and in P↑ domain decreases. Our simulation shows that as a result of this, the depolarizing field (E z fe , ) in P↓ domain increases and in P↑ domain decreases in magnitude. This implies f dep increases (decreases) in P↓ (P↑) domain. The increase (decreases) in f dep in P↓ (P ↑ ) domains is possible as it is accompanied by a decrease (increase) in f x grad , .
As the oppositely directed local E-field in FE increases (decreases) with the increase (decrease) in local P in both P↓ and P↑ domains, the effective local permittivity of the domains (ε z fe eff , ) become negative. At the same time, in the DW, the asymmetry in P distribution (due to unequal P↑ and P↓ domain sizes and P magnitudes) causes F dw (comprised of f x grad , and f x elec , ) to decrease compared to the symmetric P distribution 2 (at V app = 0 V). Such a decrease in F dw allows a further increase in average E z fe avg , (an increase in depolarization energy) in the DW, while the average-P (directed opposite to E z fe avg , ) in the DW increases (due to unequal P magnitudes in P↑ and P↓ domain). As a consequence, the permittivity of the DW region also becomes negative. These effective local (and non-homogeneous) negative permittivity ( z fe eff , ε ) of the domain and DW regions give rise to an average effective negative permittivity in the FE layer, i.e. ε < − 0 z fe eff avg , . Here, it is important to note that the appearance of negative effective permittivity is essentially an apparent phenomena of change in long range interaction of P, its gradient and/or DW energy under DW motion. In particular, the change in P in MFIM is not directly driven by the local E-field, rather, the change in P is driven by the applied E-field induced domain-wall motion. Therefore, the change in local E-field is the effect of change in P and not the opposite. In other words, the depolarizing E-field appears depending on the change in P induced by DW motion. Even though, such phenomena leads to a negative effective permittivity of the FE layer, the susceptibility of the FE layer and the whole system (MFIM stack) is positive with respect to the applied E-field. That implies, the change in polarization is always in the direction of the change in applied E-field.
As we have identified that the F dw plays a crucial role in providing negative ε − z fe eff avg , , therefore, it is intuitive that the NC effect is dependent on its components, i.e. . As the f x grad , increases with the increase in g, a higher energy modulation is achieved by displacing the DW, which further provides a higher increase (or decrease) in f dep in P↓ (or P↑) domains, leading to larger NC effect. Similarly, dP dx / increases as the number of domains and the DWs increase with the decrease in T fe (discussed before). Therefore, f x grad , increases and provides an increased NC effect with decreasing T fe (Fig. 3(d)). However, the soft-DW induced NC path does not depend of T de (Fig. 3(e)). This is because, in the MD state, the average depolarization field (which is zero at V app = 0) as well as f x grad , and f x elec , are independent of T de within the limit of soft-DW. Interestingly, the MD-NC path does depend on the relative DE permittivity ( de ε ) as shown in Fig. 3(f). This is because the in-plane E-field, E x fe , in the DW needs to satisfy the in-plane boundary condition at the FE-DE interface, which is , and E x fe , are the in-plane E-field in DE and FE, respectively. As the E x de , increases with the decrease in ε de (considering similar P difference between two consecutive domains), therefore, E x fe , also increases in FE, which further increases the f x elec , stored in the DW. Therefore, the F dw increases and hence, NC effect increases with the decrease in ε de as shown in Fig. 3(f). From this analysis, we can summarize that, (i) an FE material with higher g, (ii) T fe scaling and/or (iii) using DE materials with low ε de are key design knobs to enhance DW-induced NC effect in MFIM stack. Note that in all of the cases discussed above, the MD NC path does not coincide with the Landau path (Fig. 3(c-f)) and the MD NC effect is less ( is less) compared to the NC effect that corresponds to Landau path.
As the MD NC path is dependent on T fe , g and ε de , therefore, the charge enhancement characteristics also depend on them as shown in Fig. 4(a-c). For simplicity, we only show the charge response for V 0 app > . Now, a relation between the charge response in MFIM (Q MFIM ) and MIM (Q MIM ) can be written as the following equation. Here C MIM (=ε ε T / de de 0 www.nature.com/scientificreports www.nature.com/scientificreports/ ment shows almost similar characteristics ( Fig. 4(b) for T fe = 10 nm and 15 nm). Now, with the increase in ε de (but same C MIM , obtained by increasing T de proportionately, so the Q-V app characteristics of MIM remains same), charge enhancement increases as shown in Fig. 4(c). This is because, 1 increases with the decrease in ε de . Based on the above discussion, we emphasize that the negative effective permittivity of the MD state is not an intrinsic material parameter of FE, rather, it depends on the physical thickness of the FE film, its gradient energy coefficient and the permittivity of the underlying DE material. It is important to note that, this conclusion is different than a single domain model 1 , where the negative permittivity originates from negative slope of the Landau-Khalatnikov (LK) equation and hence, remains constant irrespective of T fe and ε de . Further, ε − z fe eff avg , is a non-local quantity and can only describe the average characteristics. Moreover, the value of z fe eff avg , ε − obtained for a particular T fe and ε de can not be used to calculate the average charge response of MFIM/MFIS stack with any other T fe and de ε due to the dependency of z fe eff avg , ε − on these parameters. However, z fe eff avg , ε − can be used to calculate the average charge response of MFIM/MFIS stack with different T de within the limit of soft DW formation. Furthermore, the local effective permittivity, z fe eff , ε of is not an intrinsic material property and hence, spatially varies within the FE layer 3,4 . As the DW moves with the applied voltage, the local effective permittivity value changes with V app . As, z fe eff , ε is not a spatially static quantity, therefore, one cannot directly use this in a capacitor equation to analyze the local charge response of the heterostructures (i.e. MFIM and MFIS stack). However, z fe eff avg , ε − can be used in a capacitor equation to calculate the average charge response of MFIM/MFIS stack (within the limit of MD state) as in Eqs. 3-4. Therefore, the reason for introducing ε z fe eff , and z fe eff avg , ε − in this work is to analyze the implication of different parameters and to emphasize that, one should not use this characteristics as an intrinsic property of FE.

MD-NC Effect in MFIS Stack
So far, we discussed how the DW-induced effective NC in FE can enhance the overall charge response of MFIM stack. Next, we turn our attention to the MFIS stack with an undoped silicon as the semiconductor layer. To compare the MFIS results with conventional MOS capacitor, we also simulate MIIS (Metal-HfO 2 -SiO 2 -Si) and MIS (Metal-SiO 2 -Si) stacks. The silicon layer thickness of 10 nm is considered for all the simulations of MFIS, MIIS, and MIS stacks. The simulated Q-V app and C-V app (capacitance, = C dQ dV / app ) responses are shown in Fig. 5(a,b), which illustrate an enhanced charge and capacitance response of MFIS compared to the MIIS and MIS stacks. We attribute this enhanced charge/capacitance response to the negative z fe eff avg , ε − of FE that we discussed earlier for MFIM. Now, to analyze the effects of T fe , Q-V app characteristics for different T fe are illustrated in Fig. 5(d) showing minor enhancement in the charge response with an increase in T fe . To understand this, a relation can be derived between the charge response in MFIS (Q MFIS ) and in MIS (Q MIS ) when z fe eff avg , ε − is negative as follows. The Q-V app characteristics (Fig. 5(c)) show that the MFIS charge response enhances with the increase in g and are www.nature.com/scientificreports www.nature.com/scientificreports/   increases) with the increase in g, as we discussed earlier in the context of MFIM.
The overall enhancement in charge/capacitance response of MFIS stack (compared to MIS and MIIS) can be easily understood from the effective negative ε − z fe eff avg , of FE. However, for FE-FET operation, it is also important to analyze the semiconductor surface potential (Ψ) in MFIS, especially, when FE is in MD state ( Fig. 6(a)). In fact, Ψ(x) in MFIS is non-homogeneous as shown in Fig. 6(b) at V app = 0 V. To understand this, let us consider the potential at the FE- is directed opposite to the local P and exhibits a non-homogeneous profile along the x-direction due to periodic P↑ and P↓ domains. Therefore, V x ( ) int becomes non-homogenous and exhibits a maxima (max-V int ) and minima (min-V int ) corresponding to the center of P↓ and P↑ domains, respectively. This non-homogeneity in V x ( ) int induces a spatially varying Ψ(x) as shown in Fig. 6(b), which, in turn exhibits a maxima (max-Ψ) and minima (min-Ψ). This further leads to local accumulation and co-existence of electrons or holes in the undoped Si layer (Fig. 6(c)). Note, such a spatially varying charge profile has been experimentally shown in ref. 16 for FE-semiconductor interface, when FE is in MD state. Now, with the increase in V app (~1.2 V), P↓ domains grow and P↑ domains shrink in size leading to an overall increase in average P (Fig. 6(d)). Simultaneously, min/max-V int increases ( Fig. 7(a)) and at the same time exhibits a differential amplification ( > dV dV / 1 int a pp ) as shown in Fig. 7(b). Here the local differential amplification in min/max-V int can be attributed to the effective local negative permittivity of FE in the P↓ and P↑ domains (discussed for MFIM). Now, as V int increases, Ψ everywhere at the Si interface increases and becomes positive (but still remains non-homogeneous, see Fig. 6(e)). Therefore, electron density (n) dominates over hole density (p) locally and globally ( Fig. 6(f)). Note that the increase in n causes the non-homogeneity in Ψ to decrease (Fig. 6(e,f)) compared to V app = 0 V ( Fig. 6(b,c)). The Ψ for MFIS, MIIS and MIS stacks for V app = 0 V and 1.2 V are shown in Fig. 6(b,e). At V V 0 app = , the max(min)-Ψ in MFIS is higher(lower) than the MIIS and MIS stacks. At V app = 1.2 V, the max-Ψ in the MFIS is higher than the MIIS and MIS and the min-Ψ in MFIS is higher than MIIS but lower than the MIS. This can be understood from the following discussion. As in the MD state, E z fe , the min(max)-V int is always lower(higher) than V app in MFIS (see Fig. 7(a)). Note that, this statement is also true for MFIM and is discussed in the Supplementary Section. Now, in the MIS stack, DE potential is directly driven by V app and hence V V int a pp = . Therefore, min-V int of MFIS is always less than V int (=V app ) of MIS as shown in Fig. 7(a). In addition, dΨ/dV int is <1 and equal for both MFIS and MIS due to the same positive capacitance of the DE layer. As a consequence, the min-Ψ of MFIS is inevitably lower than the Ψ of MIS, when the FE is in 180° MD state (Fig. 7(c)). However, in MIIS, the V int (HfO 2 -SiO 2 interface potential) is not directly driven by V app and due to the positive capacitance of the HfO 2 layer, dV dV / 1 int a pp < and < V V int a pp ( Fig. 7(b,c)). Now, considering the differential amplification of min- ) as shown in Fig. 7(b), the min-V int of MFIS becomes higher than the V int of MIIS beyond a certain V app (Fig. 7(c)). As a result, min-Ψ of MFIS becomes higher than the Ψ of MIIS at > V 1 V app ) as shown in Fig. 7(c). Briefly, in MFIS, www.nature.com/scientificreports www.nature.com/scientificreports/ the min-Ψ can exceed the Ψ in MIIS but remains lower than the MIS, while the max-Ψ is always higher than the Ψ in MIIS and MIS. Now, let us make a rough assumption that the channel current in FEFET will be mostly dependent on the min-Ψ as that determines the highest potential barrier seen by the source electrons. Then, based on the above discussion, we can expect that the OFF current (at V app = 0 V) of FEFET will be significantly less compared to the MIS/MIIS-FET, and the ON current (V app ~ 1.2 V) will be higher than the MIIS-FET but comparable to MIS-FET. As the Ψ is highly non-homogeneous in MFIS stack in the low voltage regime, calculation of SS of FEFETs needs further exploration by considering source/drain regions along with the DW-induced non-homogeneous semiconductor potential and solving the transport equations to obtain the impact of the MD FE on the FEFET characteristics.

Conclusion
In summary, by performing phase-field simulation, we show that the energy stored in FE DW can be harnessed to enhance the capacitance of the MFIM and MFIS stack, where the soft-DW displacement leads to a static and hysteresis-free negative capacitance in the MD FE. Our analysis indicates that the effective negative permittivity of the FE layer is dependent on the FE thickness, gradient energy coefficient, in-plane permittivity of the DE and is independent of DE thickness within the limit of soft-DW. Further, the DW-induced NC effect can lead to an enhanced charge/capacitance response in MFIS stack compared to MIS/MIIS stack. However, such a charge/ capacitance enhancement in MFIS does not guarantee an enhanced local Ψ in Si compared to MIS. In fact, Ψ becomes spatially varying due to the MD nature of FE and the variation is higher at low applied voltages. In addition, we discuss that the minimum Ψ in MFIS can exceed the Ψ in MIIS but remains smaller than the MIS. Nevertheless, considering the local differential amplification of V int (i.e. − > d min V dV ( ) / 1 int a pp ), the on/off current ratio of FEFET can potentially exceed the MIS/MIIS-FET. Since the non-homogeneity in Ψ is absent in conventional MOS capacitor (and MOSFET), therefore, as future work, it will be important to investigate the impact of such potential profile on the low voltage conduction of FEFETs. ρ . We solve these equations (Eqs. 7-9) in real space by employing finite difference method with grid spacing of Δx = Δz = 0.5 nm by considering Dirichlet boundary condition at the top and bottom metal contacts (i.e. by defining the voltage) and periodic boundary condition at the left and right edges of the system. For all the simulation results presented in this work, the lateral (along x-axis) dimension of the system is considered as 100 nm.