Nonvolatile multistate memories in BiFeO₃-based nanocomposites designed by phase-field simulations

INTRODUCTION

Ferroelectric (FE) materials, characterized by their switchable spontaneous polarization, are widely used in information storage applications such as ferroelectric random access memory (FeRAM)^[1,2]. Compared to other nonvolatile memory technologies, FE-based devices offer several advantages, including fast read/write operations, low power consumption, and exceptional durability^[3,4]. However, traditional FeRAM technology is constrained by limited data storage capacity, hindering their further integration in advanced computing systems^[5]. To address this challenge, one promising strategy involves increasing the number of saturated polarization states within a single unit, thereby enabling multiple memory states per cell beyond the conventional binary states of “0” and “1”^[6]. This strategy not only enhances storage density but also holds promise for neuromorphic computing applications in today’s era of big data and artificial intelligence^[4,7].

Considering various domains in FE materials, one approach to achieving multiple polarization states is to control the stepwise domain switching pathways under external electric fields^[8-12]. For example, intermediate polarization values have been achieved by modulating up/down polarization ratios in PbZr_1-xTi_xO₃ (PZT) thin film^[9]. Similarly, multiple thermodynamically stable twinning structures can form in (111)-oriented PZT films during the ferroelastic switching process^[11]. While these methods can produce abundant states under varying electric fields, precise control over these polarization values remains challenging. Another approach focuses on creating regions with distinct polarization behaviors under electric fields by modifying material composition, introducing defects, or adjusting external conditions^[13-15]. For instance, Pb_1-xSr_xTiO₃ exhibits varying coercive fields based on strontium contents; hence, a bilayer with different compositions results in an in-plane (IP) four-state switching accompanied by interlayer coupling^[14]. However, this approach may compromise long-term device reliability due to the potential degradation of interfacial couplings. The third approach leverages inherent multiple polarization states within single FE materials^[16,17]. Research on (110)-oriented BiFeO₃ (BFO) films has demonstrated a three-state switching^[16]. While this method promises improved stability and scalability, such a system remains limited for practical computational usage due to its incompatibility with the binary architecture.

In fact, many FE materials exhibit multiple symmetrical phases, such as PZT near its morphotropic phase boundary (MPB)^[18,19]. Nevertheless, stabilizing individual phases via electric fields is challenging, and the work to utilize inherent multiple phases for multistate memory applications remains largely unexplored.

BFO, a multiferroic material displaying both ferroelectricity and antiferromagnetism at room temperature^[20,21], is a prominent candidate for next-generation nanoelectronics. Its notable properties, such as large remnant polarization and strong coupling between ferroic orders, make it an ideal platform for exploring multistate memory applications^[22]. In addition to a rhombohedral (R) phase, BFO also exhibits a super-tetragonal (T^S) phase characterized by a large polarization value and a high c/a ratio^[23]. Figure 1A illustrates the free energy landscape of BFO. The color on the spherical surface indicates the minimal energy density along the corresponding directional angle. As shown in the left panel of Figure 1A, under stress-free conditions, both R and T^S phases are present. However, the T^S phase exhibits significantly higher energy than that of the R phase, rendering the T^S phase unstable and prone to transitioning to the R phase. It should also be noted that the T^S phase does not possess tetragonal symmetry but rather monoclinic symmetry^[24,25], while its out-of-plane (OOP) polarization is substantially larger (~ 1.4 C/m²) compared to the R phase (~ 0.6 C/m²)^[25]. When subjected to large IP compressive strain and/or OOP tensile stress, the T^S phase becomes energetically favorable, whereas the R phase is constrained to the distorted rhombohedral (R^D) phase and becomes metastable. The right panel of Figure 1A shows the case under IP compressive strain of -1.3% and OOP tensile stress of 9.12 GPa (equals to a tensile mismatch of 4%), and it is similar to the one under IP strain of -4%. In this scenario, strain-induced phase separation allows the coexistence of both R^D and T^S phases^[23,26,27], presenting an opportunity to stabilize four distinct polarization states.

Figure 1. Design of multistate BFO nanocomposites. (A) Polar plots of free energies for BFO under stress-free conditions (left) and under an IP compressive strain of -1.3% combined with an OOP tensile stress of 9.12 GPa (right). The right panel is similar to the case of a high in-plane compressive strain of -4%, demonstrating comparable stabilization of the super-tetragonal phase; (C) Schematic of the proposed nanocomposite structure consisting of a substrate layer, a film layer of a BFO pillar surrounded by a dielectric matrix, and an air layer; (D) Designed energy density profiles of a plain BFO film and a BFO pillar embedded in a dielectric matrix; (E) Electric hysteresis loop of the designed BFO nanocomposite film. BFO: BiFeO₃; IP: in-plane; OOP: out-of-plane.

Despite its potential, plain BFO films under large IP compressive strain cannot maintain a single R phase due to intrinsic constraints. Figure 1B presents the polarization-electric field (P-E) hysteresis loop of plain BFO films subjected to an IP strain of -4%. The resulting square-shaped loop highlights the inability to stabilize the R phase under these conditions, which impedes the direct utilization of high-strained BFO films for achieving multiple polarization states.

In recent years, self-assembled ferroelectric/dielectric nanocomposites have attracted considerable attention for their potential in advanced nanotechnologies^[28]. These nanocomposites provide enhanced control over the functions of FE phases, offering a versatile platform for the design of multistate memory systems^[29]. Specifically, vertically aligned nanocomposites have demonstrated notable improvements in the ferroelectric, piezoelectric, and dielectric properties of functional materials^[30,31], while also enhancing overall stability. These characteristics make them promising candidates for stabilizing both the R^D and T^S phases in BFO films, thereby enabling the realization of multistate memory devices.

In this study, we propose a novel structure in which a BFO circular pillar is embedded within a dielectric matrix [Figure 1C]. The tensile stress arising from the larger lattice constants of the matrix promotes the stabilization of the T^S phase^[31]. In addition, careful interface engineering may lead to similar energy minima for both R^D and T^S phases [Figure 1D], facilitating the stabilization of each single phase. Using phase-field simulations, we successfully stabilize four distinct states in this configuration [Figure 1E], thereby demonstrating its potential for high-density memory applications.

METHODS

In this study, we employed a phase-field model to explore the polarization dynamics and multistate design of BFO/dielectric nanocomposites. The polarization vector P_i (i = 1, 2, 3) was chosen as the order parameter. The total free energy can be written as^[32]:

where f_Land, f_grad, f_elas, and f_elec represent the volume densities of Landau free energy, gradient energy, elastic energy, and electrostatic energy, respectively.

The Landau free energy density is extended to eighth order to describe both the R and T^S phases:

where α₁, α₁₁, α₁₂, α₁₁₁, α₁₁₂, α₁₂₃, α₁₁₁₁, α₁₁₁₂, α₁₁₂₂, and α₁₁₂₃ are the Landau polynomial coefficients.

The gradient energy density, arising from the discontinuity of order parameters and accounting for the contribution of interfaces, is given by:

where G₁₁, G₁₂, G₄₄, and G₄₄ are the gradient energy density coefficients, and P_i,j denotes ∂P_i / ∂x_j with x_j as the spatial coordinate.

The elastic energy density, describing the interactions between strain and stress, is expressed as:

where c_ijkl is the elastic stiffness tensor, ε_ij is the total strain, and ε⁰_ij = Q_ijklP_kP_l is the eigenstrain with Q_ijkl as the electrostrictive coefficients. The Einstein notation is used.

The electrostatic energy density is:

where E_i = - ∂φ / ∂x_i is the electric field with φ as the electrostatic potential, ε₀ is the permittivity of vacuum, and κ_b is the background relative permittivity of materials.

The time evolution of the order parameters is governed by the time-dependent Ginzburg-Landau (TDGL) equations:

where L_P is the kinetic coefficient. The TDGL equations are solved using the semi-implicit Fourier spectral method^[33]. For each step, the elastic and electric driving forces are determined by solving the mechanical equilibrium equations σ_ij,j = 0 and the electrostatic equilibrium equation D_i,i = 0, where σ_ij = c_ijkl (ε_ij - ε⁰_ij) is the local stress and D_i = ε₀κ_bE_i + P_i is the electric displacement. The superposition method is employed to solve film boundary conditions^[34,35].

The mechanical boundary conditions are mixed, i.e., ε₁₁ = ε₂₂ = ε_s, ε₁₂ = ε₂₁ = 0, and σ₃₃ = σ₁₃ = σ₃₁ = σ₂₃ = σ₃₂ = 0, where ε_s is the IP strain imposed by substrate^[34]. The electric boundary conditions are closed-circuit^[36]. The potential of the bottom film-substate interface is set to zero. Under a uniform electric field, the potential on the top surface is also uniform. If applying a piezoresponse force microscopy (PFM) tip, the surface potential distribution follows a two-dimensional Lorentz-like function^[37,38]:

where φ₀ is the applied potential, i.e., the maximal potential,γ is the length of half-maximum potential, and (x₁⁰,x₂⁰) is the tip location.

In our simulations, we assumed that the substrate material is SrTiO₃ (STO) and the surrounding dielectric material is MgO. Based on the lattice constants of these materials (BFO: 3.96 Å, STO: 3.91 Å, and MgO: 4.21 Å) and setting BFO as the reference, the IP strain is -1.3%, and the mismatch strain of MgO is 4% (lower than the theoretical 6.3% based on lattice constants) to account for dislocation-induced strain relaxation at BFO-MgO interfaces^[39]. For BFO, according to previous literature^[27], α₁ = 4.184 × 10⁵ × (T - 1123.5) C^-2m²N with T = 298 (in Kelvin) as the ambient temperature, α₁₁ = 2.127 × 10⁹ C^-4m⁶N, α₁₂ = -2.049 × 10⁹ C^-4m⁶N, α₁₁₁ = -1.76 × 10⁹ C^-6m¹⁰N, α₁₁₂ = 8.298 × 10⁸ C^-6m¹⁰N, α₁₂₃ = 3.358 × 10⁹ C^-6m¹⁰N, α₁₁₁₁ = 3.92 × 10⁸ C^-8m¹⁴N, α₁₁₁₂ = 4.400 × 10⁷ C^-8m¹⁴N, α₁₁₂₂ = -3.8 × 10⁸ C^-8m¹⁴N, α₁₁₂₃ = 8 × 10⁸ C^-8m¹⁴N, c₁₁ = 228 GPa, c₁₂ = 128 GPa, c₄₄ = 65 GPa, and κ_b = 40. For MgO, ₁ = 5.0 × 10⁹ C²m²N, α₁₁ = α₁₂ = α₁₁₁ = α₁₁₂ = α₁₂₃ = α₁₁₁₁ = α₁₁₁₂ = α₁₁₂₂ = α₁₁₂₃ = 0, c₁₁ = 293 GPa, c₁₂ = 92 GPa, c₄₄ = 155 GPa, and κ_b = 10. For STO, c₁₁ = 228 GPa, c₁₂ = 128 GPa, and c₄₄ = 65 GPa. Voigt notation is used. Additionally, Q₁₁ = 0.071 C^-2m⁴, Q₁₂ = -0.030 C^-2m⁴, Q₄₄ = -0.295 C^-2m⁴, G₁₁ = 1.2G₁₁₀, G₁₂ = 0, and G₄₄ = G’₄₄ = 0.6G₁₁₀, where G₁₁₀ = 3.45410^-10 C^-2m⁴N. The spectral iterative perturbation method is used^[40]. For a PFM tip, the location (x₁⁰, x₂⁰) is at the center of the BFO pillar, and γ = 3.0 nm.

The size of one unit case is 64∆x × 64∆x × 40∆x with ∆x = 0.4 nm, including 15 layers for the substrate and 15 layers for the BFO-MgO film, as illustrated in Figure 1C. The size of the 4 × 4 cell array is 200∆x × 200∆x × 40∆x, as shown below.

RESULTS AND DISCUSSION

Stability of four states

Using the phase-field model, we first examined the stability of the four potential polarization states by varying the diameter of the BFO pillars. Figure 2A shows the P-E loops for BFO-MgO systems with different BFO pillar diameters (d) under an electric field of 2 MV/cm. For larger pillar diameters (d = 12 nm), the results resemble those of plain BFO films [Figure 1B], albeit with a significant reduction in the coercive field due to the influence of the dielectric interfaces.

Figure 2. Influence of pillar size on phase stability. (A) Electric hysteresis loops of BFO pillar films at 2 MV/cm for various pillar diameters d; (B) Three-dimensional polarization vectors corresponding to four stable states for d = 8 nm. The vector colors represent the magnitude of out-of-plane polarization P_z; (C and D) Free energy densities and polarization evolution for d = 8 nm; (C) and for d = 6 nm and 12 nm (D) from specific states in the hysteresis loops indicated by purple triangles in (A) under zero electric field. BFO: BiFeO₃； OOP: out-of-plane.

In plain strained BFO films, the T^S phase typically exhibits much lower elastic energy compared to the R^D phase, despite having higher Landau free energy. However, in BFO-MgO nanocomposites, the low permittivity of surrounding dielectrics leads to reduced polarization even under high fields. This causes substantial gradient and depolarization interactions in the vicinity of the BFO/MgO interfaces. To minimize these energies, the BFO dipoles near interfaces tend to decrease their magnitudes and align parallel to the interfaces. Similar behaviors have been observed in PbTiO₃/SrTiO₃ superlattices near FE/dielectric interface^[41]. As a result, the T^S phases with large polarization are less favored near the interfaces, thereby lowering the energy barrier between up and down polarizations of the T^S phases. Conversely, the R^D phases, which exhibit relatively smaller polarization, are more stable near the interfaces. Therefore, at d = 12 nm, intermediate R^D phases appear during the polarization switching, although this effect is subtle.

As the pillar diameter decreases, the surface-to-volume ratio increases, enhancing the stability of the R^D phase. When the diameter was reduced to 8 nm, both the R^D and T^S phases were found to be stable. At even smaller diameters, the T^S phase becomes metastable; it appears only when an electric field is applied. For pillar diameters smaller than 4 nm, the T^S phase could no longer be stabilized under the applied electric field.

Figure 2B illustrates the four distinct polarization states at d = 8 nm, with these states being referred to as +H, +L, -H, and -L for simplicity. They correspond to the T^S phase (+H and -H) and the R^D phase (+L and -L). By extracting two specific points from the P-E loop at d = 8 nm (indicated by purple triangles) and allowing the system to relax after removing the external electric field, we confirmed the stability of both phases over time [Figure 2C]. From an energy perspective, the total energy density of the two phases is comparable, consistent with our designed energy landscape shown in Figure 1D.

Further analysis revealed that the diameter of the BFO pillars significantly influences phase stability. For d = 6 nm, the T^S phase under zero electric field is unstable due to large gradient and Landau energy contributions [Figure 2D, left panel]. In contrast, for d = 12 nm, the initial R^D phase quickly transitions to the T^S phase due to the dominance of elastic energy [Figure 2D, right panel]. These findings emphasize the critical role of BFO/MgO interfaces, governed by pillar size, in balancing the free energy contributions and stabilizing desired phases. For subsequent investigations, we selected d = 8 nm as the optimal pillar diameter.

Switching dynamics among four states

Having confirmed the stability of the four polarization states, we next examined the switching dynamics among these states under pulsed electric fields. Here, we use +E^H, +E^L, -E^L, and -E^H to represent the fields of +1.8, +1, -1, and -1.8 MV/cm, respectively. Figure 3A shows the polarization evolution as a function of the pulsed electric field, indicating that all four states can be easily switched and stabilized by applying external electric fields of diferent directions and magnitudes. Figure 3B summarizes all the possible switching routes between every two states.

Figure 3. Switching dynamics between multiple states. (A) Polarization evolution as a function of external pulsed electric field; (B) Schematical illustration of all possible switching routes between the four stable states controlled by electric fields; (C) Two proposed approaches to obtaining all states from any initial state.

To assess the potential for localized control, we also conducted simulations using a PFM tip. We found that applying voltages of ± 1.2V and ± 0.6V corresponded to the electric fields of $$ {\overset{-}{+}} $$E^H and $$ {\overset{-}{+}} $$E^L, respectively. These results validate the feasibility of precise and flexible control over the polarization states.

For practical implementation, we propose two approaches to achieve transitions from any starting state to all four states, as shown in Figure 3C. For the first approach, an external +E^H field is applied to initialize the system in the +H state. The -H and +L states can then be obtained by applying the -E^H and -E^L fields, respectively. To reach the -L state, an additional +E^L field is applied to the -H state. The second approach involves applying +E^H and -E^H fields to establish the +H and -H states, respectively. The +L and -L states can then be achieved by applying the -E^L and +E^L fields, respectively.

To quantify the energy barrier between states, we then analyzed the switching fields and found that E^L ranges from 0.7 to 1.4 MV/cm, while E^H should be at least 1.5 MV/cm. Considering the polarization values of H and L states(~ 1.0 and 0.5 C/m²), the energy density barrier between two adjacent high/low states (e.g., +H to +L) is estimated as P · E^L ≈ 3.5 × 10⁷ J/m³, while transitions between opposite states (e.g., +H to -H) require P · E^H ≈ 3.0 × 10⁸ J/m³. Assuming a cell volume of ~ 300 nm³ (IP diameter d = 8 nm and film thickness h = 6 nm), the estimated switching energy per operation is ~ 10^-5 to 10^-4 pJ, suggesting ultralow power consumption compared to conventional FeRAM. This is primarily due to the nanoscale size of memory cells, though practical implementations would require further integration with transistor-based architecture.

A demonstration of cell array

To demonstrate the potential for high-density data storage, we also simulated a 4 × 4 memory cell array using the proposed multistate memory design. The layout of the proposed array is shown in Figure 4A. Each cell has a diameter of 8 nm, with a spacing of 12 nm between adjacent cells, sufficiently avoiding the interactions between neighboring cells. To switch a specific cell, a PFM tip is applied at the center of the targeted cell with an appropriate voltage.

Figure 4. Demonstration of data writing in a 4 × 4 cell array. (A) Simulated structure of a 4 × 4 cell array, each cell 20 nm × 20 nm with BFO pillar diameter of 8 nm. The switching process is performed using a PFM tip; (B) Design of data storage, where each 8-bit ASCII character is stored by four cells; (C) Phase-field simulation of the data writing process in (B), beginning from an initial poled state with all cells in the “00” state, and then writing each cell sequentially by one or two steps. BFO: BiFeO₃; PFM: piezoresponse force microscopy.

By encoding the four polarization states as binary bits (“00”, “01”, “10”, and “11” correspond to +H, +L, -L, and -H states, respectively), the storage density is effectively doubled. This allows for more data to be stored within the same physical area [Figure 4B]. For example, knowing that storing one ASCII character requires 8 bits, four ASCII characters can be stored in 16 cells with the proposed design, instead of 32 cells by traditional designs.

Figure 4C shows the successful switching process of individual cells in the 4 × 4 array to store desired information in Figure 4B. We first applied a uniform +E^H field, initializing all cells in the “00” states. Individual cells are then switched to the desired states sequentially using the PFM tip, following a switching process of zero, one, or two steps, as outlined in the left panel of Figure 3C. For example, the first binary-bit cell was applied by a voltage of 0.6 V, switching to the “01” state. The second cell was already in the desired “00” state and did not require additional manipulation. The third cell underwent a sequential voltage application of 1.2 V, followed by -0.6 V, ultimately reaching the “11” state. After a series of steps, the designed information configuration was successfully achieved. This array configuration (pillar diameter d = 8 nm and center-to-center distance of adjacent pillar is 20 nm) results in each cell occupying 400 nm². With 2-bit/cell encoding, the theoretical storage density is 2 bit / 400 nm² = 500 Gbit/cm². This value remarkably surpasses the storage capacity of traditional FeRAM devices, though practical limitations such as read/write circuitry and thermal stability may reduce the achievable densities.

CONCLUSIONS

In conclusion, we have proposed an innovative design for nonvolatile multistate memories based on BFO-dielectric nanocomposites. Using phase-field simulations, we successfully demonstrated the stabilization and switching of four distinct polarization states. By carefully optimizing the size of the BFO pillars, both the R^D and T^S phases could be stabilized through the subtle balance of free energies, especially the gradient and elastic energy terms. These polarization states can be effectively switched by applying external electric fields or localized voltages via a PFM tip. Low energy consumption and large polarization contrast ensure low-cost and reliable read/write operations. Additionally, our simulation of a 4 × 4 memory cell array illustrates the potential for high-density data storage, offering substantial advantages over traditional FeRAM devices. However, interfacial fatigue during cycling, multistep write operations, and the trade-off between applied voltage and leakage in nanocomposite films remain critical challenges, requiring materials engineering and interfacial optimization. This work not only provides a novel approach for boosting the capacity of FE-based memory devices but also opens up new possibilities for designing next-generation electronic nanotechnologies. Future research could extend this approach to other materials, such as other perovskite ferroelectrics like PbTiO₃(exhibiting T^S phases under strain), solid solutions like PZT near MPBs (enabling rhombohedral and tetragonal phase coexistence), and van der Waals materials like CuInP₂S₆ (with quadruple-well polarization potentials^[15,42,43]), to further optimize data capabilities and switching flexibility.