A desktop app for simulating and analyzing ferroelectric P-E hysteresis loops using Landau-Ginzburg-Devonshire theory — built to help researchers identify nonidealities (leakage, imprint, dead layers) in thin-film ferroelectric capacitors.
Author: Rohan S. Pankaj
Advisor: Prof. Lucas Caretta
Engineering Senior Capstone — Brown University
Ferroelectric materials maintain spontaneous electric polarization without an external field. Applying a field switches this polarization, enabling multiple stable states — the basis of technologies like FeRAM. These devices rely on ferroelectric capacitors (FeCAPs), where a ferroelectric film replaces the dielectric in a standard capacitor. Hysteresis measurements of FeCAPs visualize this switching behavior and expose key material properties: saturation polarization (Ps), remanent polarization (PR), and coercive field (Ec).
In practice, measured loops are distorted by nonidealities — leakage current, built-in fields (imprint), dead layers, and dielectric tilt. FeSim is open-source software that simulates ferroelectric hysteresis and these nonidealities, designed to help researchers identify sources of distortion in their devices and inform material growth and iteration.
- Forward simulation — generate P-E and P-V loops from Landau-Ginzburg-Devonshire theory for HZO, PZT, or BaTiO₃
- Parasitic effects — dielectric response, ohmic leakage, dead layers, and imprint
- Experimental analysis — load a probe-station CSV (PIEC-compatible), extract P_r, V_c, loop area, imprint, and effective permittivity
- Landau model overlay — superimpose a physics model on experimental data with adjustable dead layer and imprint parameters
- Material database — literature-calibrated parameters for HZO, PZT, BTO; substrates SRO, STO, Si, MgO, TiN; electrodes Pt, TiN, W, YBCO, SRO
cd ferroelectric_app
pip install -r requirements.txt
python main.pyThe simulation follows the methodology outlined in Chandra et al. [1].
The free energy density of the ferroelectric is modeled as a Landau-Ginzburg-Devonshire expansion in the out-of-plane polarization P:
where E is the applied electric field. Minimizing with respect to P (i.e. setting ∂𝒻_P/∂P = 0) yields the constitutive relation used in the simulation:
This is a fifth-degree polynomial in P, solved numerically at each voltage step using scipy.fsolve. The double-well shape of the free energy produces two stable polarization states — the physical origin of hysteresis.
The three Landau coefficients each capture a distinct physical effect.
Coefficient a — temperature sensitivity:
where T₀ is the Curie-Weiss temperature and a₀ is a material-specific proportionality constant. For T < T₀, a < 0, giving a double-well (ferroelectric phase). For T > T₀, a > 0, giving a single-well (paraelectric phase).
Coefficient b — fourth-order stiffness:
Controls the magnitude of spontaneous polarization. Its sign determines the order of the phase transition: b < 0 is first-order (BaTiO₃), b > 0 is second-order (HZO, PZT).
Coefficient c — stabilization term:
Ensures the free energy is bounded from below. Always positive.
Bulk Landau coefficients must be corrected for epitaxial thin films, where biaxial strain from lattice mismatch with the substrate alters the ferroelectric behavior. The misfit strain is:
where
The coefficient c is unaffected by strain. The misfit strain also shifts the effective Curie-Weiss temperature:
The final equation solved by the simulation is:
Equation (9) is multivalued in the switching region — for a given E there are three solutions for P, of which two are stable (the upper and lower branches of the loop). The simulation traces the physically realized metastable curve:
- Upper branch — the solver follows the high-P root until V falls below the negative coercive voltage Vc⁻
- Lower branch — the solver follows the low-P root until V rises above the positive coercive voltage Vc⁺
At each step scipy.fsolve solves V(P) - V_target = 0 with the current branch as the initial guess. When the voltage crosses a coercive point, the initial guess is switched to the opposite branch and the solver finds the new stable solution.
The coercive voltages are found analytically from the spinodal condition dV/dP = 0:
which is a quadratic in P² with roots at the switching polarizations.
Real loops are distorted by currents beyond the ideal switching current. These are modeled by augmenting the total current:
The non-switching dielectric background of the ferroelectric contributes:
where d is the film thickness and εr is the background permittivity. Integrating this current into polarization produces a linear tilt of the P-V loop.
Leakage current is calculated from the applied voltage and a user-specified leakage resistance:
Integrating this into the polarization fattens the loop and inflates the apparent loop area without contributing to true switchable polarization.
A non-ferroelectric interfacial layer (thickness tdl, permittivity εdl) acts as a series capacitor, taking a fraction of the applied voltage:
This tilts and narrows the loop, reducing apparent Pr and Ec.
Asymmetric charge trapping at one interface generates a built-in field that shifts the loop horizontally. In the simulation, this is applied as:
Imprint is extracted from data as the average of the positive and negative coercive voltages:
Users can manually correct leakage in the analysis tab by specifying a leakage resistance. The estimated leakage current is calculated and subtracted from the current response. Automatic leakage extraction from multi-frequency measurements (following Meyer et al. [2]) is a planned future feature.
From a measured or simulated P-V loop:
| Parameter | Definition | Physical meaning |
|---|---|---|
| Pr⁺, Pr⁻ | Polarization at V = 0 | Remanent (stored) polarization |
| Vc⁺, Vc⁻ | Voltage at P = 0 | Coercive voltages |
| Ps | (Pmax − Pmin) / 2 | Saturation polarization estimate |
| Vimprint | (Vc⁺ + Vc⁻) / 2 | Built-in bias from interface asymmetry |
| Loop area | ∮ P dV | Energy dissipated per cycle (μJ/cm²) |
| ε_eff | (dP/dV) · t / ε₀ at Vc | Effective permittivity near switching |
Coefficients are taken from literature on bulk properties and corrected for thin films per equations (5)–(8).
| Parameter | HZO | PZT (53/47) | BaTiO₃ | Units |
|---|---|---|---|---|
| a₀ | 4.0×10⁶ | 3.1×10⁵ | 4.1×10⁵ | J·m/(C²·K) |
| b | 4.74×10¹⁰ | 1.91×10⁸ | −2.1×10⁸ | J·m⁵/C⁴ |
| c | 1.0×10¹⁰ | 1.40×10⁹ | 1.3×10⁹ | J·m⁹/C⁶ |
| T₀ | 773 K | 673 K | 393 K | K |
| Q₁₂ | −0.02 | −0.046 | −0.034 | m⁴/C² |
| εr | 25 | 500 | 1500 | — |
| Film thickness | 10 nm | 40 nm | 100 nm | — |
BTO has b < 0, making it a first-order (discontinuous) ferroelectric phase transition. HZO coefficients are calibrated to Tian et al. [9].
Python was chosen for its scientific computing ecosystem (NumPy, SciPy, Matplotlib, Pandas) and suitability for rapid prototyping in research environments. scipy.fsolve is used to numerically solve equation (9), avoiding the need for a closed-form analytical solution and preserving computational resources.
LGD model over alternatives — the fifth-order polynomial is computationally lightweight, well-suited for a desktop app, and physically intuitive. The key limitation is that it cannot capture frequency-dependent switching dynamics.
Fitting is deliberately constrained to a₀, b, c, and T₀ only, to prevent overfitting when dozens of parameters are present.
PIEC compatibility — data files accepted by FeSim follow the CSV structure produced by PIEC, an open-source ferroelectric testing tool. The shared Python/Tkinter foundation enables future integration.
ferroelectric_app/
├── main.py # Entry point
├── fe_model.py # LGD physics engine
├── gui.py # CustomTkinter GUI (Simulation + Analysis tabs)
├── material_parameters.py # Material, substrate, and electrode databases
└── requirements.txt
- Automatic fitting of Landau coefficients to experimental data
- Multi-file upload for frequency-dependent leakage extraction (following Meyer et al. [2])
Portions of this codebase were drafted with the assistance of generative AI tools (Claude, GitHub Copilot). All generated code was reviewed, tested, and validated by the author. All physics, design decisions, and scientific content are the author's own work.
[1] P. Chandra and P. B. Littlewood, "A Landau Primer for Ferroelectrics," in Physics of Ferroelectrics, Springer, Berlin, 2007.
[2] R. Meyer and R. Waser, "Dynamic leakage current compensation in ferroelectric thin-film capacitor structures," Applied Physics Letters, vol. 86, 2005.
[3] C. R. Harris et al., "Array programming with NumPy," Nature, vol. 585, pp. 357–362, 2020.
[4] J. D. Hunter, "Matplotlib: A 2D Graphics Environment," Computing in Science & Engineering, vol. 9, no. 3, pp. 90–95, 2007.
[5] P. Virtanen et al., "SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python," Nature Methods, vol. 17, pp. 261–272, 2020.
[6] T. Schimansky, "CustomTkinter," 2024. https://github.com/TomSchimansky/CustomTkinter
[7] A. Fratian et al., "piec," GitHub, 2025. https://github.com/ElPsyKurisu/piec
[9] Tian et al., Phys. Rev. Applied 20, 054007 (2023) — HZO Landau coefficients.