Thermodynamic equilibrium theory revealing increased hysteresis in ferroelectric field-effect transistors with free charge accumulation

At the core of the theoretical framework of the ferroelectric field-effect transistor (FeFET) is the thermodynamic principle that one can determine the equilibrium behavior of ferroelectric (FERRO) systems using the appropriate thermodynamic potential. In literature, it is often implicitly assumed, without formal justification, that the Gibbs free energy is the appropriate potential and that the impact of free charge accumulation can be neglected. In this Article, we first formally demonstrate that the Grand Potential is the appropriate thermodynamic potential to analyze the equilibrium behavior of perfectly coherent and uniform FERRO-systems. We demonstrate that the Grand Potential only reduces to the Gibbs free energy for perfectly non-conductive FERRO-systems. Consequently, the Grand Potential is always required for free charge-conducting FERRO-systems. We demonstrate that free charge accumulation at the FERRO interface increases the hysteretic device characteristics. Lastly, a theoretical best-case upper limit for the interface defect density DFI is identified. The ferroelectric field-effect transistor, which has attracted much attention for application as both a highly energy-efficient logic device and a non-volatile memory device, has often been studied within the framework of equilibrium thermodynamics. Here, the authors theoretically demonstrate the importance of utilizing the correct thermodynamic potential and investigate the impact of free charge accumulation on the equilibrium performance of ferroelectric-based systems.

The FeFET has also attracted much attention for application as a non-volatile memory device. Owing to its complimentary metal-oxide-semiconductor compatibility in combination with the potential for faster and more power-efficient device operation, the FeFET promises key advantages over conventional nonvolatile memory devices [32][33][34][35][36] . Note that, to focus the discussion and with no loss of generality, the remainder of this Article concentrates on the former application as SS-FeFET.
Salahuddin and Datta's pioneering paper 14 lays the foundation of the stabilized SS-FeFET by describing how a ferroelectric (FERRO) can be stabilized in an intrinsically unstable lowpolarization state by stacking it on top of an insulator (INS) in a metal-FERRO-INS-metal (MFIM) structure. By casting the physics of the MFIM structure in the form of an equivalent electrical circuit of two capacitors in series, it is then argued that FERRO stabilization is achieved if the total capacitance of the system remains positive 14,37 .
A complementary, but more formal, approach to the stabilized SS-FeFET is based on a thermodynamic-equilibrium analysis of ferroelectric systems (FERRO-systems). It revolves around the notion that thermodynamically unstable states of the isolated FERRO can become thermodynamically stable states of the MFIM structure. In the framework of equilibrium thermodynamics, one ascertains the equilibrium states and their stability from the appropriate thermodynamic potential of the system. The importance of analyzing the appropriate thermodynamic potential cannot be understated: a thermodynamic equilibrium analysis based on an inappropriate potential is at risk of producing unphysical results. In the vast majority of works published today, the Gibbs free energy 21,[38][39][40][41][42][43][44] or (though less frequently) the Helmholtz free energy 27,45,46 is used. However, a formal justification has never been provided.
So far, both the equivalent electrical circuit approach and thermodynamic equilibrium approach to the stabilized SS-FeFET assume the absence of free charge accumulation at the FERRO interface. If we envision a stack of a FERRO in series with a regular INS, with the stack covered with a top and bottom metal plate, then the charge on each metal plate is typically different, because of unwanted free charge accumulation at the ferroelectric-insulator (FERRO-INS) interface. One very wellknown source of charge in the semiconductor industry originates from interface traps, mostly related to dangling bonds or other interface defects. Unfortunately, to date, the topic of the impact of interface charge has not been rigorously studied.
In this Article, we present a rigorous theoretical study of the impact of free charge accumulation on the thermodynamic equilibrium theory of stabilized FeFETs, while first providing a formal justification for the appropriate thermodynamic potential. In particular, we demonstrate that for the idealized MFIM structure without any internal accumulation of free charge, under the condition of fixed applied voltage, a Gibbsbased formalism is justified, while a Helmholtz-based approach is not. With the inclusion of free charge accumulation at the FERRO-INS interface because of e.g., charged interface defects, the Gibbs-based formalism breaks down. We show that the Grand Potential is required to correctly describe the thermodynamic equilibrium behavior of MFIM structures or metalferroelectric-semiconductor (MFS) structures with internal free charge accumulation. We demonstrate within the framework of equilibrium thermodynamics that the accumulation of free charge reduces the stability of both ferroelectric MFIM and MFS systems. Lastly, in the case of MFIM systems, we identify a theoretical maximum critical interface trap density above which the system becomes destabilized, giving rise to (increased) hysteretic device operation.

Results and discussion
Coherent uniform ferroelectric systems. In equilibrium thermodynamics 47,48 , a physical system is characterized by a set of intensive state variables, which do not scale with system size (temperature T, pressure P, chemical potential energies μ j and μ 0 α of the j th -mobile and α th -immobile particle, respectively, etc.), and extensive state variables, which scale with system size (volume V, numbers N j and N 0 α of the j th -mobile and α thimmobile particle, respectively, etc.). A physical system is in thermodynamic equilibrium when the state of the system is independent of time and macroscopic "flows" of energy or matter do not exist within the system or at its boundaries with the surrounding. Therefore, the internal production (represented as d i ) of entropy S has vanished (i.e., d i S = 0) and the physical system has evolved to an equilibrium state that minimizes a specific thermodynamic potential Φ that is appropriate for the system under the considered constraints: dΦ = 0. Applying a different set of constraints on the system results to the minimization of a different thermodynamic potential, such as, for example, the internal energy U (¼ R UdV), the Helmholtz free energy F (¼ R F dV), or the Gibbs free energy G (¼ R GdV). For the Reader unfamiliar with the notion of thermodynamic potentials, we provide a short and accessible summary of the key concepts in Supplementary Note 1.
In the FERRO-systems considered in this Article, the FERRO layer is assumed to be a perfectly coherent (=one continuous FERRO material system without grain boundaries or regions of non-FERRO phase interrupting the FERRO phase) and uniform (=without defects or non-uniform external forces) material, characterized by a macroscopic spontaneous polarization density vector P F 49 . The Gibbs free energy density functional (G) of a coherent uniform FERRO is (see Supplementary Note 2): with κ G (∇P F ) defined as: where the summation indices (i and j) run over the x-, y-, and zdirection, ∂ i = ∂/∂x i , ∇P F is a tensor with components ∂ i P F,j , and where the field-free (E F ¼ 0) reference energy functional G 0 has been expanded in terms of P F and ∇P F in the form of the Landau free energy functional 50-54 with Landau parameters α, β, β 0 , κ, κ 0 , and κ 00 , and where E F and ε F are the electric field strength and the background permittivity, respectively. When P F is a scalar, Eq. (2) reduces to κ∇ðP F Þ 2 . Note that the total polarization of the FERRO is P T F ¼ P F þ ðε F À ε 0 ÞE F and that the spatial gradient of P F in the x-direction, which runs perpendicular to the FERRO-INS interface, is assumed to be zero: ∂ x P F = 0. As shown in Supplementary Note 3, the thermodynamic equilibrium behavior of an isolated coherent uniform FERRO at constant electric field E F is determined by the Gibbs free energy, and the characteristic set of P F -E F solutions of a single-domain FERRO ( Supplementary  Fig. S1) is readily found from the first derivative of Eq. (1) over P F : Note that Eq. (3) is temperature-dependent, as illustrated in Supplementary Fig. S2, since α = a 0 (T − T C,0 ), where a 0 is a constant and T C,0 is the Curie temperature.
In the numerical calculations and when deriving the stabilization condition of the FERRO-system, we omit the polarization gradient term ∇P F and we assume that the Landau expansion of the reference free energy functional (such as G 0 in Eq. (1)) remains valid for all values of P F (including around P F ≈ 0). Though often implicitly assumed to be the case in literature, we point out that these are non-trivial assumptions and require uniform single-domain (=without domain wall formation) polarization switching dynamics in a uniform (no stray fields, temperature fluctuations, etc.) FERRO. FERRO-systems, however, can display a mixture of uniform and non-uniform [55][56][57][58][59] (such as domain nucleation and propagation) polarization switching dynamics, including non-coherent switching. As a consequence, the stabilization criteria, obtained in this Article, represent the best-case minimum requirements and are, in principle, only strictly valid for idealized perfectly uniform single-domain FERRO-systems. Qualitatively, however, the results of our Article, demonstrating increased hysteretic behavior with increased free (interface) charge accumulation, are valid beyond these restrictions.
Lastly, throughout this Article, the single-domain ferroelectric is assumed to have the following material parameters at room temperature (300 K): α = −1.1 × 10 9 mF −1 , β = 2.5 × 10 10 m 5 C −2 F −1 , and ε F = 35, which correspond most closely to HfO 2 . However, note that the choice of material parameters is rather arbitrary as there currently does not exist a ferroelectric, including HfO 2 , that can be processed to behave as a perfectly uniform and coherent material.
Intrinsic MFIM system. The first single-domain FERRO-system we consider is an intrinsic (=no defects or free charge carriers) MFIM (i-MFIM) structure with FERRO layer thickness t F = 3.5 nm and INS layer thickness t I = 1 nm and surface area A = L × W ( Supplementary Fig. S3). We ignore variations in the z-direction, such that the physics become two-dimensional, and assume that the voltage (V TB ) applied between the top and bottom electrodes varies sufficiently slow with time, such that the i-MFIM system is studied in the quasi-static approximation. The following conditions then apply:  (4) and (5) enters through the FERRO polarization P F (i.e., we allow for nonzero ∂ y P F in the derivations). Because a charge-and defect-free i-MFIM system is considered here, we have σ FI = 0. From Eqs. (4) and (5), it then directly follows: with C F = ε F /t F and C I = ε I /t I . The i-MFIM structure has been studied in literature using both a Helmholtz or Gibbs free energy approach. However, we demonstrate (see Supplementary Note 4 for a full derivation) that, at constant (non-zero) applied voltage, the total Helmholtz free energy F MFIM is not the appropriate thermodynamic potential of the i-MFIM structure.
In a Gibbs-based approach, the total Gibbs free energy of the system is taken as: We now first provide a formal justification that the total Gibbs free energy G MFIM is the appropriate thermodynamic potential of the i-MFIM structure at constant applied voltage. In Supplementary Note 5, we demonstrate this for the general case of a coherent uniform ferroelectric. For this case, we consider the i-MFIM to consist of N y thin sheets (discretized in the y-direction). The width of the λ th -sheet is Δ λ y , such that: We assume that the FERRO and INS are locally uniform (in the x-direction) within each sheet (i.e., for every state variable u: ∂ x u = 0), and show that the differential form of the Gibbs free energy becomes: Because V TB is constrained in the MFIM structure (dV TB = 0), at which demonstrates that the i-MFIM system evolves to an equilibrium state that minimizes: See Supplementary Note 5 for full details of the derivation.
In the remainder of Supplementary Note 5, we perform a thermodynamic equilibrium analysis of the i-MFIM system in the case of a FERRO layer with uniform spontaneous polarization P F (i.e., ∇P F = 0). The analysis shows that the equation of state (dG MFIM =dP F ¼ 0) reproduces the same equilibrium states (i.e., P F -E F solutions, see Eq. (3)) as for the isolated FERRO layer. Furthermore, it shows that the inherently unstable equilibrium states of the isolated FERRO can only become stable states of the i-MFIM system if the following stabilization condition is satisfied 14,27 : where t F,c is the critical FERRO thickness. Though strictly valid for the i-MFIM system with a coherent uniform FERRO, in practice, Eq. (10) represents a theoretically necessary, but not necessarily sufficient, condition for a hysteresis-free device operation.
Intrinsic MFIM system with fixed charge at FERRO-INS interface. We consider now an equivalent MFIM structure with a fixed (immobile) surface charge density σ FI = σ 0 at the FERRO-INS interface. In Supplementary Note 6, we theoretically demonstrate that this system remains equivalent to the i-MFIM system and that a Gibbs-based approach is still valid, provided one properly accounts for the fixed surface charge density σ 0 when establishing the reference free energy functional of the system (G 0 ).
Free-charge-conducting MFIM system. In practice, the FERRO-INS interface is never completely pristine and, due to the presence of interface traps, some free charge accumulation occurs at the FERRO-INS interface (Fig. 1a). We, therefore, consider a freecharge-conducting MFIM structure (c-MFIM), with leakage through the INS and with identical device dimensions (t F = 3.5 nm and t I = 1 nm), but with surface charge density σ FI . The applied sawtooth voltage V TB between the electrodes has a peak voltage amplitude V P (Fig. 1b).
As per the usual treatment of interface traps 60-62 , the surface charge density is expressed as: ð Þ is the Fermi energy level at position x (Fig. 1c). This expression reflects how a change in the occupation probability of the trap states at the FERRO-INS interface results in a build-up of σ FI . Since macroscopic flows of matter cannot exist in a system in thermodynamic equilibrium, we assume the conductance of the FERRO layer to be negligible such that there is no steady-state leakage current. As there exists no steady-state leakage current, the Fermi level is constant throughout the INS layer ( Fig. 1c), and setting σ FI,0 ≡ qD FI ΔE N , we obtain: For the c-MFIM system the following conditions still apply: (4)) and D I = D F + σ FI (Eq. (5)) with σ FI given by Eq. (11). From this directly follows: Note that the electrostatics of the MFIM structure (Eqs. (12) and (13)) reduce to that of a MFM structure for large D FI . In contrast to the i-MFIM system, at constant applied voltage, the thermodynamic equilibrium behavior of the c-MFIM system is no longer described by the total Gibbs free energy G MFIM (Eq. (9)), as we show in Supplementary Note 7. In Supplementary Note 8 and Supplementary Fig. S5, we demonstrate that an approach based on the Gibbs free energy G MFIM , as expressed by Eq. (9), produces unphysical equilibrium states (Eq. (S75)) and an inaccurate stabilization condition (Eq. (S78)).
Instead of the Gibbs free energy G MFIM , we identify that, at constant applied voltage, the c-MFIM system adopts an equilibrium state that minimizes the Grand Potential Ω MFIM of the system, which is defined as 63 : with the following associated differential form (Ω = ∫ωdV): where the first, second and last summation in Eq. (15) runs over the components of the electric field and displacement vector, all fixed particles (N 0 α ¼ R n 0 α dV), and all mobile particles (N j = ∫n j dV) in the system.
In the c-MFIM structure, the collection of traps at the FERRO-INS interface results in an additional energetic contribution to the thermodynamic potential of the system. In particular, the Grand Potential Ω MFIM is considered to be of the form: Here, Ω FI and ω FI represent the free energy and surface energy density contribution, respectively, to Ω MFIM due to the subsystem at the FERRO-INS interface. As the FERRO-INS interface subsystem only consists of the immobile interface traps and the mobile particles that make up the surface charge σ FI , the differential form of dω FI is of the following form: where s is the entropy surface density, and n 0 α and n j are the surface density of the α th -immobile and j th -mobile particle, respectively.
The differential form of Ω MFIM is shown to reduce to (see Supplementary Note 9 for a detailed derivation): In contrast to the applied voltage, the electric field in the INS (E I ) and the chemical potential at the FERRO-INS interface (φ I t F À Á ) are not explicitly constrained in the c-MFIM system. As a result, it is not self-evident that Eq. (17) becomes zero at equilibrium (d i s λ F ¼ d i s λ I ¼ d i s λ FI ¼ 0 and dV TB = 0). However, because the steady-state leakage current through the FERRO-system has to be negligible (otherwise the system can not be considered to be in a thermodynamic equilibrium state in the first place and has to be studied within the framework of nonequilibrium thermodynamics 47,64 ), the Fermi level in the INS layer is constrained to be constant. As a consequence of the constrained Fermi level, the chemical potential at the FERRO-INS interface (φ t F À Á ) and the electric field in the INS (E I ) become jointly constrained as well and are, therefore, no longer independent from each other (also see Supplementary Fig. S4): This allows us to rewrite Eq. (17) as: As a consequence, we have formally demonstrated that, at fixed applied voltage, the Grand Potential is the appropriate thermodynamic potential and the c-MFIM system evolves to an equilibrium state that minimizes:  (20) is equivalent to Eq. (9) in the limit of D FI = 0, which demonstrates that a Grand Potential-based approach is indistinguishable from a Gibbs-based approach for the i-MFIM system. Next, we utilize Ω MFIM to analyze the equilibrium behavior of the c-MFIM system in the case of a coherent uniform FERRO layer with uniform P F (i.e., ∇P F = 0) and constant interface state density (D FI ) across the entire FERRO-INS interface surface area.
The equilibrium states of the FERRO-system are ascertained from the equation of state (dΩ MFIM =dP F ¼ 0): Using Eq. (12), we retrieve Eq. (3) from Eq. (21), which, once again, demonstrates that the P F -E F characteristics of the uniform single-domain FERRO layer remain unaffected by the adjacent material properties. The stability of the equilibrium states is determined by requiring that d 2 Ω MFIM =dP 2 F > 0: The impact of the interface trap density D FI on the stability of the system and, therefore, the hysteretic device operation is made most clear by using Eq. (22) to establish the equivalent stabilization condition for the c-MFIM system: Note that, for D FI = 0, Eq. (23) reduces to Eq. (10). Figure 3 illustrates the impact of an increasing D FI on the quasi-static simulation results for the "time"-dependent signal of the polarization P F t ð Þ (Fig. 3b), the P F -E F characteristics (Fig. 3c), and the dependence of the electrostatic potential difference over the INS ΔΨ I on the applied voltage V TB , which acts as a proxy for the desired SS-steepening effect in SS-FeFETs (Fig. 3d).
The quasi-static simulations demonstrate that the c-MFIM system reduces to the i-MFIM system in the limit of D FI = 0, and to the MFM system in the limit of D FI → ∞). In addition, the results shown in Fig. 3b-d show that, all else being equal, an increase in D FI results in a (worsened) hysteretic device operation and an internal voltage amplification effect that derives from transient, instead of quasi-static, polarization switching dynamics. Note that our finding of increased hysteretic device operation with D FI is in agreement with theoretical predictions using a kinetics-based model for the polarization switching 65 , and the experimental observation that the memory window of a device increases with the trap density 66 .
In Fig. 4, the critical FERRO thickness t F,c is plotted as a function of D FI and clearly shows how the stabilization condition of the c-MFIM system is negatively impacted by an increasing D FI . The sensitivity of t F to values of α, C I , and ε F is also shown in Fig. 4, demonstrating that, all else being equal, t F,c decreases with increasing values of α, C I , and ε F , which is consistent with earlier reports 27 . Note that, for D FI → ∞, Eq. (23) tends to which shows that the impact of C I becomes negligible for larger values of D FI and reflects the transition from a MFIM to a MFM system. Lastly, note that t F,c also depends on the temperature, as demonstrated in Supplementary Fig. S6.
The origin of the destabilizing effect of a surface charge density σ FI is found in the distortion of the Grand Potential energy landscape (see Eq. (20) and (Fig. 2a, b), which is caused by the dependence of Ω MFIM on D FI . The dependence of Ω MFIM on D FI is traced to the impact of D FI on E F (Eq. (12)) and E I (Eq. (13)). An increasing D FI more strongly affects E I (and therefore Ω I ) than E F (and therefore Ω F ), resulting in a proportionally reduced stabilizing contribution of Ω I to Ω MFIM . As a consequence, a smaller C I (=ε I /t I ) or larger C F (=ε F /t F ), which, in turn, disproportionately affects E F over E I , is required to prevent the  Fig. 3 Surface charge density at the ferroelectric-insulator interface results in increased hysteretic device operation of the free-charge-conducting metal-ferroelectric-insulator-metal (c-MFIM) system. a The Grand Potential energy (Ω MFIM ) landscape of the c-MFIM system (see (Fig. 1a) is plotted along the trajectory (aa) of constant applied voltage V TB = 0 V for different values of the interface state density D FI . The constant-V TB trajectories for D FI = 1 × C I /q 2 and D FI = 5 × C I /q 2 , where C I is the capacitance of the insulator layer and q is the electron charge, are indicated in (Fig. 2. For each case, the stable ferroelectric polarization states are indicated with solid circles. For D FI = 0, the c-MFIM system has a single stable state with ferroelectric polarization P F = 0. However, increasing the interface state density D FI deforms the Grand Potential energy landscape. As a consequence, for larger D FI , the c-MFIM system is characterized by two distinct stable states with non-zero polarization (P 1 F;0 and P 2 F;0 ). b Quasi-static P F -t simulation results show that increasing D FI results in abrupt switching of the ferroelectric polarization P F . (The dashed and solid curves correspond to quasi-static simulations with different initial polarization states (P 1 F;0 or P 2 F;0 ) of the system). c The P F -E F characteristics, where E F is the electric field strength in the ferroelectric layer, and d ΔΨ I -V TB curves, where ΔΨ I is the electrostatic potential difference over the insulator layer, demonstrate the appearance and subsequent worsening of hysteretic device operation with increasing interface state density D FI . The case of D FI = 0 corresponds to the intrinsic MFIM system (Supplementary Fig. S3) and D FI → ∞ to the metal-ferroelectric-metal (MFM) system ( Supplementary Fig. S1).
formation of unstable regimes in the Grand Potential energy landscape.
Lastly, Eq. (23) also allows us to establish the theoretical maximum value of D FI for which a hysteresis-free device operation of a MFIM system with a uniform single-domain ferroelectric layer is possible: For the c-MFIM system considered in this Article (Fig. 1a), D FI,c ≈ 2 × 10 12 eV −1 cm −2 .
Intrinsic MFS system. Next, we consider an intrinsic (i.e., without defect states) metal-ferroelectric-semiconductor (i-MFS) system with semiconductor (SEMI) layer thickness t S (and area A), and with pristine FERRO-SEMI interface (D FS = 0). In contrast to the INS and FERRO layer in the MFIM system, at finite temperatures, it is no longer acceptable to neglect the thermallygenerated electron (n S ) and hole (p S ) free carrier concentrations in the SEMI layer. The local free charge density, ρ S ¼ q p S À n S À Á , affects the local electric displacement field (D S ) in the SEMI through Gauss' law: the local electric field strength and ε S the permittivity of the SEMI layer. Because of the position dependent (x-direction) electric field E S and electric displacement field D S , the SEMI is, inherently, a non-uniform system.
As a consequence, the thermodynamic equilibrium behavior of the i-MFS system is not equivalent to that of the i-MFIM system. In particular, we first demonstrate that, at constant applied voltage V TB , similarly to the c-MFIM system, the Grand Potential (Ω MFS ) is the appropriate thermodynamic potential of the i-MFS system. To allow for a thermodynamic description of the SEMI layer, we use the local equilibrium approximation 47,64 and consider, within each λ th thin sheet, the SEMI to consist of N x elemental volumes dV ι;λ S , discretized in the x-direction with index ι, that are locally uniform (see Supplementary Note 11 and Supplementary  Fig. S7). We ignore variations in the z-direction, such that: Each elemental volume is then characterized by a locally uniform electric field E ι;λ S , electric displacement field D ι;λ S , and free charge density ρ ι;λ S ¼ qðp ι;λ S À n ι;λ S Þ. As a result, Eq. (15) is still valid for each elemental volume dV ι;λ S of the SEMI, such that the differential form of the Grand Potential of the i-MFS system (Ω MFS = Ω F + Ω S ) is found as: With the requirement of a negligible steady-state leakage current through the system, we show that, at equilibrium which demonstrates that, at constant applied voltage (dV TB = 0), the MFS system evolves to an equilibrium state that minimizes the total Grand Potential Ω MFS of the system: where φ S,0 is the zero-field (E S ¼ 0) reference chemical potential and φ S is the local chemical potential of the SEMI layer. Comparing Eq. (27) with the expression for the total Grand Potential of the c-MFIM system (Eq. (20)), it is clear that the thermodynamic equilibrium behavior of the MFS system aligns more closely with the free-charge-conducting MFIM system than the intrinsic MFIM system. We demonstrate in Supplementary Note 12 that, as a first-order approximation, the MFS system is equivalent to the c-MFIM system with mobile surface charge density σ FI . We further show that a) the accumulation of free charge inherently destabilizes the MFS system, and b) that the destabilizing effect becomes more severe with increasing free charge accumulation in the system.
Free-charge-conducting MFIS system. The last single-domain FERRO-system that we discuss is a metal-ferroelectric-insulatorsemiconductor (MFIS) system. Similarly as for the c-MFIM system, we allow for free surface charge accumulation at the FERRO-INS interface (σ FI ) due to a (uniform) interface trap density (D FI ). In addition, besides the intrinsic charge accumulation due to the free carriers in the SEMI, we also allow for free surface charge accumulation (σ IS ) due to a (uniform) interface trap density (D IS ) at the INS-SEMI interface. In the MFIS system, the following conditions apply: free surface charge accumulation due the free (electron and hole) carriers in the SEMI. The surface charge accumulation at the FERRO-INS and INS-SEMI interface is treated similarly as in the case of the c-MFIM system: σ FI = qD FI (E N,FI − E F (t F )) and σ IS = qD IS (E N,IS − E F (t F + t I )). We again assume a negligible conductance of the FERRO layer, such that there exists no steady-state leakage current and the Fermi level is constant throughout the INS and SEMI. We then have for σ FI and σ IS : with σ FI,0 = qD FI ΔE N,FI and σ IS,0 = qD IS ΔE N,IS . We demonstrate that the Grand Potential (Ω MFIS ) is the appropriate thermodynamic potential of the MFIS system at constant applied voltage V TB . Analogous to the MFS system, we use the local equilibrium approximation 47,64 and consider the SEMI layer, within the λ th thin sheet, to consist of N x locally uniform elemental volumes, discretized in the x-direction.
The Grand Potential of the MFIS system is of the form: As the SEMI is in local equilibrium, Eq. (15) for dω S remains valid for each elemental volume. The differential form of Ω MFIS is, therefore, expressed as: At equilibrium (d i s λ , we again require a negligible steady-state leakage current through the FERRO layer and, therefore, the MFIS system. Using Eq. (15) for dω F , dω I , and dω S , and Eq. (16) for dω FI and dω IS , we demonstrate that, at equilibrium, dΩ MFIS reduces to (see Supplementary Note 13 for a detailed derivation): which demonstrates that, at constant applied voltage (dV TB = 0), the MFIS system evolves to an equilibrium state that minimizes the total Grand Potential Ω MFIS , which is expressed as follows (see Supplementary Note 13): Lastly, in Supplementary Note 13, we present a first-order approximation for the thermodynamic equilibrium analysis of the MFIS system based on the Grand Potential Ω MFIS . The equilibrium analysis reveals that, similar to the c-MFIM and MFS system, the MFIS system is destabilized by the free charge accumulation at the FERRO-INS and INS-SEMI interface, and within the SEMI. This becomes clear when examining the (firstorder approximation) stabilization condition of the MFIS system: where C y Σ is defined as: The stabilization condition (Eq. (36)) indicates that the critical FERRO thickness decreases with increasing interface trap density at both the FERRO-INS interface (D FI ) and INS-SEMI interface (D IS ), and with increasing free charge density in the SEMI (n i ).
Note that our formalism remains valid even when the MFIS system is further extended with additional material or contact regions (e.g., the source and drain regions in a FeFET), as long as the steady-state current through the system is negligible.

Conclusion
In conclusion, we have theoretically reconciled the thermodynamic equilibrium theory of perfectly coherent uniform FERRO-systems with free charge accumulation. We have presented the first formal demonstration that the Grand Potential is the appropriate thermodynamic potential to analyze the equilibrium behavior of FERRO-systems, at constant applied voltage.
We have demonstrated that, only in the limiting case of perfectly non-conductive MFIM systems, a thermodynamic equilibrium analysis based on the Grand Potential becomes equivalent to an approach based on the Gibbs free energy. However, we have shown that this is no longer the case in free-charge-conducting FERRO-systems, as the Gibbs free energy approach becomes invalid. As a consequence, the thermodynamic equilibrium behavior of FERRO-systems with free charge accumulation has to be studied using the Grand Potential.
Using the Grand Potential, we have demonstrated that an interface trap density D FI at the FERRO-INS interface destabilizes MFIM systems. Furthermore, we have identified the theoretical upper limit for the critical interface trap density D FI,c in MFIM structures, above which (increased) hysteretic device operation becomes unavoidable. As D FI,c has been derived assuming perfectly uniform single-domain FERRO-systems, D FI,c is expected to further decrease once non-coherent polarization switching dynamics, such as domain nucleation and propagation, are properly taken into account in our analysis.
Lastly, we have demonstrated that the inherent accumulation of free charge carriers in the SEMI layer destabilizes defect-free MFS systems and MFIS systems, which might explain, at least partially, the apparent lack of reported hysteretic-free device operation in MFS systems.
Our findings motivate further research efforts to delineate the theoretical limitations for obtaining a hysteretic-free device operation in advanced FERRO-systems. In particular, future work should focus on extending the presented framework in three steps, as described in Supplementary Note 14. In a first step, the Grand Potential should be extended to allow for a fully spatially non-uniform polarization vector P F . In a second step, the framework should be extended by incorporating a model, inspired by, for example, molecular dynamics or ab-initio calculations, for transient and incoherent polarization switching to describe the time-dependent evolution towards equilibrium. Lastly, an extension towards multi-domain ferroelectrics can be achieved by parallelizing the formalism for independent single-domains.

Methods
Quasi-static simulations. For the FERRO-systems considered in this Article, the quasi-static simulations are performed for the case of an applied voltage V TB with quasi-static sawtooth profile, as shown in Supplementary Fig. S4b, with peak voltage amplitude V P = 0.5 V. We discretize V TB in equidistant steps with increments of ΔV (typically 0.5 mV). For each value of V TB , the quasi-static FERRO polarization P F V TB À Á is determined from the equation of state of the FERROsystem: dΦ/dP F = 0, where Φ represents the appropriate thermodynamic potential of the system. In the case of the i-MFIM system, Φ is the Gibbs free energy per unit area G MFIM (see Eq. (9)), and, in the case of the c-MFIM system, Φ is the Grand Potential Ω MFIM (see Eq. (20)). In a first step, we use MATLAB's fzero algorithm to solve for all existing real-valued roots P F;1 ≤ i ≤ N n o N ≤ 3 to the cubic equation of state (dΦ/dP F = 0). In a second step, we establish the stability of each real-valued root using the corresponding stability condition: d 2 Φ=dP 2 F > 0 (see Eq. (S52) and (22) for the i-MFIM and c-MFIM system, respectively), and only withhold the stable or metastable roots. The quasi-static solution for the FERRO polarization is then determined as: a) in the case one stable root (P F,1 ) remains: and b) in case of both a stable (P F,1 ) and metastable (P F,2 ) root: with P F V TB À ΔV À Á the FERRO polarization of the previous (time) step. At t = 0 s, P F V TB À ΔV À Á corresponds to the assumed initial polarization P F,0 of the FERRO-system (see, for example, Fig. 3a). As a last step, the electric field strength in the FERRO (E F ) and the INS (E I ) are calculated using the corresponding electrostatics equations of the FERRO-system (Eqs. (6) and (7) in case of the i-MFIM system, and Eqs. (12) and (13) in case of the c-MFIM system).

Data availability
The computational scripts that were used to generate the quasi-static simulation results in this work can be obtained from the corresponding author upon reasonable request.