Introduction

Multilayer ceramic capacitors (MLCCs) are extensively utilized in electronic power systems owing to their remarkable charge/discharge capabilities and good stability across a wide temperature range1,2,3. However, their low energy storage density Wrec and/or low efficiency η have seriously hindered their development in the miniaturization and integration of electronic devices4,5. In general, Wrec and η can be calculated by the following equations: Wrec = \({\int }_{{P}_{r}}^{{P}_{m}}{EdP}\) and η = Wrec/(Wrec+Wloss), where E, Pr, and Pm denote the electric field, remnant polarization, and maximum polarization, respectively (Supplementary Fig. S1). Therefore, the simultaneous realization of a small Pr, a large Pm, and a high E is essential for obtaining ultra-high Wrec and η.

Recently, different types of lead-free dielectric materials have been used to prepare MLCCs to improve the energy storage performance (ESP), including antiferroelectrics (AFEs)6,7,8, relaxor ferroelectrics (RFEs)9,10,11, normal ferroelectrics (FEs)12,13,14, and linear dielectric materials15. Among these dielectrics, AFEs dielectrics, characterized by double hysteresis loops, have garnered considerable attention for energy storage applications due to relatively high Pm and low Pr compared to other dielectrics. NaNbO3-based ceramic materials, as representatives of the lead-free antiferroelectric system, show very great potential for energy storage due to their wide bandgap (~3.45 eV), high polarization strength (~40 μC▪cm−2) and small bulk density (~4.55 g▪cm−3)16. Zhang et al.17 improved the energy storage efficiency from 30% to 90% in NN-based ceramic materials with tailored functionality from antiferroelectric to relaxation states through local structural modifications and changes in defect chemistry. However, the energy storage density is low at 1.7 J▪cm−3. Lv et al.7 obtained high Wrec = 12.56 J▪cm−3 in Bi(Mg0.5Hf0.5)O3 doped NaNbO3-Bi0.5Na0.5TiO3 MLCCs system by minimizing the ___domain size to generate polar nanoregions (PNRs) and refining the grain structure to boost the insulating performance. However, the lower energy storage efficiency (<90%), in which more than 10% of the energy is lost in heat, hampers their application in electronic components. Therefore, simultaneously obtaining high Wrec (>14 J▪cm−3) and ultra-high η (>95%) remains a serious challenge in NN-based ceramics/MLCCs. In addition, it is widely recognized that the formation of PNRs in NN-based antiferroelectric material systems is a reason for obtaining excellent ESP due to the breakdown of long-range ordered antiferroelectrics. However, people seem to have overlooked the question: Is antiferroelectric NN a nonpolar matrix? The answer is negative. Therefore, it is critical to gain a deeper understanding of ESP enhancement in NN-based systems.

In this work, we strategically design a method to achieve superior ESP using phase field simulations. Due to the doping of guest phase Bi(Ni0.5Hf0.5)O3, the long-range ordered NaNbO3 antiferroelectric microdomains are broken into nanodomains embedded in a polar orthorhombic matrix with a certain orientation, which is called directional slush-like polar (DSP) structure design for simplicity. As a result, a gigantic Wrec = 14.1 J▪cm−3 accompanied by an ultra-high η = 96.8% at a critical breakdown electric field of 1215 kV▪cm−1 was realized in NN-0.15BNH MLCCs. Furthermore, the relationship between structural evolution and ESP was thoroughly investigated.

Results and discussion

Phase-field simulations of ___domain structure evolution

Phase-field simulation, a powerful computational tool, is essential for characterizing the evolution of microstructures and their associated physical properties, particularly in energy storage dielectric materials, where it is applied to simulate phenomena such as solidification, crack propagation, and breakdown18,19. Here, we use the NaNbO3 (NN) system as the main phase and Bi(Ni0.5Hf0.5)O3 (BNH) as the guest phase for constructing a phase-field simulation. It is widely acknowledged that NN undergoes an irreversible antiferroelectric-ferroelectric phase transition (AFE-FE) at room temperature, resulting in the emergence of ferroelectric behavior upon repeated application of an electric field17. The computational parameters were determined based on the NN-0.05BNH system since the doping of a small amount of BNH (0.05 mol) stabilizes the AFE orthorhombic (Pbcm) phase (this will be confirmed in the next section).

The phase field simulation results as a function of BNH guest phase content (Fig. 1a) indicate that the NN-0.05BNH system is an O phase, confirmed by previous results20. We use vectors of different colors to express different polarization directions. However, distinguishing AFE domains from other ___domain structures is very difficult. Therefore, we further use blue, yellow, and gray as background colors to assist in identifying different ___domain structures. It is worth noting that blue represents AFE domains, yellow reflects a state in which long-range ordered AFE domains are broken by a guest phase, which we call non-antiferroelectric (NAFE) domains, and gray represents the paraelectric phase (Supplementary Fig. S2). It is clear that the long-range ordered AFE domains are disrupted by nanodomains and that the NAFE ___domain content gradually increases, allowing the material to have pronounced relaxation properties when x increases. The effect of external electric field and temperature on the polarization direction for different guest phase doping was calculated and displayed in Fig. 1a and Supplementary Fig. S2. The x = 0.05 system possesses conventional AFE states accompanied by large AFE domains. The direction of polarization is nearly aligned with the direction of the electric field, and the initial state cannot be completely recovered once the electric field is removed, implying that some degree of polarization distribution persists, which leads to higher remanent polarization (Pr) (Fig. 1b). As the doping content increases, the polarization direction becomes more disordered, the domains exhibit a stronger response to the external electric field, and the polarization distribution can quickly return to its initial state once the electric field is removed. These microstructural changes make the P-E loop smaller (Fig. 1b), which facilitates the ESP. The reason we chose x = 0.15 as the target component is that it has a reasonable polarization ordered-disordered arrangement. Calculated P-E loops and microstructural evolutions related to the x = 0.15 component (Fig. 1c) indicate that the polarization structure evolution (O1-P1-N-P2-O2) is very sensitive to the external electric field, while the polarization structures at O1 and O2 are essentially the same, implying a very high activity, leading to a low Pr.

Fig. 1: The phase-field simulation results of (1-x)NN-xBNH.
figure 1

a Calculated microstructure evolutions under given fields with vector contours. b Calculated P-E loop of (1-x)NN-xBNH ceramics. c Calculated P-E hysteresis loop and microstructure evolution of x = 0.15 ceramic at room temperature. The blue color represents the antiferroelectric ___domain, and yellow color denotes the non-antiferroelectric ___domain.

Full size image

Energy storage capability, stability, and charge/discharge properties

Based on the results of the phase field simulations, x = 0.15 was chosen as the target component for further investigation. NN-based MLCCs of the x = 0.15 component with six dielectric layers and each layer thickness of approximately 6 μm were successfully prepared by a tape casting process. We realized a dense dielectric layer with a continuous electrode layer, while there was no obvious elemental diffusion in the dielectric and electrode layers, as demonstrated in Fig. 2a and Supplementary Fig. S3. The impact of breakdown field strength Eb on the ESP of MLCCs is critical. Figure 2b illustrates the room temperature P-E loop for x = 0.15 (MLCC) tested at 700–1215 kV▪cm−1. It is noteworthy that Pr continues to remain at a low level as the electric field increases, even when the electric field reaches 1250 kV▪cm−1, which implies ultra-high η. The Weibull distribution experiments represent a crucial instrument for assessing the stability of the breakdown field strength5,21. The Weibull modulus β value of 34.25 was obtained by linear fitting, reflecting a shorter distribution while improving the stability of the Eb data (Fig. 2c). As a result, the giant Wrec = 14.1 J▪cm−3 and ultra-high η = 96.8% were simultaneously obtained for the x = 0.15 (MLCC) component at the critical breakdown electric field (Fig. 2d). We use the figure of merit UF = Wrec/(1-η) to reflect the overall ESP, reached a high value of 440.6 owing to the combined high Wrec and ultrahigh η, which outperform the major of recently reported state-of-the-art MLCCs1,2,6,7,9,10,12,13,22,23,24,25,26, see Fig. 2e.

Fig. 2: Energy-storage performance of NN-based MLCC devices.
figure 2

a Digital image of the MLCCs and SEM image of the cross-section area with corresponding element distribution. b The P-E loop of the MLCCs. c The Weibull distribution analysis of the breakdown strength of the MLCCs. d Energy density and efficiency as a function of electric field up to Eb. e Comparations of UF (unit. J▪cm−3) versus η of our NN-based MLCCs with state-of-the-art dielectric MLCCs. f Temperature stability of the energy-storage performance at 700 kV▪cm−1. g Time dependence of discharge energy density under different fields (R = 1000 Ω).

Full size image

Temperature-induced, frequency-induced, and cycling-induced stabilization of the ESP must also be considered to ensure stable operation in complex environments. We tested the thermal stability of the ESP concerning x = 0.15 (MLCC) samples in the range of from 20 °C to 150 °C (Fig. 2f and Supplementary Fig. S4) and found that it demonstrated good performance stability at 720 kV▪cm−1, possessing high Wrec (6.11 ± 0.25 J▪cm−3) and η (93.15 ± 6.15%), which is attributed to the polarization distribution state and ___domain size being unaffected by temperature changes (Supplementary Fig. S2). Frequency-induced ESP is also evaluated (Supplementary Fig. S5a, c), demonstrating that the samples exhibit excellent frequency stability possessing a deterioration of less than 4% for Wrec and η over the entire test frequency range (1–100 Hz). The test frequency affects the response of the material. The hysteresis loss (Wloss) with change in frequency is provided in Supplementary Fig. S6. The nanodomains flip with the change in electric field, and the flip is accompanied by ___domain wall motion. As the test frequency increases, the flip period of the electric field decreases, and a part of the nanodomains cannot flip, showing a decrease in the residual polarization, leading to a reduction in the Wloss. In addition, the ___domain wall motion during nanodomain flip is resisted by viscous-like forces between the nanodomains, and this resistance increases with the frequency, making a larger electric field required for a complete flip of the nanodomains, exhibiting an increase in the coercive field and leading to a further rise in Wloss. During accelerated charging and discharging at 1 kHz, it can be seen that x = 0.15 (MLCC) exhibits a small performance deterioration (ΔWrec ~ 2%, Δη ~ 1%) after performing 1 × 105 cycles (Supplementary Fig. S5b, d). The elimination of macroscopic ___domain walls and the high dynamics of nanodomains suppress the ___domain wall pinning effect, which is beneficial for x = 0.15 (MLCC) to achieve good cycling stability27. The excellent temperature stability, frequency stability, and cycling stability have contributed to the widespread application of MLCCs devices.

For practical dielectric capacitor applications, the charge/discharge performance of MLCCs devices is also critical. The overdamped discharge performance for x = 0.15 (MLCC) under an electric field of 300-700 kV▪cm−1 was characterized using 1000 Ω as a load resistor while regular overdamped oscillatory waveforms were provided in Fig. 2f and Supplementary Fig. S7. In the meanwhile, it can be seen that the peak current Imax increases from 0.27 A to 0.63 A, and the actual maximum discharge energy density Wd increases from 1.26 J▪cm−3 to 5.94 J▪cm−3 with the increase of electric field (Supplementary Fig. S7). It is noteworthy that the t0.9 (90% of Wd is released) is only 1.34 µs at 700 kV▪cm−1, which suggests that the x = 0.15 (MLCC) demonstrates good charging and discharging properties, making it a highly promising material for pulsed power applications.

Average and local structural features

The evolution of the average and local structure is crucial for an in-depth elaboration of the reasons for the excellent ESP. The XRD spectra at room temperature (Supplementary Fig. S8) show that the samples display a characteristic perovskite structure when the BNH content is below 0.15 mol, which proves that BNH enters the NN main lattice and forms a solid solution. However, the second phase appears when the doping content is 0.20 mol, which indicates that BNH possesses a limited solid solubility in the NN matrix due to the different ionic radii and valence states28. XRD refinement for the different fractions was processed by FullProf software and the results are provided in Supplementary Fig. S9. It is widely acknowledged that NN exhibits both the AFE P phase and the FE Q phase in coexistence, and as the solid solubility of BNH increases, the phase structure gradually transitions from AFE P to AFE R, following previous research findings16. Temperature-dependent XRD (Supplementary Fig. S10) for x = 0.15 samples show no splitting peaks present or phase structure changes, implying excellent temperature stability. The HR-TEM was used to analyze the average structure in representative samples and demonstrated in Supplementary Fig. S11. 1/4 and 1/2 superlattice diffraction in selected area electron diffraction (SAED) concerning NN samples indicates the coexistence of AFE P (Pbcm) and FE Q (P21ma). Notably only the 1/4 superlattice diffraction patterns were found in the x = 0.05 sample, indicating that the AFE P phase was stabilized, which corresponds well to the XRD refinement results (Supplementary Fig. S9a, b). Interestingly, the weak stripe feature is found in the x = 0.15 sample, while no 1/6 superlattice diffraction pattern is observed, and therefore does not exhibit an AFE long-range ordered state (Supplementary Fig. S11h), which is in agreement with Zhang et al.17. In addition, an overview of bright-field images shows that the introduction of 0.05 mol BNH reduces the ___domain width of the ordered bands from ~500 nm to ~150 nm. However, we were unable to detect the nanoscale microregions in the x = 0.15 sample due to the limitation of the device resolution.

The high-angle annular dark-field scanning transmission electron microscopy (HAADF-STEM) was employed to analyze local structural properties, which can reveal the microstructural origins of strong local polarization. Precise atomic positions were captured and fitted by a two-dimensional Gaussian function, and polarization vectors were plotted according to the projected displacements of the B-site cations (Nb/Ni/Hf with weak intensity contrasts) relative to the centers of the crystal lattice of the four nearest A-site cations (Na/Bi with strong intensity contrasts), as indicated by the blue arrows and shown in Fig. 3a. Vectors with the same polarization direction are marked out through the red dashed box, and enlarged plots are provided in Fig. 3b–e, respectively. To better assess the distribution of polarization vectors, the displacement and angle of the polarization vectors are plotted and provided in Fig. 3f. Interestingly, we see unique O-type nanodomains with scales ranging from ~1 nm to ~3 nm. The presence of nanodomains breaks the long-range ordered antiferroelectric ordering while generating extensive nanoscale domains, which corresponds to the results of the phase-field calculation. These nanodomains are not segregated by a non-polar matrix, instead, they are interconnected by a polar matrix. A large number of polarized disordered matrices characterized by multiple localized disorders contribute to the transition between O-type nanodomains of different sizes, similar to slush-like polar structures29,30. In addition, the polar matrix with directionality, embodied by the average displacement of approximately 19.3 pm (Supplementary Fig. S12a) and a polarization vector angle of roughly 0–120° (Supplementary Fig. S12b), leads to the formation of directional slush-like polar structures (DSP), thus contributing to the simultaneous achievement of high Pmax and small Pr, which facilitates high Wrec and ultra-high η at the same electric field.

Fig. 3: The local structure analysis of the x = 0 and x = 0.15 ceramic.
figure 3

a Atomic-resolution HAADF-STEM polarization vector image of x = 0.15 along [001] direction. b–e Magnification of the marked areas in (a). f Projections of the ionic displacement polarization in (a).23Na MAS NMR spectra of g x = 0 and j x = 0.15. i, l Summation of the signals along δ1. h, k Summation of the signals along δ2.

Full size image

Solid-state nuclear magnetic resonance (NMR) spectroscopy equipped with a 9.4 T magnetic field was used to further analyze the local structural features31. The 23Na MAS NMR spectra about the x = 0 sample exhibit two distinct sodium sites, which are reflected in the positions of the two readily recognizable resonance peaks in the indirect dimension (δ1) at +5.91 ppm and −3.96 ppm, respectively (Fig. 3g, i), which are characteristic of the AFE P phase. Interestingly, a tiny shoulder peak (+5.22 ppm) appears next to the main resonance peak at the Na (1) site, which is the FE Q phase feature. As a contrast, only one signal was found in the x = 0.15 sample along the δ1 direction at −1.94 ppm (Fig. 3j, l), which suggests that a complete phase transition occurs. The position of the signal along the δ1 direction for the sample with x = 0.15 is similar to the 23Na NMR spectrum concerning NaNbO3 with the Pbnm space group17, which agrees with the experimental results of XRD and TEM. Alternatively, the distribution of the local structure for the NN matrix can be expressed by NMR parameters, e.g. the quadrupole coupling parameter (CQ) was used to characterize the degree of structural distortion to the NaO12 cuboctahedra32, and isotropic chemical shifts (δiso) were used to express the average Na-O distances33. Fitting the Na(1) signal of the x = 0 sample and the single signal of the x = 0.15 sample along the δ2 direction reveals that the CQ value decreases from 1043 kHz to 813 kHz, which suggests a decrease in the aberration of the NaO12 cuboctahedra (Fig. 3h, k, Supplementary Fig. S13, Supplementary Table S1), and the δiso decreases from −4.50 ppm to −6.42 ppm, which implies an increase in the Na-O bond distance. Thus, the entry of the guest-phase component 0.15 mol BNH into the main-phase NN matrix leads to a more disordered localized structure with smaller distortions at the Na (1) site. As CQ and δiso decrease, both γ and Wrec of the material increase substantially for x = 0 and x = 0.15 ceramics (Supplementary Fig. S14), which agrees with most NMR results for relaxing perovskite oxides possessing disordered localized structures17,31, further demonstrating that x = 0.15 component is indeed in a relaxed state, suggesting that more disordered localized structures possessing smaller distortions at the Na(1) site are favorable for energy storage performance. The localized disorder of the structure is more favorable for the flip of the ___domain structure under the action of an external electric field, which can reduce the Pr, achieving ultra-high η in energy storage applications.

Composition, field, and temperature drive characteristics

To further understand the component, electric field, and temperature-induced energy storage properties, (1-x)NN-xBNH and x = 0.15 (RRP) ceramics were prepared via traditional solid-phase sintering. We tested the uniaxial P-E loop at 250 kV▪cm−1 (Fig. 4a) and found that pure NN exhibits an irreversible AFE-FE transition16,20, leading to ferroelectric behavior under repeated electric field application. The hysteresis loss decreases dramatically, accompanied by the curve becoming elongated with the doping of the BNH guest phase, which corresponds to its calculated results from the phase field simulation method (Fig. 1b). To assess the overall ESP of the ceramics, we tested the uniaxial P-E loops about these ceramics under their breakdown electric field at 10 Hz (Supplementary Fig. S15). Interestingly, the shape of the P-E loop shifts from nearly square to thin stripes as the doping content increases, implying a transition of the system from the antiferroelectric to the relaxation state, so the optimum Wrec = 5.42 J▪cm−3 and η = 83.3% are obtained when x = 0.15 component. To further suppress the generation of holes during the sintering process, x = 0.15 powder was used to prepare ceramic tapes by repeating rolling process (RRP) to improve the ceramic densities and breakdown field strengths (Eb), which are conducive to the optimization of ESP. Figure 4b shows a P-E loop about x = 0.15 (RRP) ceramics from 100 to 820 kV▪cm−1 at room temperature. Strikingly, the P-E loop demonstrates a near-zero Pr during charging and discharging up to the breakdown field of 820 kV▪cm−1. Consequently, a high Wrec = 11 J▪cm−3 accompanied by an ultra-high η = 94.2% was simultaneously achieved for the x = 0.15 (RRP) ceramic at the highest electric field, outperforming most recently reported lead-free ceramics with excellent ESP (Supplementary Fig. S16). In addition, Eb first increases and then decreases with the doping of the BNH guest phase, and the Weibull distribution modulus β is greater than 15 (Fig. 4c), indicating improved reliability and uniformity regarding the ceramic samples. The enhancement of Eb can be attributed to two factors. On the one hand, the grain size gradually decreases with increasing doping content (Supplementary Fig. S17) for (1-x)NN-xBNH ceramics while Eb is inversely related to the grain size of the ceramics, which can be described as Eb ∝ (G)−α, where G is the grain size34. The x = 0.15 (MLCC) possesses the smallest grain size (1.5 µm) obtained by the linear intercept method (Supplementary Fig. S18), which is the main reason for its highest breakdown field strength. On the other hand, the significant lattice distortion caused by BNH doping in the guest phase increases the likelihood of electron collisions with lattice atoms, which can enhance electron scattering and raise electrical resistivity1. The grain boundary resistance and activation energy (Ea) gradually increase with increasing doping content, as demonstrated by the complex impedance spectra (Supplementary Fig. S19) and fitted plots (Fig. 4d), with the conductivity decreasing from 5.36 × 10−7 to 2.23 × 10−7 S▪cm−1 (Supplementary Fig. S20), which effectively reduces the conductive losses while favoring high efficiency and Eb of dielectrics.

Fig. 4: Ferroelectric properties, breakdown analysis, and complex impedance properties, along with relaxation behavior of NN and (1-x)NN-xBNH ceramics.
figure 4

a Unipolar P-E loops at 250 kV▪cm−1. b Room temperature P-E loop measured till the critical electric field of x = 0.15 (RRP) ceramic. c Weibull distribution of the breakdown strength on ceramic samples. d Arrhenius plots of total electric conductivity. e Temperature and frequency-dependence dielectric spectra. f Raman spectra at room temperature.

Full size image

Temperature- and frequency-dependent dielectric properties related to (1-x)NN-xBNH ceramics, as demonstrated in Fig. 4e and Supplementary Fig. S21, were characterized to probe the intrinsic mechanism of the polarization relaxation behavior. Pure NN ceramics possess a typical temperature-dependent dielectric anomaly peak (~390 °C) labeled Tm, representing the transition from the AFE P phase to the AFE R phase16. The Tm peak shifts to a position of ~195 °C when x = 0.05, exhibiting antiferroelectric characteristics at room temperature. However, the Tm peak moves rapidly to lower temperatures and a progressively broader temperature plateau is observed when x ≥ 0.10, which is characteristic of the relaxation state (Fig. 4e). Meanwhile, the calculated γ values for the (1-x)NN-xBNH ceramics lie in the range of 1.07–1.92 (Supplementary Fig. S22), implying an obvious relaxation behavior. In addition, the increasing full width at half-maximum (FWHM) and decreasing peak positions for ν5 and ν1 based on Raman spectra with increasing BNH content indicate a gradual elevation of localized component disorder concomitant with an enhancement of relaxation behavior, as shown in Fig. 4f and Supplementary Fig. S23. Our piezoelectric force microscopy (PFM) results (Supplementary Fig. S24) show that the ___domain size gradually decreases with increasing doping content, coinciding with localized compositional disorder, implying the fact that negligible Pr and hysteresis losses are present for x = 0.15 sample, thus increasing energy storage efficiency.

In summary, a rational phase-field simulation guided by the dielectric material design methodology is demonstrated to achieve nanoscale domains with overall superior ESP. The DSP strategy that we propose in the antiferroelectric system, forming highly dynamic slush-like polar structures with a certain degree of directionality, can significantly reduce the hysteresis loss and increase the breakdown field strength. Therefore, high Wrec = 14.1 J▪cm−3, ultra-high η = 96.8%, and figure of merit (UF = 440.6) accompanied by excellent temperature and cycling stability are realized in MLCCs. Guided by the roadmap for the design of directional slush-like polar structures, the DSP strategy provides a valuable research idea for broadening the exploration of high-performance relaxation ferroelectrics/antiferroelectrics and other related functionalities, accelerating their application to multilayer ceramic capacitors.

Methods

Sample preparation

(1-x)NaNbO3-xBi(Ni0.5Hf0.5)O3 with 0 ≤ x ≤ 0.20 ceramics were prepared by conventional sintering. Na2CO3 (>99.8%), Nb2O5 (>99.9%), Bi2O3 (>99.0%), NiO (>99.99%), and HfO2 (>99.9%) were weighed in specific proportions as starting materials. The initial powders and milling media were mixed and then ball-milled for 12 h. Then, the slurry was dried (4 h at 80 °C), calcined (4 h at 850 °C), and subjected to a secondary ball milling process under the same conditions as the initial ball milling. After drying the obtained slurry and passing it through a 200-mesh sieve, the resulting powder was pressed into green pellets with φ12 × 1 mm. Subsequently, the samples were sintered at 1190–1250 °C for 4 h. The sintered ceramic samples were polished to a thickness of 50–100 μm. Additionally, the x = 0.15 ceramic powder was mixed with 25 wt. % polyvinyl alcohol (PVA) and repeatedly rolled. After removing the PVA, the ceramic thick film was sintered at 1120–1190 °C for 2 h to obtain the x = 0.15 (RRP) samples.

MLCCs preparation

The NN-based MLCCs with matrix powder (x = 0.15) using 0.4 wt.% MnO2 (>99.95%) as a sintering aid were fabricated using the same preparation process as that used for bulk ceramics. The matrix powder was mixed with benzyl butyl phthalate and polypropylene carbonate using an ultra-high-speed mixer at a speed of 1500 r▪min−1 for 2 h, with methyl ethyl ketone as the milling medium, to obtain the initial slurry. The slurry was uniformly cast onto the film tape using a tape-casting machine. After the film tape was dried, it was coated with Pt electrodes, which served as internal electrodes, and then cut into pieces that were stacked into six dielectric layers using a mold, followed by hot pressing at 80 °C for 20 min. The MLCCs were heated to 360 °C for 3 h to remove the binder, followed by sintering at a temperature range of 1100–1150 °C for 2 h to obtain x = 0.15 (MLCC) samples. MLCCs with Pt as the inner electrode have a single-layer electrode area of about 3 mm2. The silver paste was used as the termination electrode for electrical measurements.

Structural characterization

The phase structure was analyzed using an X-ray diffractometer (XRD, SmartLab-3kW, Tokyo, Japan). The microstructure was examined using a scanning electron microscope (SEM, FEI Quanta 250 FEG, Hillsboro, USA). The ___domain structure was analyzed using transmission electron microscopy (TEM, JEM-2100CX, Japan) and piezoelectric force microscopy (PFM, MDP-3D, USA). The HAADF atomic-scale images were acquired by STEM and processed by 2D Gaussian fitting. The vibrational modes were characterized using Raman spectroscopy (Invia, Renishaw) with a 532 nm excitation source.

Electrical properties

The dielectric properties and impedance spectra were measured using a custom-built inductor-capacitance-resistance LCR meter (E4980A, Agilent, USA). The P-E loops were measured using a ferroelectric analyzer (TF Analyzer 2000, Aachen, Germany). The dielectric breakdown behavior was analyzed using the Weibull distribution. The charge and discharge performance was evaluated using a charge-discharge testing instrument (Tongguo Technology, China).

NMR characterization

Solid-state NMR spectra were characterized using a Bruker 400 Avance III spectrometer, equipped with a wide-bore 9.4 T magnet. The powder sample was placed in a conventional 4 mm ZrO2 rotor and subjected to magic angle spinning (MAS) at 10 kHz. The 23Na MAS NMR spectra were obtained using a single-pulse experiment under a 9.4 T magnetic field, with a typical pulse length of 3.19 µs.

Phase-field simulations

In phase-field models, the phase transition from paraelectrics to antiferroelectrics in the NaNbO3-based system has been established. The total free energy of the whole system can be described as follows:

$$F={\int }_{V}( \, \, {f}_{{lan}}+{f}_{{grad}}){dV}+{\int }_{V}( \, \, {f}_{{elas}}+{f}_{{elec}}){dV}+{\int }_{V}{f}_{{localelec}}{dV}$$
(1)

where V is the volume of the system.

The landau free energy flan can be expressed in terms of P as:

$${f}_{lan} ={A}_{1}({{{{\rm{P}}}}}^{2})-{A}_{11}{({{{{\rm{P}}}}}^{2})}^{2}+{A}_{12}\left({\sum }_{i,j,i\ne j}^{3}{P}_{i}^{2}{P}_{j}^{2}\right)\\ \quad+ {A}_{111}{({{{{\rm{P}}}}}^{2})}^{3}+{A}_{112}\left({\sum }_{i,j,i\ne j}^{3}{P}_{i}^{4}{P}_{j}^{2}\right)+{A}_{113}\left({P}_{1}^{2}{P}_{2}^{2}{P}_{3}^{2}\right)\\ ={A}_{1}\left({P}_{1}^{2}+{P}_{2}^{2}+{P}_{3}^{2}\right)-{A}_{11}{\left({P}_{1}^{2}+{P}_{2}^{2}+{P}_{3}^{2}\right)}^{2} \\ \quad+{A}_{12}\left({P}_{1}^{2}{P}_{2}^{2}+{P}_{2}^{2}{P}_{3}^{2}+{P}_{1}^{2}{P}_{3}^{2}\right) +{A}_{111}{\left({P}_{1}^{2}+{P}_{2}^{2}+{P}_{3}^{2}\right)}^{3} \\ \quad+{A}_{112}\left({P}_{1}^{4}{P}_{2}^{2}+{P}_{2}^{4}{P}_{3}^{2}+{P}_{1}^{4}{P}_{3}^{2}+{P}_{1}^{2}{P}_{2}^{4}+{P}_{2}^{2}{P}_{3}^{4}+{P}_{1}^{2}{P}_{3}^{4}\right)\\ \quad+{A}_{113}\left({P}_{1}^{2}{P}_{2}^{2}{P}_{3}^{2}\right)$$
(2)

where Aij is the Landau coefficient related with the dopant’s concentration and temperature.

The gradient-free energy fgrad representing the energy from polarization inhomogeneity is described as:

$${f}_{{grad}}={\sum }_{i}\left\{\begin{array}{c}{\alpha }_{1}\left[{\left(\frac{\partial {P}_{i}}{\partial x}\right)}^{2}+{\left(\frac{\partial {P}_{i}}{\partial y}\right)}^{2}+{\left(\frac{\partial {P}_{i}}{\partial z}\right)}^{2}\right]\\+{\alpha }_{2}\left[{\left(\frac{{\partial }^{2}{P}_{i}}{\partial {x}^{2}}\right)}^{2}+{\left(\frac{{\partial }^{2}{P}_{i}}{\partial {y}^{2}}\right)}^{2}+{\left(\frac{{\partial }^{2}{P}_{i}}{\partial {z}^{2}}\right)}^{2}\right]\end{array}\right\},\,(i=1,2,3)$$
(3)

where α1 is a negative coefficient representing the stability of the antiferroelectric phase, α2 is a positive coefficient representing the stability of the ferroelectric phase.

The elastic energy felas can be written as

$${f}_{{elas}}=\frac{1}{2}{C}_{{ijkl}}{e}_{{ij}}{e}_{{kl}}=\frac{1}{2}{C}_{{ijkl}}({\varepsilon }_{{ij}}-{\varepsilon }_{{ij}}^{0})({\varepsilon }_{{kl}}-{\varepsilon }_{{kl}}^{0})$$
(4)

where Cijkl is the elastic stiffness tensor, εij is the total strain, ε0kl is the electrostrictive stress-free strain, i.e., ε0kl=QijklPkPl.

The electrostatic energy felec can be expressed as:

$${f}_{elec} ={f}_{dipole}+{f}_{depol}+{f}_{appl}\\ =-\frac{1}{2}{E}_{i}{P}_{i}-\frac{1}{2}{E}_{i,depol}\overline{{P}_{i}}-\frac{1}{2}{E}_{i,appl}\overline{{P}_{i}}$$
(5)

where fdipole, fdepol, and fappl are the dipole-dipole interaction energy from polarization, the depolarization energy, and the energy from the applied electric field respectively. Ei represents the inhomogeneous electric field caused by dipole-dipole interactions, Ei, depol is the mean depolarization field arising from the surface charges, Ei, appl is expressed by a sine function.

The local random distribution electric field energy flocalelec related with dopant valence and concentration can be expressed as:

$${f}_{{localelec}}={E}_{i,{localelec}}\bar{{P}_{i}}$$
(6)

where Ei,localelec is the random distribution electric field, \(\bar{{P}_{i}}\) is the spatial average polarization.

The temporal evolution of the polarization of the antiferroelectric system can be calculated by solving the time-dependent Ginzburg-Landau (TDGL) equation:

$$\frac{\partial {P}_{i}(x,t)}{\partial t}=-L\frac{\delta F}{\delta {P}_{i}(x,t)},\,(i=1,2,3)$$
(7)

where L is the kinetic coefficient related to the ___domain mobility, t is time, Pi(x, t) is the polarization and F is the total free energy of the whole system.

The phase-field simulations were carried out in two dimensions with cell sizes 128Δx × 128 Δy (Δx and Δy are the number of grid points which are equal to 2 nm) with periodic boundary conditions.