The effective Hamiltonian used to study energy storage in some prototypical ferroelectrics (namely, BaTiO₃, PbTiO₃, and KNbO₃) can be expressed as (Nishimatsu et al., 2008; Luo et al., 2016) H eff = M dipole ∗ 2 ∑ i u ˙ i 2 + M acoustic ∗ 2 ∑ i w ˙ i 2 + V s e l f ( { u i } ) + V d p l ( { u i } ) + V s h o r t ( { u i } ) + V e l a s , hom o ( η 1 ,   ⋯ , η 6   ) + V e l a s , i n h o ( { w i } ) + V c o u p , hom o ( { u i } , η 1 ,   ⋯ , η 6 )         + V c o u p , i n h o ( { u i } , { w i } ) − Z ∗ ∑ i E . u i (1) where u i is the local soft-mode amplitude located at the unit cell i, w i is the local acoustic displacement vector, η 1 ,   ⋯ , η 6 represent the component of the homogeneous strain tensor in Voigt notation (Nye, 1985), and M d i p o l e ∗ and M a c o u s t i c ∗ correspond to the effective masses for local soft modes and acoustic displacements, respectively. The different energies appearing in Eq. 1 are the two kinetic energies associated with the local soft modes and acoustic displacements, local mode self-energy, long-range dipole–dipole interaction, short-range interactions associated with the local soft modes, elastic energies stemming from homogeneous and inhomogeneous strains, coupling between local soft modes and both homogeneous in inhomogeneous strains, and effect of an external electric field, E.

The effective Hamiltonian parameters for BaTiO₃, PbTiO₃, and KNbO₃ FE systems were determined from first-principles calculations [see details in References (Waghmare and Rabe, 1997; Waghmare et al., 1998; Nishimatsu et al., 2010)].

An effective Hamiltonian (H _eff) approach was developed in Refs (Xu et al., 2015; Jiang et al., 2021b) for Bi_1−x Nd _x FeO₃ (BNFO) solid solutions. The total internal energy E _int of this effective Hamiltonian can be expressed as E int = E B F O ( { u i } , { η H } , { η I } , { ω i } , { m i } ) + E a l l o y ( { u i } , { ω i } , { m i } , { η l o c } ) (2) where E B F O is the effective Hamiltonian of pure BiFeO₃ (Kornev et al., 2007), while E a l l o y represents the effect of substituting Bi by Nd ions in the A sublattice of this perovskite compound. Technically, the effective Hamiltonian of BNFO comprises four different types of degrees of freedom: i) the local soft mode { u i } centered on the A site i   (Bi or Nd ion), which is directly related to the local electric dipole moment centered on unit cell i (Zhong et al., 1995); ii) strain tensor that includes the contributions of both homogeneous { η H } and inhomogeneous { η I } strains; iii) { ω i } that characterizes the oxygen octahedral tilting about the Fe ion located at the i site of the B sublattice (Kornev et al., 2006); and iv) magnetic moment { m i } centered on the B site of Fe ions (Kornev et al., 2007). Note that a local quantity   η l o c ( i ) = δ R i o n i c 8 ∑ j σ j centered on the B site i of Fe ion, where σ j represents the A sublattice atomic configuration distribution in the BNFO solid solutions and the sum over j runs over the 8 A nearest neighbors of Fe site i ; δ R i o n i c   represents the A sublattice relative difference of the ionic radius between the Nd and Bi ions. Note also that one can add to E _int an energy contribution given by minus the dot product between polarization and applied electric field, in order to simulate the effect of such a field on physical properties (Xu et al., 2017).

BNFO systems are typically modeled by using 12 × 12 × 12 supercells (containing 8,640 atoms) within this H _eff scheme and inside, in which the A site of Bi and Nd ions is randomly distributed. Practically, the total energy of this H _eff approach can be used in Monte Carlo (MC) simulations to compute the static properties of BNFO solid solutions at a finite temperature. A total of 20,000 MC sweeps are typically used for equilibration, and an additional 20,000 MC sweeps are employed to compute statistical thermal averages at a finite temperature and an electric field, to obtain converged results (Xu et al., 2017).

The first-principles-based effective Hamiltonian approach has also been recently developed and used to simulate the static and dynamical properties of bulks and films made of Ba(Zr,Ti)O₃ solid solutions (Akbarzadeh et al., 2012; Prosandeev et al., 2013a; Prosandeev et al., 2013b; Wang et al., 2016; Jiang et al., 2017; Jiang et al., 2022). The total internal energy of the effective Hamiltonian in Ba(Zr,Ti)O₃ contains two main terms too: E int ( { u i } , { v i } , η H , { σ j } ) = E a v e ( { u i } , { v i } , η H ) + E l o c ( { u i } , { v i } , { σ j } ) (3) where { u i } is the local soft mode in unit cell i (which is proportional to the electric dipole moment centered on the B site of Zr or Ti ions), { v i } are the variables related to the inhomogeneous strain inside cell i, η H is the homogeneous strain tensor, and { σ j } describes the B-sublattice distribution (that is, Zr or Ti ion located at the j site) in the Ba(Zr,Ti)O₃ solid solutions. Note that E a v e contains five different contributions: i) the local soft mode self-energy; ii) long-range dipole–dipole interaction; iii) energy as a result of short-range interactions between local soft modes; iv) elastic energy; and v) energy due to the interaction between local soft modes and strains (Zhong et al., 1995). The second energy term E l o c describes how the actual distribution of Zr and Ti ions affects the energetics involving the local soft modes u i and the local strain variables (which depend on the { σ j } distribution).

This effective Hamiltonian has been implemented within MC simulations on 12 × 12 × 12 supercells (8,640 atoms), in order to determine and understand energy storage in disordered Ba(Zr,Ti)O₃ (BZT) relaxor ferroelectrics, for the fixed composition of 50% of both Zr and Ti, in both bulks and epitaxial films (Jiang et al., 2022). Note that this H _eff successfully predicted the existence of three characteristic temperatures in relaxor ferroelectrics. For instance, for BZT bulks with 50% of Zr and Ti ions, 1) the temperature T _m ≅ 130 K, at which the dielectric constant can exhibit a peak (Cross, 1994; Akbarzadeh et al., 2012); 2) the so-called novel critical temperature (T ^∗ ≅ 240 K) at which static polar nanoregions (PNRs) typically occur and that has been recently observed in relaxors (Svitelskiy et al., 2005; Dkhil et al., 2009; Akbarzadeh et al., 2012); and 3) Burns temperature (T _b ≅ 450 K) that marks the dynamical PNRs at which the finite lifetime of polar fluctuations becomes prominent (Burns and Dacol, 1983; Akbarzadeh et al., 2012).

In Ref. (Jiang et al., 2021a), first-principles calculations were performed on (001) epitaxial 1 × 1 AlN/ScN superlattices within the local density approximation to density functional theory (DFT) and using norm-conserving pseudopotentials (Hamann, 2013), as implemented in the ABINIT package (Gonze et al., 2002). The epitaxial strain can be expressed as η i n = ( a − a e q ) / a e q , where a e q corresponds to the in-plane lattice constant of the equilibrium structure of the nonpolar P 6 ¯ m 2 phase (Jiang et al., 2019a). A 6 × 6 × 4 grid of special k-point and plane-wave kinetic energy cutoff of 50 Hartrees were employed. The effects of dc electric fields applied along the pseudocubic [001] direction on structural properties were taken into account by using finite electric-field methods (Nunes and Vanderbilt, 1994; Nunes and Gonze, 2001; Souza et al., 2002; Zwanziger et al., 2012). The electrical polarization P was computed from the Berry phase approach (King-Smith and Vanderbilt, 1993; Resta, 1994). Note that for each considered strain and magnitude of the applied field, the in-plane lattice vectors were kept fixed while the out-of-plane lattice vector and atomic positions were fully relaxed until all the forces acting on the atoms have a value smaller than 10⁻⁶ Hartrees/Bohr, in order to mimic epitaxial films being under electric fields.

Let us first present the energy storage results of the prototypical perovskite ferroelectrics BaTiO₃, PbTiO₃, and KNbO₃. A first-principles-based effective Hamiltonian within molecular dynamics simulations was used for these perovskite ferroelectrics and is described in Section 2.1 and Ref. (Luo et al., 2016).

The temperature dependence of the P-E hysteresis and energy density properties for BaTiO₃, PbTiO₃, and KNbO₃ are shown in Figure 1. The energy density of BaTiO₃ slowly increases with the temperature below the Curie point at 380 K. In contrast, the energy density rapidly increases with a temperature above the Curie point. Furthermore, the energy density behavior of KNbO₃ was numerically found to be very similar to BaTiO₃ because of the similarity in the structural phase transition and field-induced phase transition. In the case of PbTiO₃, the energy density varies nonmonotonically with the temperature, that is, it first reaches a peak around 820 K and then decreases. This can be explained by the fact that PbTiO₃ exhibits a rather weak polarization for temperatures above 820 K. Note that ferroelectrics are not ideal for energy storage applications due to the square shape of their hysteresis loops—which typically gives rise to low energy density and efficiency.

As compared to typical FE systems, antiferroelectric (AFE) materials are very promising for high-power energy storage applications because of their characteristic P-E double hysteresis loops, which is schematized in Figure 2A (Xu et al., 2017). There, E _up describes the critical field at which the AFE-to-FE transition appears upon charging when increasing field, while E _down denotes the critical field for the FE-to-AFE phase transition upon discharging (decreasing field). The green area in Figure 2A represents the energy density.

As indicated in Section 2.2, a first-principles-based effective Hamiltonian method was employed to investigate the energy storage properties in the rare earth-substituted BiFeO₃ multiferroic systems, in general, and disordered BNFO solid solutions, in particular (Xu et al., 2015). This effective Hamiltonian method successfully reproduced the temperature-versus-compositional phase diagram of BNFO (Levin et al., 2010; Levin et al., 2011), showing that the Pnma phase is the equilibrium structure at room temperature for a moderate level of doping Nd. That is the main reason why BNFO was chosen to investigate energy storage properties at room temperature in Ref. (Xu et al., 2017).

Due to the fact that energy density relies on the applied electric field and polarization, the electric fields considered for BNFO solid solutions have been rescaled by a factor of 1/23, so that the computational P-E loop is in good agreement with measurements (Xu et al., 2017). Figure 3 displays the P-E hysteresis curves for four different applied directions of the field (namely, the pseudocubic [001], [100], [110], and [111] directions) for Nd composition changing from 0.4 to 1.0, at room temperature. Note that all the considered compositions adopt the Pnma phase at 300 K under zero field. The polarization increases smoothly within the AFE phase under a small electric field and then abruptly jumps up at the AFE-FE phase transition for fields applied along the pseudocubic [001], [100], [110], and [111] directions, respectively.

Let us now focus on the energy density and efficiency, which can be obtained from the P-E curves from Figure 3. Figures 4A,B show the results obtained for different compositions and electric field orientations, using a maximal applied electric field of E _max = 2.6 MV/cm. The energy density for fields applied along the [001] direction gives the largest values in magnitude. In contrast, the smallest energy densities in Figure 4A correspond to fields applied along the [111] direction. Figure 4B further shows that the highest efficiency corresponds to the electric field applied along the [110] direction.

Three compositions (x = 0.5, 0.7, and 1.0) were considered for the electric field applied along the [001] direction, and one composition of x = 0.5 has been selected for the field along the [110] direction in BNFO solid solutions. The energy storage-related results are shown in Figure 5 [note that the intrinsic breakdown field of E _max = 4.37 MV/cm was considered, which is estimated based on an empirical relation (Wang, 2006) that takes into account the experimental band gap of BiFiO₃ (Ihlefeld et al., 2008)]. More precisely, the related P-E hysteresis curves are shown in Figure 5A, while Figure 5B displays the energy density of BNFO obtained from the P-E hysteresis curves of Figure 5A. The energy densities of other experimentally reported materials are shown for comparison and are smaller than that of BNFO systems. The predicted large energy densities (efficiencies) are 164, 191, and 213 J/cm³ (76%, 88%, and 91%), respectively for x = 0.5, 0.7, and 1, respectively, for the electric field applied along the [001] direction. Similarly, both energy density and efficiency (161 J/cm³ and 91%) of BNFO with x = 0.5 are very large for the field applied along the [110] direction as well. Note also that BNFO superlattices were not numerically found to have different energy storage properties when compared with disordered BNFO systems. These large energy densities and efficiencies therefore indicate that BNFO systems are promising for energy storage.

In order to understand the energy storage results of Figure 4, a simple and general model was proposed for AFE materials, which gives W = P FE 0 E d o w n + 1 2 ϵ 0 χ e E max 2 (4) η = P F E 0 E d o w n + 1 2 ϵ 0 χ e E max 2 P F E 0 E u p + 1 2 ϵ 0 χ e E max 2 (5) where P F E 0 is the polarization of the FE phase along the field’s direction as extrapolated to zero field, and χ e is the dielectric susceptibility for both FE and AFE phases.

Figure 2B shows that the critical field E _C and FE polarization can be expressed in terms of the parameters characterizing the relevant energy landscape of the AFE system, which can be obtained as P F E 0 = ( U F E 0 − U A F E 0 ) / E C (6)

The energy difference U F E 0 − U A F E 0 defines the relative stability of the AFE and FE phases, where the superscript corresponds to E = 0. We also have E d o w n = ( δ E d o w n − δ 0 δ E C − δ 0 ) E C (7) E u p = ( 2 δ E C − δ E d o w n − δ 0 δ E C − δ 0 ) E C (8) where δ 0 , δ E d o w n , and δ E C are the FE-to-AFE barrier at zero field, E d o w n , and E C , respectively. Consequently, combining Eqs 4–8 leads to W = ( δ E d o w n − δ 0 δ E C − δ 0 ) ( U F E 0 − U A F E 0 ) + 1 2 ϵ 0 χ e E max 2 (9) η = ( δ E d o w n − δ 0 δ E C − δ 0 ) ( U F E 0 − U A F E 0 ) + 1 2 ϵ 0 χ e E max 2 ( 2 δ E C − δ E d o w n − δ 0 δ E C − δ 0 ) ( U F E 0 − U A F E 0 ) + 1 2 ϵ 0 χ e E max 2 (10)

Interestingly, the energy density and efficiency data obtained by the effective Hamiltonian MC simulations and shown in Figures 4A,B can be well fitted by Eq. 9 and Eq. 10, respectively (the related parameters are shown in Figures 4C–E). These good fits of Figures 4A,B testify the validity of the aforementioned simple model, which can thus be used to understand and analyze energy storage properties for other AFE materials.

Relaxor ferroelectrics have also attracted special attention because of their high energy densities and efficiencies, which is highly promising for energy storage applications. Despite the fact that the highest recoverable energy density was achieved in Ba(Zr _x Ti_1−x )O₃ relaxor FE thin films at a sustained high electric field of 3.0 MV/cm (Instan et al., 2017), many questions remain to be addressed. For instance, can the first-principles-based effective Hamiltonian method reproduce the experimental finding of ultrahigh energy density in the lead-free relaxor Ba(Zr _x Ti_1−x )O₃ system? What will happen if one considers the epitaxial strain on energy storage properties? In particular, could such ultrahigh energy density be understood by a simple phenomenological model? To the best of our knowledge, the effect of the direction of the applied electric field on energy density and efficiency is also presently unknown in Ba(Zr _x Ti_1−x )O₃ relaxor ferrolectrics, both in their bulk and epitaxial films.

This section is to answer all these questions by using the first-principles-based effective Hamiltonian approach in Ba(Zr_0.5Ti_0.5)O₃ bulk and epitaxial films (Jiang et al., 2022). The effective Hamiltonian is described in Section 2.3.

In order to improve the energy storage performance, it is timely and important to wonder if there are some multifunctional materials awaiting to be discovered/revealed that have 1) ultrahigh energy storage density; 2) optimal 100% energy efficiency; and 3) giant strain levels when under electric fields. Note that a 100% energy efficiency automatically implies that the transition from the paraelectric/AFE to FE states is fully continuous and reversible. To maximize both energy density and efficiency, one can, for example, imagine a nonlinear dielectric material that exhibits large polarization under feasible high electric fields while being initially in a nonpolar phase at a zero field, and with the charging and discharging processes being completely reversible. The III–V semiconductor-based systems made by mixing AlN and ScN to form Al_1−x Sc _x N solid solutions or AlN/ScN superlattices are attractive as these systems can have very large polarization (on the order of 1.0 C/m²) (Jiang et al., 2019b; Fichtner et al., 2019; Noor-A-Alam et al., 2019; Yasuoka et al., 2020; Yazawa et al., 2021), and an FE phase is energetically close to a nonpolar phase (such that the ground state may be tuned by a physical handle) (Jiang et al., 2019b), which is highly promising for the energy storage performance. The aim of Ref (Jiang et al., 2021a) was to demonstrate that ultrahigh energy storage performance can be achieved in 1 × 1 AlN/ScN superlattices, by using first-principles calculations (the method is described in Section 2.4).

Let us now focus on some physical properties of 1 × 1 AlN/ScN superlattices under different epitaxial strains. Previous studies have predicted the existence of different strain-induced regions in 1 × 1 AlN/ScN superlattices (Jiang et al., 2019b), including a wurtzite-derived structure (polar P3m1 space group with a polarization along the c-axis) and a hexagonal-derived phase (paraelectric P 6 ¯ m 2 space group). They are related to each other by a continuous change in the internal parameter (u) and the axial ratio (c/a) (Levin et al., 2010). Figure 12 shows the P-E curves for different compressive strains smaller than −1%, all having the wurtzite-derived structure for ground state (see Figure 13C for such structure). A similar information is displayed in Figure 13A but for compressive strains and tensile strains larger than −0.5%, which result in a hexagonal-derived ground state (see Figure 13B for such phase). The values of the electric fields considered here are the theoretical ones divided by three, in order that the P-E loop of the 1 × 1 AlN/ScN system under a −3% strain becomes very similar to the experimental one corresponding to the Al_0.57Sc_0.43N film (Fichtner et al., 2019)—as demonstrated in Figure 12B. Note that ab-initio electric fields are known to be typically larger than their experimental counterparts (Lu et al., 2019; Chen et al., 2019; Jiang et al., 2020). Note also that the polarizations under zero field computed in the 1 × 1 AlN/ScN superlattice experiencing a compressive strain of −3% and the ones measured in the Al_0.57Sc_0.43N system were found to be very close to each other (see Figure 12A), which explains why we decided to determine the renormalization factor of the theoretical fields by comparing the P-E loops of these two compounds. Regarding the magnitude of the (renormalized) fields to be considered here, it is first worthwhile to indicate that a recent experiment reported a feasible electric field as high as ∼5 MV/cm in Al_1−x Sc _x N films (Fichtner et al., 2019). Moreover, one can also estimate the intrinsic breakdown field of such systems via the universal expression (Wang, 2006), that is E B I = 1.36 × 10 9 ( E g 4.0 ) α (V/m), where E _g is the band gap, and α = 1 for insulators, and α = 3 for semiconductors. In such expression, a theoretical corrected band gap of ∼3.1 eV is considered for the 1 × 1 AlN/ScN superlattice ground state (Jiang et al., 2019b). An intrinsic breakdown field electric field, E _break, of about 6.3 MV/cm is obtained, which is the maximal value of applied electric fields that will be considered below.

Several features are remarkable in Figure 12B and Figure 13A. First of all, the P-E loops of Figure 12B are those typical of ferroelectrics with first-order transitions (between two wurtzite-derived structures) for some values of polarization changing of sign under the application of initially opposite electric fields, but with the magnitude of the critical fields of such transitions being very strongly strain dependent. For example, the magnitude of these critical fields is 5.8 MV/cm versus 0.2 MV/cm for strains of −5% and −1%, respectively. Such a strong dependency of these fields has also been observed but with respect to the Sc composition in Al_1−x Sc _x N solid solutions (Fichtner et al., 2019). Consequently, varying composition in these solid solutions or changing the strain for a fixed alloy or superlattice should yield similar physics and results. Secondly and strikingly, Figure 13A reveals a nonlinear behavior of the polarization for small fields (starting from the paraelectric P 6 ¯ m 2 hexagonal-derived phase at a zero field), with the transition between the paraelectric hexagonal-like structure to an FE wurtzite-type phase being continuous and fully reversible upon increasing and decreasing electric fields. Note that such type of continuous and reversible transition has been mentioned in relaxor ferroelectrics (Jiang et al., 2022). As we are going to see, such latter features are indeed promising for applications in electronics and electric power systems. Moreover, larger strain makes the polarization less nonlinear at small fields and also results in a smaller out-of-plane polarization at large fields (Figure 13A), which are behaviors that can affect energy storage density (as we also discuss below).

Let us now focus on energy storage properties for the 1 × 1 AlN/ScN superlattices for misfit strains changing from −0.5% to +1%. Figure 13A shows that the energy efficiency is 100% because the charging and discharging processes are completely reversible. Regarding the calculated energy density, Figure 13D reports it as a function of epitaxial strain ranging between −0.5% and +1%, for electric fields up to two different values, namely, E _max = 5 MV/cm, which has been experimentally realized (Fichtner et al., 2019), and E _break = 6.3 MV/cm which has been estimated as the intrinsic breakdown field. The energy density can reach very large values, that is, it varies from 127 to 135 J/cm³ and from 187 to 200 J/cm³ when using 5 and 6.3 MV/cm for the maximal applied field, respectively. All these large energy densities with the ideal efficiency of 100% thus indicate that the 1 × 1 AlN/ScN superlattices are highly promising for high-power energy storage applications.

To understand these energy storage results of Figure 13D, we can also use Eqs 11–13 to describe the behaviors in initially nonpolar 1 × 1 AlN/ScN superlattices. Figure 14A displays the E-P data for electric field applied along the pseudocubic [001] direction at a 0.25% strain. The yellow area represents the energy density. The E-P data can be well fitted by Eq. 12 for all the considered strains ranging from −0.5 to +1%.

Figure 14B displays the a and b fitting parameters versus strain for fields up to the intrinsic breakdown field E _break = 6.3 MV/cm. It is worth noting that the a parameter linearly increases with strain ranging from −0.5 to +1%. In contrast, the b parameter linearly decreases with strain. Note that the change in a is much larger in percentage than the change in b. Figure 14C shows the maximum polarization P _max as a function of strain obtained from DFT calculations and the Landau model at E _max = 6.3 MV/cm. For all considered strains, the DFT and Landau model provide nearly identical results. In order to understand the results of Figure 13D and Eq. 13 is used to obtain the energy density at E _max = 6.3 MV/cm for strain ranging between −0.5 and +1%. Figure 14D shows that the Landau model and DFT-obtained energy density behave in a similar qualitative and even quantitative way with the strain, implying that the Landau model can be trusted.

Once again, Eq. 13 indicates that the energy density can be decomposed into two terms. The first term relies on the product of the fitting a parameter and P max 2 , while the second term depends on the product between the b parameter and P max 4 . Figures 15A,B display the total and decomposed energy densities versus strain obtained from Eq. 13 for two different maximum fields E _max = 5.0 and 6.3 MV/cm, respectively. Figure 15 shows that the Landau model predicted total energy density exhibits a maximum at 0 and 0.25% misfit strain when the largest applied fields are 5.0 and 6.3 MV/cm, respectively. Such numerical findings agree very well with the DFT-obtained results of Figure 13D. The strain dependency of 1 2 a P max 2 and 1 4 b P max 4 eventually give rise to a maximum value of the total energy density occurring at different specific strain for each considered maximal electric field.