Introduction

The miniaturization of the metal-oxide-semiconductor field-effect transistor (FET) has pushed transistor technology to the point where thermally excited carriers prevent any further reduction of the supply voltage due to the related leakage current increase1. Multiple strategies are pursued to go beyond the limit of a thermally-restricted subthreshold slope. A first strategy explores device concepts based on novel charge injection mechanisms, such as tunnel FETs2,3,4,5,6,7,8 and impact ionization FETs9,10. A second strategy exploits a positive feedback mechanism that amplifies the electrostatic control of the gate over the channel surface potential, such as in the nanoelectromechanical FET11,12,13 and the steep-slope ferroelectric FET (SS-FeFET)14,15,16,17. In a SS-FeFET, the internal voltage amplification is realized either by means of abrupt (non-stabilized SS-FeFET)18,19,20,21,22,23, or gradual polarization switching (stabilized SS-FeFET)24,25,26,27,28,29,30,31 of the ferroelectric.

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 non-volatile memory devices32,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 paper14 lays the foundation of the stabilized SS-FeFET by describing how a ferroelectric (FERRO) can be stabilized in an intrinsically unstable low-polarization 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 positive14,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 energy21,38,39,40,41,42,43,44 or (though less frequently) the Helmholtz free energy27,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 well-known 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 Gibbs-based 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 metal-ferroelectric-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 thermodynamics47,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 \({\mu }_{\alpha }^{\prime}\) of the jth-mobile and αth-immobile particle, respectively, etc.), and extensive state variables, which scale with system size (volume V, numbers Nj and \({N}_{{\rm{\alpha }}}^{\prime}\) of the jth-mobile and αth-immobile 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 di) of entropy S has vanished (i.e., diS = 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 (\(=\ \int {\mathcal{U}}{\rm{d}}V\)), the Helmholtz free energy F (\(=\ \int {\mathcal{F}}{\rm{d}}V\)), or the Gibbs free energy G (\(=\ \int {\mathcal{G}}{\rm{d}}V\)). 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 PF49. The Gibbs free energy density functional (\({\mathcal{G}}\)) of a coherent uniform FERRO is (see Supplementary Note 2):

$${{\mathcal{G}}}_{{\rm{F}}}\left({{\bf{P}}}_{{\rm{F}}},{{\mathcal{E}}}_{{\rm{F}}}\right)= \, \alpha {\left|{{\bf{P}}}_{{\rm{F}}}\right|}^{2}\,+\,\beta {\left|{{\bf{P}}}_{{\rm{F}}}\right|}^{4}\,+\,{\beta }^{\prime}\mathop{\sum}_{i\,\ne\, j}{P}_{{\rm{F}},i}^{2}{P}_{{\rm{F}},j}^{2}\\ +{\kappa }_{{\rm{G}}}\left(\nabla {{\bf{P}}}_{{\rm{F}}}\right)\,-\,\mathop{\sum} _{i}{P}_{{\rm{F}},i}{{\mathcal{E}}}_{{\rm{F}},i}\,-\,\frac{{\varepsilon }_{{\rm{F}}}}{2}{\left|{{\bf{{\mathcal{E}}}}}_{{\rm{F}}}\right|}^{2},$$
( 1)

with κG(∇PF) defined as:

$${\kappa }_{{\rm{G}}}\left(\nabla {{\bf{P}}}_{{\rm{F}}}\right)= \, \kappa \mathop{\sum} _{i}{\left({\partial }_{i}{P}_{{\rm{F}},i}\right)}^{2}\,+\,{\kappa }^{\prime}\mathop{\sum} _{i\,\ne\, j}\left({\partial }_{i}{P}_{{\rm{F}},i}\right)\left({\partial }_{j}{P}_{{\rm{F}},j}\right)\\ +{\kappa }{^{\prime\prime} }\mathop{\sum} _{i\ <\ j}{\left({\partial }_{i}{P}_{{\rm{F}},j}\,+\,{\partial }_{j}{P}_{{\rm{F}},i}\right)}^{2},$$
( 2)

where the summation indices (i and j) run over the x-, y-, and z-direction, ∂i = ∂/∂xi, ∇PF is a tensor with components ∂iPF,j, and where the field-free (\({{\bf{{\mathcal{E}}}}}_{{\rm{F}}}\,=\,0\)) reference energy functional \({{\mathcal{G}}}_{{\rm{0}}}\) has been expanded in terms of PF and ∇PF in the form of the Landau free energy functional50,51,52,53,54 with Landau parameters α, β, \({\beta }^{\prime}\), κ, \({\kappa }^{\prime}\), and \({\kappa }^{\prime\prime}\), and where \({{\mathcal{E}}}_{{\rm{F}}}\) and εF are the electric field strength and the background permittivity, respectively. When PF is a scalar, Eq. (2) reduces to \(\kappa \nabla {({P}_{{\rm{F}}})}^{2}\). Note that the total polarization of the FERRO is \({{\bf{P}}}_{{\rm{F}}}^{{\rm{T}}}\,=\,{{\bf{P}}}_{{\rm{F}}}\,+\,({\varepsilon }_{{\rm{F}}}\,-\,{\varepsilon }_{0}){{\bf{{\mathcal{E}}}}}_{{\rm{F}}}\) and that the spatial gradient of PF in the x-direction, which runs perpendicular to the FERRO-INS interface, is assumed to be zero: ∂xPF = 0. As shown in Supplementary Note 3, the thermodynamic equilibrium behavior of an isolated coherent uniform FERRO at constant electric field \({{\mathcal{E}}}_{{\rm{F}}}\) is determined by the Gibbs free energy, and the characteristic set of PF-\({{\mathcal{E}}}_{{\rm{F}}}\) solutions of a single-domain FERRO (Supplementary Fig. S1) is readily found from the first derivative of Eq. (1) over PF:

$$2\alpha {P}_{{\rm{F}}}\,+\,4\beta {P}_{{\rm{F}}}^{3}\,=\,{{\mathcal{E}}}_{{\rm{F}}}$$
( 3)

Note that Eq. (3) is temperature-dependent, as illustrated in Supplementary Fig. S2, since α = a0(T − TC,0), where a0 is a constant and TC,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 ∇PF and we assume that the Landau expansion of the reference free energy functional (such as \({{\mathcal{G}}}_{{\rm{0}}}\) in Eq. (1)) remains valid for all values of PF (including around PF ≈ 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-uniform55,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 × 109 mF−1, β = 2.5 × 1010 m5C−2F−1, and εF = 35, which correspond most closely to HfO2. However, note that the choice of material parameters is rather arbitrary as there currently does not exist a ferroelectric, including HfO2, 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 tF = 3.5 nm and INS layer thickness tI = 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 (VTB) 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:

$${V}_{{\rm{TB}}}\,=\,{{\mathcal{E}}}_{{\rm{F}}}\left(y\right){t}_{{\rm{F}}}\,+\,{{\mathcal{E}}}_{{\rm{I}}}\left(y\right){t}_{{\rm{I}}}$$
( 4)
$${D}_{{\rm{I}}}\left(y\right)\,=\,{D}_{{\rm{F}}}\left(y\right)\,+\,{\sigma }_{{\rm{FI}}}\left(y\right),$$
( 5)

with Eq. (4) the Maxwell-Faraday equation and Eq. (5) Gauss’s law applied at the FERRO-INS interface. In Eq. (4), we assume no work function difference ΔϕTB between top and bottom metal (if ΔϕTB ≠ 0: VTB → VTB − ΔϕTB). In Eq. (5), \({D}_{{\rm{F}}}\,=\,{\varepsilon }_{{\rm{F}}}{{\mathcal{E}}}_{{\rm{F}}}\,+\,{P}_{{\rm{F}}}\) and \({D}_{{\rm{I}}}\,=\,{\varepsilon }_{{\rm{I}}}{{\mathcal{E}}}_{{\rm{I}}}\) are the electric displacement field inside the FERRO and INS layer, with \({{\mathcal{E}}}_{{\rm{I}}}\) and εI the electric field strength and permittivity of the INS, and σFI the surface charge density at the FERRO-INS interface. The y-dependence in Eqs. (4) and (5) enters through the FERRO polarization PF (i.e., we allow for non-zero ∂yPF 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:

$${{\mathcal{E}}}_{{\rm{F}}}\left(y\right)\,=\,\frac{{V}_{{\rm{TB}}}{C}_{{\rm{I}}}\,-\,{P}_{{\rm{F}}}\left(y\right)}{{t}_{{\rm{F}}}\left({C}_{{\rm{F}}}\,+\,{C}_{{\rm{I}}}\right)},$$
( 6)
$${{\mathcal{E}}}_{{\rm{I}}}\left(y\right)\,=\,\frac{{V}_{{\rm{TB}}}{C}_{{\rm{F}}}\,+\,{P}_{{\rm{F}}}\left(y\right)}{{t}_{{\rm{I}}}\left({C}_{{\rm{F}}}\,+\,{C}_{{\rm{I}}}\right)},$$
( 7)

with CF = εF/tF and CI = εI/tI.

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}_{{\rm{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: \({G}_{{\rm{MFIM}}}\,=\,{G}_{{\rm{F}}}\,+\,{G}_{{\rm{I}}}\,=\,{\int}_{{V}_{{\rm{F}}}}{{\mathcal{G}}}_{{\rm{F}}}{\rm{d}}V\,+\,{\int}_{{V}_{{\rm{I}}}}{{\mathcal{G}}}_{{\rm{I}}}{\rm{d}}V\). We now first provide a formal justification that the total Gibbs free energy \({G}_{{\rm{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 Ny thin sheets (discretized in the y-direction). The width of the λth-sheet is \({{{\Delta }}}_{{\rm{y}}}^{\lambda }\), such that: \(L\,=\,{\sum }_{\lambda }{{{\Delta }}}_{{\rm{y}}}^{\lambda }\). We assume that the FERRO and INS are locally uniform (in the x-direction) within each sheet (i.e., for every state variable u: ∂xu = 0), and show that the differential form of the Gibbs free energy becomes:

$${\rm{d}}{G}_{{\rm{MFIM}}} = \, W\mathop{\sum }\limits_{\lambda \,=\,0}^{{N}_{{\rm{y}}}\,-\,1}{{{\Delta }}}_{{\rm{y}}}^{\lambda }\left[{t}_{{\rm{F}}}{\rm{d}}{{\mathcal{G}}}_{{\rm{F}}}^{\lambda }\,+\,{t}_{{\rm{I}}}{\rm{d}}{{\mathcal{G}}}_{{\rm{I}}}^{\lambda }\right],\\ = \, W\mathop{\sum }\limits_{\lambda \,=\,0}^{{N}_{{\rm{y}}}\,-\,1}{{{\Delta }}}_{{\rm{y}}}^{\lambda }\left[-{T}_{{\rm{F}}}\left({t}_{{\rm{F}}}{{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{F}}}^{\lambda }\,+\,{t}_{{\rm{I}}}{{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{I}}}^{\lambda }\right)\,-\,{D}_{{\rm{F}}}^{\lambda }{\rm{d}}{V}_{{\rm{TB}}}\right]$$
( 8)

Because VTB is constrained in the MFIM structure (dVTB = 0), at equilibrium (\({{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{F}}}^{\lambda }\,=\,{{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{I}}}^{\lambda }\,=\,0\)) we have:

$${\rm{d}}{G}_{{\rm{MFIM}}}\,=\,-W\left(\mathop{\sum }\limits_{\lambda \,=\,0}^{{N}_{{\rm{y}}}\,-\,1}{{{\Delta }}}_{{\rm{y}}}^{\lambda }{D}_{{\rm{F}}}^{\lambda }\right){\rm{d}}{V}_{{\rm{TB}}}\,=\,0,$$

which demonstrates that the i-MFIM system evolves to an equilibrium state that minimizes:

$${G}_{{\rm{MFIM}}}= \, W\mathop{\int }\nolimits_{0}^{L}\left({t}_{{\rm{F}}}{{\mathcal{G}}}_{{\rm{F}}}\,+\,{t}_{{\rm{I}}}{{\mathcal{G}}}_{{\rm{I}}}\right){\rm{d}}y,\\ = \, W\mathop{\int }\nolimits_{0}^{L}\left[{t}_{{\rm{F}}}\left(\alpha {\left|{{\bf{P}}}_{{\rm{F}}}\right|}^{2}\,+\,\beta {\left|{{\bf{P}}}_{{\rm{F}}}\right|}^{4}\,+\,{\beta }^{\prime}\mathop{\sum} _{i\ne j}{P}_{{\rm{F}},i}^{2}{P}_{{\rm{F}},j}^{2}\right.\right.\\ +{\kappa }_{{\rm{G}}}\left(\nabla {{\bf{P}}}_{{\rm{F}}}\right)\,-\,{P}_{{\rm{F}},x}{{\mathcal{E}}}_{{\rm{F}}}\,-\,\frac{{\varepsilon }_{{\rm{F}}}}{2}{{\mathcal{E}}}_{{\rm{F}}}^{2}\bigg)\,+\,{t}_{{\rm{I}}}\left(-\frac{{\varepsilon }_{{\rm{I}}}}{2}{{\mathcal{E}}}_{{\rm{I}}}^{2}\right)\bigg]{\rm{d}}y.$$
( 9)

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 PF (i.e., ∇PF = 0). The analysis shows that the equation of state (\({\rm{d}}{G}_{{\rm{MFIM}}}/{\rm{d}}{P}_{{\rm{F}}}=0\)) reproduces the same equilibrium states (i.e., PF-\({{\mathcal{E}}}_{{\rm{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 satisfied14,27:

$${t}_{{\rm{F}}}\,<\,{t}_{{\rm{F,c}}}\,=\,\frac{-1}{2\alpha {C}_{{\rm{I}}}}\left(2\alpha {\varepsilon }_{{\rm{F}}}\,+\,1\right),$$
( 10)

where tF,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 (\({{\mathcal{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 free-charge-conducting MFIM structure (c-MFIM), with leakage through the INS and with identical device dimensions (tF = 3.5 nm and tI = 1 nm), but with surface charge density σFI. The applied sawtooth voltage VTB between the electrodes has a peak voltage amplitude VP (Fig. 1b).

Fig. 1: Free-charge-conducting metal-ferroelectric-insulator-metal (c-MFIM) system with surface charge density σFI at the ferroelectric-insulator (FERRO-INS) interface.
figure 1

a Schematic of the MFIM structure with indicated location of σFI. The FERRO layer, with length L and thickness tF = 3.5 nm, is characterized by a spontaneous polarization density PF and electric field strength \({{\mathcal{E}}}_{{\rm{F}}}\). Similarly, the INS layer, with thickness tI = 1 nm, is characterized by an electric field strength \({{\mathcal{E}}}_{{\rm{I}}}\) and electric permittivity εI = 3.9. b Time signal of the applied input voltage VTB with peak magnitude Vp. c Energy band diagram (VTB = 1 V) showing the Fermi energy level Ef at the top (Ef(0)) and bottom (Ef(tΣ = tF + tI)) metal contact, and the interface charge-neutrality level EN. The interface state density DFI is assumed to be constant across the entire FERRO-INS interface. When the Fermi energy at the FERRO-INS interface is located above (below) EN, σFI is negative (positive).

Full size image

As per the usual treatment of interface traps60,61,62, the surface charge density is expressed as: \({\sigma }_{{\rm{FI}}}\,=\,q{D}_{{\rm{FI}}}\left({E}_{{\rm{N}}}\,-\,{E}_{{\rm{F}}}\left({t}_{{\rm{F}}}\right)\right)\), where DFI is an interface state density in units of eV−1 cm−2, EN is the charge-neutrality level at the FERRO-INS interface, \({E}_{{\rm{F}}}\left(x\right)\) 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 (\({E}_{{\rm{F}}}\left({t}_{{\rm{F}}}\right)\,\approx\, {E}_{{\rm{F}}}\left({t}_{{{\Sigma }}}\right)\)), and we have: \({\sigma }_{{\rm{FI}}}\,\approx\, q{D}_{{\rm{FI}}}\left({E}_{{\rm{N}}}\ \,-\,\ {E}_{{\rm{F}}}\left({t}_{{{\Sigma }}}\right)\right)\). By using \({E}_{{\rm{N}}}\,=\,{{\Delta }}{E}_{{\rm{N}}}\,+\,{E}_{{\rm{F}}}\left({t}_{{{\Sigma }}}\right)\,-\,q{{\mathcal{E}}}_{{\rm{I}}}{t}_{{\rm{I}}}\) (see Fig. 1c), and setting σFI,0 ≡ qDFIΔEN, we obtain:

$${\sigma }_{{\rm{FI}}}\,=\,q{D}_{{\rm{FI}}}\left({{\Delta }}{E}_{{\rm{N}}}\,-\,q{{\mathcal{E}}}_{{\rm{I}}}{t}_{{\rm{I}}}\right)\,=\,{\sigma }_{{\rm{FI,0}}}\,-\,{q}^{2}{D}_{{\rm{FI}}}{{\mathcal{E}}}_{{\rm{I}}}{t}_{{\rm{I}}}.$$
( 11)

For the c-MFIM system the following conditions still apply: \({V}_{{\rm{TB}}}\,=\,{{\mathcal{E}}}_{{\rm{F}}}{t}_{{\rm{F}}}\,+\,{{\mathcal{E}}}_{{\rm{I}}}{t}_{{\rm{I}}}\) (Eq. (4)) and DI = DF + σFI (Eq. (5)) with σFI given by Eq. (11). From this directly follows:

$${{\mathcal{E}}}_{{\rm{F}}}\,=\,\frac{{V}_{{\rm{TB}}}\left({C}_{{\rm{I}}}\,+\,{q}^{2}{D}_{{\rm{FI}}}\right)\,-\,{\sigma }_{{\rm{FI,0}}}\,-\,{P}_{{\rm{F}}}}{{t}_{{\rm{F}}}\left({C}_{{\rm{F}}}\,+\,{C}_{{\rm{I}}}\,+\,{q}^{2}{D}_{{\rm{FI}}}\right)},$$
( 12)
$${{\mathcal{E}}}_{{\rm{I}}}\,=\,\frac{{V}_{{\rm{TB}}}{C}_{{\rm{F}}}\,+\,{\sigma }_{{\rm{FI,0}}}\,+\,{P}_{{\rm{F}}}}{{t}_{{\rm{I}}}\left({C}_{{\rm{F}}}\,+\,{C}_{{\rm{I}}}\,+\,{q}^{2}{D}_{{\rm{FI}}}\right)},$$
( 13)

Note that the electrostatics of the MFIM structure (Eqs. (12) and (13)) reduce to that of a MFM structure for large DFI.

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}_{{\rm{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}_{{\rm{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}_{{\rm{MFIM}}}\), we identify that, at constant applied voltage, the c-MFIM system adopts an equilibrium state that minimizes the Grand Potential \({{{\Omega }}}_{{\rm{MFIM}}}\) of the system, which is defined as63:

$${{\Omega }}\,=\,U\,-\,TS\,-\,\int _{V}\mathop{\sum} _{\nu }{{\mathcal{E}}}_{\nu }{D}_{\nu }{\rm{d}}V\,-\,\mathop{\sum} _{{\rm{j}}}{\mu }_{{\rm{j}}}{N}_{{\rm{j}}},$$
( 14)

with the following associated differential form (Ω = ∫ωdV):

$${\rm{d}}\omega \,=\,-s{\rm{d}}T\,-\,T{{\rm{d}}}_{{\rm{i}}}s\,-\,\mathop{\sum} _{\nu }{D}_{\nu }{\rm{d}}{{\mathcal{E}}}_{\nu }\,+\,\mathop{\sum} _{\alpha }{\mu }_{\alpha }^{\prime}{\rm{d}}{n}_{\alpha }^{\prime}\,-\,\mathop{\sum} _{{\rm{j}}}{n}_{{\rm{j}}}{\rm{d}}{\mu }_{{\rm{j}}},$$
( 15)

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}_{\alpha }^{\prime}\,=\,\int {n}_{\alpha }^{\prime}{\rm{d}}V\)), and all mobile particles (Nj = ∫njdV) 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 \({{{\Omega }}}_{{\rm{MFIM}}}\) is considered to be of the form: \({{{\Omega }}}_{{\rm{MFIM}}}\,=\,{{{\Omega }}}_{{\rm{F}}}\,+\,{{{\Omega }}}_{{\rm{I}}}\,+\,{{{\Omega }}}_{{\rm{FI}}}\,=\,{\int}_{{V}_{{\rm{F}}}}{\omega }_{{\rm{F}}}{\rm{d}}V\,+\,{\int}_{{V}_{{\rm{I}}}}{\omega }_{{\rm{I}}}{\rm{d}}V\,+\,{\int}_{{A}_{{\rm{FI}}}}{\underline{\omega }}_{{\rm{FI}}}{\rm{d}}A\). Here, ΩFI and \({\underline{{{\omega }}}}_{{\rm{FI}}}\) represent the free energy and surface energy density contribution, respectively, to \({{{\Omega }}}_{{\rm{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 \({\rm{d}}{\underline{\omega }}_{FI}\) is of the following form:

$${\rm{d}}\underline{\omega }\,=\,-\underline{s}{\rm{d}}T\,-\,T{{\rm{d}}}_{{\rm{i}}}\underline{s}\,+\,\mathop{\sum} _{\alpha }{\mu }_{\alpha }^{\prime}{\rm{d}}{\underline{n}}_{\alpha }^{\prime}\,-\,\mathop{\sum} _{{\rm{j}}}{\underline{n}}_{{\rm{j}}}{\rm{d}}{\mu }_{{\rm{j}}},$$
( 16)

where \(\underline{s}\) is the entropy surface density, and \({\underline{n}}_{\alpha }^{\prime}\) and \({\underline{n}}_{{\rm{j}}}\) are the surface density of the αth-immobile and jth-mobile particle, respectively.

The differential form of \({{{\Omega }}}_{{\rm{MFIM}}}\) is shown to reduce to (see Supplementary Note 9 for a detailed derivation):

$${\rm{d}}{{{\Omega }}}_{{\rm{MFIM}}}= \,W\mathop{\sum }\limits_{\lambda \,=\,0}^{{N}_{{\rm{y}}}\,-\,1}{{{\Delta }}}_{{\rm{y}}}^{\lambda }\left[\!\!\right.-\!\!{T}_{{\rm{F}}}\left({t}_{{\rm{F}}}{{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{F}}}^{\lambda }\,+\,{t}_{{\rm{I}}}{{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{I}}}^{\lambda }\,+\,{{\rm{d}}}_{{\rm{i}}}{\underline{s}}_{{\rm{FI}}}^{\lambda }\right)\\ -{D}_{{\rm{F}}}^{\lambda }{\rm{d}}{V}_{{\rm{TB}}}\,-\,{t}_{{\rm{I}}}{\sigma }_{{\rm{FI}}}^{\lambda }{\rm{d}}{{\mathcal{E}}}_{{\rm{I}}}^{\lambda }\,-\,{\sigma }_{{\rm{FI}}}^{\lambda }{\rm{d}}{\varphi }_{{\rm{I}}}^{\lambda }\left({t}_{{\rm{F}}}\right)\left]\right..$$
(17)

In contrast to the applied voltage, the electric field in the INS (\({{\mathcal{E}}}_{{\rm{I}}}\)) and the chemical potential at the FERRO-INS interface (\({\varphi }_{{\rm{I}}}\left({t}_{{\rm{F}}}\right)\)) are not explicitly constrained in the c-MFIM system. As a result, it is not self-evident that Eq. (17) becomes zero at equilibrium (\({{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{F}}}^{\lambda }\ \,=\,\ {{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{I}}}^{\lambda }\ \,=\,\ {{\rm{d}}}_{{\rm{i}}}{\underline{s}}_{{\rm{FI}}}^{\lambda }\ \,=\,\ 0\) and dVTB = 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 non-equilibrium thermodynamics47,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 (\(\varphi \left({t}_{{\rm{F}}}\right)\)) and the electric field in the INS (\({{\mathcal{E}}}_{{\rm{I}}}\)) become jointly constrained as well and are, therefore, no longer independent from each other (also see Supplementary Fig. S4):

$${\varphi }_{{\rm{I}}}^{\lambda }\left({t}_{{\rm{F}}}\right)\,=\,{\varphi }_{{\rm{I}}}^{\lambda }\left({t}_{{{\Sigma }}}\right)\,-\,{{\mathcal{E}}}_{{\rm{I}}}^{\lambda }{t}_{{\rm{I}}}.$$
( 18)

This allows us to rewrite Eq. (17) as:

$${\rm{d}}{{{\Omega }}}_{{\rm{MFIM}}} = \, W\mathop{\sum }\limits_{\lambda \,=\,0}^{{N}_{{\rm{y}}}\,-\,1}{{{\Delta }}}_{{\rm{y}}}^{\lambda }\left[-{t}_{{\rm{I}}}{\sigma }_{{\rm{FI}}}^{\lambda }{\rm{d}}{{\mathcal{E}}}_{{\rm{I}}}^{\lambda }\,-\,{\sigma }_{{\rm{FI}}}^{\lambda }{\rm{d}}\left({\varphi }_{{\rm{I}}}^{\lambda }\left({t}_{{{\Sigma }}}\right)\,-\,{{\mathcal{E}}}_{{\rm{I}}}^{\lambda }{t}_{{\rm{I}}}\right)\right]\\ = \, W\mathop{\sum }\limits_{\lambda \,=\,0}^{{N}_{{\rm{y}}}\,-\,1}{{{\Delta }}}_{{\rm{y}}}^{\lambda }\left[-{\sigma }_{{\rm{FI}}}^{\lambda }{\rm{d}}{\varphi }_{{\rm{I}}}^{\lambda }\left({t}_{{{\Sigma }}}\right)\right]\,=\,0.$$
(19)

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:

$${{{\Omega }}}_{{\rm{MFIM}}}= \, W\mathop{\int }\nolimits_{0}^{L}\left({t}_{{\rm{F}}}{\omega }_{{\rm{F}}}\,+\,{t}_{{\rm{I}}}{\omega }_{{\rm{I}}}\,+\,{\underline{\omega }}_{{\rm{FI}}}\right){\rm{d}}y,\\ = \, W\mathop{\int }\nolimits_{0}^{L}\left[{t}_{{\rm{F}}}\left(\alpha {\left|{{\bf{P}}}_{{\rm{F}}}\right|}^{2}\,+\,\beta {\left|{{\bf{P}}}_{{\rm{F}}}\right|}^{4}\,+\,{\beta }^{\prime}\mathop{\sum} _{i\ne j}{P}_{{\rm{F}},i}^{2}{P}_{{\rm{F}},j}^{2} + {\kappa }_{{\rm{G}}}\left(\nabla {{\bf{P}}}_{{\rm{F}}}\right)\right.\right.\\ -\left.\left.{P}_{{\rm{F}},x}{{\mathcal{E}}}_{{\rm{F}}}\,-\,\frac{{\varepsilon }_{{\rm{F}}}}{2}{{\mathcal{E}}}_{{\rm{F}}}^{2}\right)\,\!-\,\!\frac{{\varepsilon }_{{\rm{I}}}\,+\,{q}^{2}{D}_{{\rm{FI}}}{t}_{{\rm{I}}}}{2}{t}_{{\rm{I}}}{{\mathcal{E}}}_{{\rm{I}}}^{2}\,+\,{\sigma }_{{\rm{FI,0}}}{t}_{{\rm{I}}}{{\mathcal{E}}}_{{\rm{I}}}\right]{\rm{d}}y.$$
(20)

See Supplementary Note 10 for a detailed derivation of the expression above. Observe that Eq. (20) explicitly depends on the interface state density DFI: in Figs. 2a, b, 3a we show that increasing DFI deforms the Grand Potential energy density landscape of the system. Note that Eq. (20) is equivalent to Eq. (9) in the limit of DFI = 0, which demonstrates that a Grand Potential-based approach is indistinguishable from a Gibbs-based approach for the i-MFIM system.

Fig. 2: Surface charge density at the ferroelectric-insulator interface destabilizes the Grand Potential energy landscape of the free-charge-conducting metal-ferroelectric-insulator-metal (c-MFIM) system.
figure 2

The Grand Potential energy (\({\underline{{{\Omega }}}}_{{\rm{MFIM}}}\,=\,{{{\Omega }}}_{{\rm{MFIM}}}/A\)) landscape is plotted as a function of the ferroelectric polarization PF and the applied voltage VTB for an interface state density (a) DFI = 1 × CI/q2 and (b) DFI = 5 × CI/q2, where CI is the capacitance of the insulator layer and q is the electron charge. (Note that the free energy landscape shown in Supplementary Fig. S3a corresponds to the case of DFI = 0). a, b The stable, metastable, and unstable equilibrium states of the c-MFIM system (at constant applied voltage VTB) are indicated for both cases. In addition, the trajectory (a\(\overline{{\rm{a}}}\)) along the energy landscape corresponding to a constant VTB = 0 V is indicated. Comparing the energy landscapes of a, b demonstrates that increasing the interface state density DFI deforms the Grand Potential energy landscape, resulting in an increased hysteretic device operation of the c-MFIM system.

Full size image
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.
figure 3

a The Grand Potential energy (\({\underline{{{\Omega }}}}_{{\rm{MFIM}}}\)) landscape of the c-MFIM system (see (Fig. 1a) is plotted along the trajectory (a\(\overline{{\rm{a}}}\)) of constant applied voltage VTB = 0 V for different values of the interface state density DFI. The constant-VTB trajectories for DFI = 1 × CI/q2 and DFI = 5 × CI/q2, where CI 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 DFI = 0, the c-MFIM system has a single stable state with ferroelectric polarization PF = 0. However, increasing the interface state density DFI deforms the Grand Potential energy landscape. As a consequence, for larger DFI, the c-MFIM system is characterized by two distinct stable states with non-zero polarization (\({P}_{{\rm{F,0}}}^{1}\) and \({P}_{{\rm{F,0}}}^{2}\)). b Quasi-static PF-t simulation results show that increasing DFI results in abrupt switching of the ferroelectric polarization PF. (The dashed and solid curves correspond to quasi-static simulations with different initial polarization states (\({P}_{{\rm{F,0}}}^{1}\) or \({P}_{{\rm{F,0}}}^{2}\)) of the system). c The PF-\({{\mathcal{E}}}_{{\rm{F}}}\) characteristics, where \({{\mathcal{E}}}_{{\rm{F}}}\) is the electric field strength in the ferroelectric layer, and d ΔΨI-VTB 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 DFI. The case of DFI = 0 corresponds to the intrinsic MFIM system (Supplementary Fig. S3) and DFI → ∞ to the metal-ferroelectric-metal (MFM) system (Supplementary Fig. S1).

Full size image

Next, we utilize \({{{\Omega }}}_{{\rm{MFIM}}}\) to analyze the equilibrium behavior of the c-MFIM system in the case of a coherent uniform FERRO layer with uniform PF (i.e., ∇PF = 0) and constant interface state density (DFI) across the entire FERRO-INS interface surface area. The equilibrium states of the FERRO-system are ascertained from the equation of state (\({\rm{d}}{{{\Omega }}}_{{\rm{MFIM}}}/{\rm{d}}{P}_{{\rm{F}}}\ =0\)):

$$ 2 \alpha {t}_{{\rm{F}}}{P}_{{\rm{F}}}\,+\,4\beta {t}_{{\rm{F}}}{P}_{{\rm{F}}}^{3}\,+\,\frac{{P}_{{\rm{F}}}}{\left({C}_{{\rm{F}}}\,+\,{C}_{{\rm{I}}}\,+\,{q}^{2}{D}_{{\rm{FI}}}\right)}\\ -\,\frac{{V}_{{\rm{TB}}}\left({C}_{{\rm{I}}}\,+\,{q}^{2}{D}_{{\rm{FI}}}\right)\,-\,{\sigma }_{{\rm{FI,0}}}}{\left({C}_{{\rm{F}}}\,+\,{C}_{{\rm{I}}}\,+\,{q}^{2}{D}_{{\rm{FI}}}\right)}\,=\,0.$$
( 21)

Using Eq. (12), we retrieve Eq. (3) from Eq. (21), which, once again, demonstrates that the PF-\({{\mathcal{E}}}_{{\rm{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 \({{\rm{d}}}^{2}{{{\Omega }}}_{{\rm{MFIM}}}/{\rm{d}}{P}_{{\rm{F}}}^{2}\,> \,0\):

$$12\beta {t}_{{\rm{F}}}{P}_{{\rm{F}}}^{2}\,+\,2\alpha {t}_{{\rm{F}}}\,+\,\frac{1}{\left({C}_{{\rm{F}}}\,+\,{C}_{{\rm{I}}}\,+\,{q}^{2}{D}_{{\rm{FI}}}\right)}\,> \,0.$$
( 22)

The impact of the interface trap density DFI 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:

$${t}_{{\rm{F}}}\,<\,{t}_{{\rm{F,c}}}\,=\,\frac{-1}{2\alpha \left({C}_{{\rm{I}}}\,+\,{q}^{2}{D}_{{\rm{FI}}}\right)}\left(2\alpha {\varepsilon }_{{\rm{F}}}\,+\,1\right).$$
( 23)

Note that, for DFI = 0, Eq. (23) reduces to Eq. (10).

Figure 3 illustrates the impact of an increasing DFI on the quasi-static simulation results for the “time”-dependent signal of the polarization \({P}_{{\rm{F}}}\left(t\right)\) (Fig. 3b), the PF-\({{\mathcal{E}}}_{{\rm{F}}}\) characteristics (Fig. 3c), and the dependence of the electrostatic potential difference over the INS ΔΨI on the applied voltage VTB, 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 DFI = 0, and to the MFM system in the limit of DFI → ∞). In addition, the results shown in Fig. 3b–d show that, all else being equal, an increase in DFI 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 DFI is in agreement with theoretical predictions using a kinetics-based model for the polarization switching65, and the experimental observation that the memory window of a device increases with the trap density66.

In Fig. 4, the critical FERRO thickness tF,c is plotted as a function of DFI and clearly shows how the stabilization condition of the c-MFIM system is negatively impacted by an increasing DFI. The sensitivity of tF to values of α, CI, and εF is also shown in Fig. 4, demonstrating that, all else being equal, tF,c decreases with increasing values of α, CI, and εF, which is consistent with earlier reports27. Note that, for DFI → ∞, Eq. (23) tends to \({t}_{{\rm{F}}}\,<\,{t}_{{\rm{F,c}}}\,=\,-\left(1\,+\,2{\varepsilon }_{{\rm{F}}}\alpha \,+\,\sqrt{1\,+\,2{\varepsilon }_{{\rm{F}}}\alpha }\right)/\left(2\alpha {q}^{2}{D}_{{\rm{FI}}}\right)\), which shows that the impact of CI becomes negligible for larger values of DFI and reflects the transition from a MFIM to a MFM system. Lastly, note that tF,c also depends on the temperature, as demonstrated in Supplementary Fig. S6.

Fig. 4: The critical ferroelectric thickness decreases with increasing interface state density.
figure 4

The critical ferroelectric thickness tF,c of the free-charge-conducting metal-ferroelectric-insulator-metal (c-MFIM) system (Fig. 1a) is plotted as a function of the interface state density DFI for different values of the Landau parameter α (yellow lines), the background permittivity of the ferroelectric εF (blue lines), and the capacitance of the insulator layer CI = εI/tI (red lines). The baseline curve (black line) corresponds to the parameter values: α = −1.1 × 109 mF−1, εF = 35, and CI = 3.45 μF cm−2. The c-MFIM system is characterized by a hysteretic device operation if the thickness tF of the (perfectly uniform and coherent) ferroelectric layer is larger than the critical thickness tF,c (tF > tF,c). Note that the inflection point in the curves occurs around DFI = CI/q2.

Full size image

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 \({{{\Omega }}}_{{\rm{MFIM}}}\) on DFI. The dependence of \({{{\Omega }}}_{{\rm{MFIM}}}\) on DFI is traced to the impact of DFI on \({{\mathcal{E}}}_{{\rm{F}}}\) (Eq. (12)) and \({{\mathcal{E}}}_{{\rm{I}}}\) (Eq. (13)). An increasing DFI more strongly affects \({{\mathcal{E}}}_{{\rm{I}}}\) (and therefore ΩI) than \({{\mathcal{E}}}_{{\rm{F}}}\) (and therefore ΩF), resulting in a proportionally reduced stabilizing contribution of ΩI to \({{{\Omega }}}_{{\rm{MFIM}}}\). As a consequence, a smaller CI(=εI/tI) or larger CF(=εF/tF), which, in turn, disproportionately affects \({{\mathcal{E}}}_{{\rm{F}}}\) over \({{\mathcal{E}}}_{{\rm{I}}}\), is required to prevent the formation of unstable regimes in the Grand Potential energy landscape.

Lastly, Eq. (23) also allows us to establish the theoretical maximum value of DFI for which a hysteresis-free device operation of a MFIM system with a uniform single-domain ferroelectric layer is possible:

$${D}_{{\rm{FI}}}\,<\,{D}_{{\rm{FI,c}}}\,=\,\frac{-1}{2\alpha {q}^{2}{t}_{{\rm{F}}}}\left[2\alpha {t}_{{\rm{F}}}\left({C}_{{\rm{F}}}\,+\,{C}_{{\rm{I}}}\right)\,+\,1\right]$$
( 24)

For the c-MFIM system considered in this Article (Fig. 1a), DFI,c ≈ 2 × 1012 eV−1cm−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 tS (and area A), and with pristine FERRO-SEMI interface (DFS = 0). In contrast to the INS and FERRO layer in the MFIM system, at finite temperatures, it is no longer acceptable to neglect the thermally-generated electron (nS) and hole (pS) free carrier concentrations in the SEMI layer. The local free charge density, \({\rho }_{{\rm{S}}}\,=\,q\left({p}_{{\rm{S}}}\,-\,{n}_{{\rm{S}}}\right)\), affects the local electric displacement field (DS) in the SEMI through Gauss’ law: \(\partial {D}_{{\rm{S}}}\left(x\right)/\partial x\,=\,{\rho }_{{\rm{S}}}\left(x\right)\), where \({D}_{{\rm{S}}}\,=\,{\varepsilon }_{{\rm{S}}}{{\mathcal{E}}}_{{\rm{S}}}\) with \({{\mathcal{E}}}_{{\rm{S}}}\) the local electric field strength and εS the permittivity of the SEMI layer. Because of the position dependent (x-direction) electric field \({{\mathcal{E}}}_{{\rm{S}}}\) and electric displacement field DS, 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 VTB, 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 approximation47,64 and consider, within each λth thin sheet, the SEMI to consist of Nx elemental volumes \({\rm{d}}{V}_{{\rm{S}}}^{\iota ,\lambda }\), 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: \({\rm{d}}{V}_{{\rm{S}}}^{\iota ,\lambda }\,=\,{{{\Delta }}}_{{\rm{x}}}^{\iota }\,\times\, {{{\Delta }}}_{{\rm{y}}}^{\lambda }\,\times\, W\). Each elemental volume is then characterized by a locally uniform electric field \({{\mathcal{E}}}_{{\rm{S}}}^{\iota ,\lambda }\), electric displacement field \({D}_{{\rm{S}}}^{\iota ,\lambda }\), and free charge density \({\rho }_{{\rm{S}}}^{\iota ,\lambda }\,=\,q({p}_{{\rm{S}}}^{\iota ,\lambda }\,-\,{n}_{{\rm{S}}}^{\iota ,\lambda })\). As a result, Eq. (15) is still valid for each elemental volume \({\rm{d}}{V}_{{\rm{S}}}^{\iota ,\lambda }\) of the SEMI, such that the differential form of the Grand Potential of the i-MFS system (ΩMFS = ΩF + ΩS) is found as:

$${\rm{d}}{{{\Omega }}}_{{\rm{MFS}}}\,=\,W\mathop{\sum }\limits_{\lambda \,=\,0}^{{N}_{{\rm{y}}}\,-\,1}{{{\Delta }}}_{{\rm{y}}}^{\lambda }\left[{t}_{{\rm{F}}}{\rm{d}}{\omega }_{{\rm{F}}}^{\lambda }\,+\,\mathop{\sum }\limits_{\iota \,=\,0}^{{N}_{{\rm{x}}}\,-\,1}{{{\Delta }}}_{{\rm{x}}}^{\iota }{\rm{d}}{\omega }_{{\rm{S}}}^{\iota ,\lambda }\right].$$
( 25)

With the requirement of a negligible steady-state leakage current through the system, we show that, at equilibrium (\({{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{F}}}^{\lambda }\,=\,{{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{S}}}^{\iota ,\lambda }\,=\,0\)), Eq. (25) reduces to (see Supplementary Note 11):

$${\rm{d}}{{{\Omega }}}_{{\rm{MFS}}}\,=\,W\mathop{\sum }\limits_{\lambda \,=\,0}^{{N}_{{\rm{y}}}\,-\,1}{{{\Delta }}}_{{\rm{y}}}^{\lambda }\left[-{D}_{{\rm{F}}}^{\lambda }{\rm{d}}{V}_{{\rm{TB}}}\right],$$
( 26)

which demonstrates that, at constant applied voltage (dVTB = 0), the MFS system evolves to an equilibrium state that minimizes the total Grand Potential ΩMFS of the system:

$${{{\Omega }}}_{{\rm{MFS}}}= \, W\mathop{\int }\nolimits_{0}^{L}\left[\right.{t}_{{\rm{F}}}\left(\right.\alpha {\left|{{\bf{P}}}_{{\rm{F}}}\right|}^{2}\,+\,\beta {\left|{{\bf{P}}}_{{\rm{F}}}\right|}^{4}\,+\,{\beta }^{\prime}\mathop{\sum} _{i\,\ne\, j}{P}_{{\rm{F}},i}^{2}{P}_{{\rm{F}},j}^{2}\\ +{\kappa }_{{\rm{G}}}\left(\nabla {{\bf{P}}}_{{\rm{F}}}\right)\,-\,{P}_{{\rm{F}},x}{{\mathcal{E}}}_{{\rm{F}}}\,-\,\frac{{\varepsilon }_{{\rm{F}}}}{2}{{\mathcal{E}}}_{{\rm{F}}}^{2}\left)\right.\\ +\mathop{\int }\nolimits_{0}^{{t}_{{\rm{S}}}}\left(-\frac{{\varepsilon }_{{\rm{S}}}}{2}{{\mathcal{E}}}_{{\rm{S}}}^{2}\,-\,\mathop{\int }\nolimits_{{\varphi }_{{\rm{S,0}}}}^{{\varphi }_{{\rm{S}}}}{\rho }_{{\rm{S}}}{\rm{d}}{\varphi }_{{\rm{S}}}\right){\rm{d}}x\left]\right.{\rm{d}}y,$$
(27)

where φS,0 is the zero-field (\({{\mathcal{E}}}_{{\rm{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-insulator-semiconductor (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 (DFI). 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 (DIS) at the INS-SEMI interface. In the MFIS system, the following conditions apply:

$${V}_{{\rm{TB}}}\,=\,{{\mathcal{E}}}_{{\rm{F}}}(y){t}_{{\rm{F}}}\,+\,{{\mathcal{E}}}_{{\rm{I}}}(y){t}_{{\rm{I}}}\,+\,{{\Delta }}{\Psi }_{{\rm{S}}},$$
( 28)
$${D}_{{\rm{I}}}(y)\,=\,{D}_{{\rm{F}}}(y)\,+\,{\sigma }_{{\rm{FI}}}(y),$$
( 29)
$${D}_{{\rm{S}}}(y)\,=\,{D}_{{\rm{I}}}(y)\,+\,{\sigma }_{{\rm{IS}}}(y)\,+\,{\sigma }_{{\rm{S}}}(y).$$
( 30)

Where Eq. (28) again corresponds to the Maxwell-Faraday equation (note that the left-hand side becomes: VTB − ΔϕTS, for non-zero work function difference, ΔϕTS, between top metal and SEMI), and Eqs. (29) and (30) are Gauss’s law applied at the FERRO-INS and INS-SEMI interface, respectively, with σS the 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 = qDFI(EN,FI − EF(tF)) and σIS = qDIS(EN,IS − EF(tF + tI)). 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:

$${\sigma }_{{\rm{FI}}}\,=\,{\sigma }_{{\rm{FI,0}}}\,-\,{q}^{2}{D}_{{\rm{FI}}}\left({{\mathcal{E}}}_{{\rm{I}}}{t}_{{\rm{I}}}\,+\,{{\Delta }}{\Psi }_{{\rm{S}}}\right),$$
( 31)
$${\sigma }_{{\rm{IS}}}\,=\,{\sigma }_{{\rm{IS,0}}}\,-\,{q}^{2}{D}_{{\rm{IS}}}{{\Delta }}{\Psi }_{{\rm{S}}},$$
( 32)

with σFI,0 = qDFIΔEN,FI and σIS,0 = qDISΔEN,IS.

We demonstrate that the Grand Potential (ΩMFIS) is the appropriate thermodynamic potential of the MFIS system at constant applied voltage VTB. Analogous to the MFS system, we use the local equilibrium approximation47,64 and consider the SEMI layer, within the λth thin sheet, to consist of Nx locally uniform elemental volumes, discretized in the x-direction.

The Grand Potential of the MFIS system is of the form: \({{{\Omega }}}_{{\rm{MFIS}}}\,=\,{{{\Omega }}}_{{\rm{F}}}\,+\,{{{\Omega }}}_{{\rm{I}}}\,+\,{{{\Omega }}}_{{\rm{FI}}}\,+\,{{{\Omega }}}_{{\rm{IS}}}\,=\,{\int}_{{V}_{{\rm{F}}}}{\omega }_{{\rm{F}}}{\rm{d}}V\,+\,{\int}_{{V}_{{\rm{I}}}}{\omega }_{{\rm{I}}}{\rm{d}}V\,+{\int}_{{A}_{{\rm{FI}}}}{\underline{\omega }}_{{\rm{FI}}}{\rm{d}}A\,+\,{\int}_{{A}_{{\rm{IS}}}}{\underline{\omega }}_{{\rm{IS}}}{\rm{d}}A\). 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:

$${\rm{d}}{{{\Omega }}}_{{\rm{MFIS}}}\,=\,W\mathop{\sum }\limits_{\lambda \,=\,0}^{{N}_{{\rm{y}}}\,-\,1}{{{\Delta }}}_{{\rm{y}}}^{\lambda }\left[{t}_{{\rm{F}}}{\rm{d}}{\omega }_{{\rm{F}}}^{\lambda }\,+\,{t}_{{\rm{I}}}{\rm{d}}{\omega }_{{\rm{I}}}^{\lambda }\,+\,\mathop{\sum }\limits_{\iota \,=\,0}^{{N}_{{\rm{x}}}\,-\,1}\left(\right.{{{\Delta }}}_{{\rm{x}}}^{\iota }{\rm{d}}{\omega }_{{\rm{S}}}^{\iota ,\lambda }\left)\right.\,+\,{\rm{d}}{\underline{\omega }}_{{\rm{FI}}}^{\lambda }\,+\,{\rm{d}}{\underline{\omega }}_{{\rm{IS}}}^{\lambda }\right].$$
( 33)

At equilibrium (\({{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{F}}}^{\lambda }\,=\,{{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{I}}}^{\lambda }={{\rm{d}}}_{{\rm{i}}}{s}_{{\rm{S}}}^{\iota ,\lambda }\,=\,{{\rm{d}}}_{{\rm{i}}}{\underline{s}}_{{\rm{FI}}}^{\lambda }\,=\,{{\rm{d}}}_{{\rm{i}}}{\underline{s}}_{{\rm{IS}}}^{\lambda }\,=\,0\)), 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 \({\rm{d}}{\underline{\omega }}_{{\rm{FI}}}\) and \({\rm{d}}{\underline{\omega }}_{{\rm{IS}}}\), we demonstrate that, at equilibrium, dΩMFIS reduces to (see Supplementary Note 13 for a detailed derivation):

$${\rm{d}}{{{\Omega }}}_{{\rm{MFIS}}}\,=\,W\mathop{\sum }\limits_{\lambda \,=\,0}^{{N}_{{\rm{y}}}\,-\,1}{{{\Delta }}}_{{\rm{y}}}^{\lambda }\left[\right.-{D}_{{\rm{F}}}^{\lambda }{\rm{d}}{V}_{{\rm{TB}}}\left]\right..$$
( 34)

which demonstrates that, at constant applied voltage (dVTB = 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):

$${{{\Omega }}}_{{\rm{MFS}}}= W\mathop{\int }\nolimits_{0}^{L}\left[\right.{t}_{{\rm{F}}}\left(\right.\alpha {\left|{{\bf{P}}}_{{\rm{F}}}\right|}^{2}\,+\,\beta {\left|{{\bf{P}}}_{{\rm{F}}}\right|}^{4}\,+\,{\beta }^{\prime}\mathop{\sum} _{i\ne j}{P}_{{\rm{F}},i}^{2}{P}_{{\rm{F}},j}^{2}\,+\,{\kappa }_{{\rm{G}}}\left(\nabla {{\bf{P}}}_{{\rm{F}}}\right)\\ -{P}_{{\rm{F}},x}{{\mathcal{E}}}_{{\rm{F}}}\,-\,\frac{{\varepsilon }_{{\rm{F}}}}{2}{{\mathcal{E}}}_{{\rm{F}}}^{2}\left)\right.\,-\,\frac{{\varepsilon }_{{\rm{I}}}{t}_{{\rm{I}}}}{2}{{\mathcal{E}}}_{{\rm{I}}}^{2}\,-\,\mathop{\int }\nolimits_{0}^{{t}_{{\rm{S}}}}\left(\frac{{\varepsilon }_{{\rm{S}}}}{2}{{\mathcal{E}}}_{{\rm{S}}}^{2}\,+\,\mathop{\int }\nolimits_{{\varphi }_{{\rm{S,0}}}}^{{\varphi }_{{\rm{S}}}}{\rho }_{{\rm{S}}}{\rm{d}}{\varphi }_{{\rm{S}}}\right){\rm{d}}x\\ +\mathop{\int }\nolimits_{0}^{{{\mathcal{E}}}_{{\rm{I}}}}{\sigma }_{{\rm{FI}}}{\rm{d}}\left({{\mathcal{E}}}_{{\rm{I}}}{t}_{{\rm{I}}}\right)\,+\,\mathop{\int }\nolimits_{0}^{{{\Delta }}{\Psi }_{{\rm{S}}}}{\sigma }_{{\rm{FI}}}{\rm{d}}\left({{\Delta }}{\Psi }_{{\rm{S}}}\right)\,+\,\mathop{\int }\nolimits_{0}^{{{\Delta }}{\Psi }_{{\rm{S}}}}{\sigma }_{{\rm{IS}}}{\rm{d}}\left({{\Delta }}{\Psi }_{{\rm{S}}}\right)\left]\right.{\rm{d}}y.$$
(35)

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 (first-order approximation) stabilization condition of the MFIS system:

$${t}_{{\rm{F}}}\,<\,{t}_{{\rm{F,c}}}\,=\,\frac{-1}{2\alpha {C}_{{{\Sigma }}}^{\dagger }}\left(2\alpha {\varepsilon }_{{\rm{F}}}\,+\,1\right),$$
( 36)

where \({C}_{{{\Sigma }}}^{\dagger }\) is defined as:

$${C}_{{{\Sigma }}}^{\dagger }\,=\,{q}^{2}{D}_{{\rm{FI}}}\,+\,{C}_{{\rm{I}}}\frac{{C}_{{\rm{S}}}\,+\,2q\ {n}_{{\rm{i}}}{t}_{{\rm{S}}}/{V}_{{\rm{th}}}\,+\,{q}^{2}{D}_{{\rm{IS}}}}{{C}_{{\rm{I}}}\,+\,{C}_{{\rm{S}}}\,+\,2q\ {n}_{{\rm{i}}}{t}_{{\rm{S}}}/{V}_{{\rm{th}}}\,+\,{q}^{2}{D}_{{\rm{IS}}}}.$$
( 37)

The stabilization condition (Eq. (36)) indicates that the critical FERRO thickness decreases with increasing interface trap density at both the FERRO-INS interface (DFI) and INS-SEMI interface (DIS), and with increasing free charge density in the SEMI (ni).

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 DFI at the FERRO-INS interface destabilizes MFIM systems. Furthermore, we have identified the theoretical upper limit for the critical interface trap density DFI,c in MFIM structures, above which (increased) hysteretic device operation becomes unavoidable. As DFI,c has been derived assuming perfectly uniform single-domain FERRO-systems, DFI,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 \({\bf{P}}_{\mathrm{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 VTB with quasi-static sawtooth profile, as shown in Supplementary Fig. S4b, with peak voltage amplitude VP = 0.5 V. We discretize VTB in equidistant steps with increments of ΔV (typically 0.5 mV). For each value of VTB, the quasi-static FERRO polarization \({P}_{{\rm{F}}}\left({V}_{{\rm{TB}}}\right)\) is determined from the equation of state of the FERRO-system: dΦ/dPF = 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 \({\underline{G}}_{{\rm{MFIM}}}\) (see Eq. (9)), and, in the case of the c-MFIM system, Φ is the Grand Potential \({\underline{{{\Omega }}}}_{{\rm{MFIM}}}\) (see Eq. (20)). In a first step, we use MATLAB’s fzero algorithm to solve for all existing real-valued roots \({\left\{{P}_{{\rm{F}},1\le i\le N}\right\}}_{N\le3}\) to the cubic equation of state (dΦ/dPF = 0). In a second step, we establish the stability of each real-valued root using the corresponding stability condition: \({{\rm{d}}}^{2}{{\Phi }}/{\rm{d}}{P}_{{\rm{F}}}^{2}\;> \;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 (PF,1) remains:

$${P}_{{\rm{F}}}\left({V}_{{\rm{TB}}}\right)\,=\,{P}_{{\rm{F,1}}},$$

and b) in case of both a stable (PF,1) and metastable (PF,2) root:

$${P}_{{\rm{F}}}\left({V}_{{\rm{TB}}}\right)\,=\,\min \left(\left|{P}_{{\rm{F}}}\left({V}_{{\rm{TB}}}\,-\,{{\Delta }}V\right)\,-\,{P}_{{\rm{F,1}}}\right|,\left|{P}_{{\rm{F}}}\left({V}_{{\rm{TB}}}\,-\,{{\Delta }}V\right)\,-\,{P}_{{\rm{F}},2}\right|\right),$$

with \({P}_{{\rm{F}}}\left({V}_{{\rm{TB}}}\,-\,{{\Delta }}V\right)\) the FERRO polarization of the previous (time) step. At t = 0 s, \({P}_{{\rm{F}}}\left({V}_{{\rm{TB}}}\,-\,{{\Delta }}V\right)\) corresponds to the assumed initial polarization PF,0 of the FERRO-system (see, for example, Fig. 3a). As a last step, the electric field strength in the FERRO (\({{\mathcal{E}}}_{{\rm{F}}}\)) and the INS (\({{\mathcal{E}}}_{{\rm{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).