Figure 1(a) shows the room-temperature XRD patterns of the ceramics with different *x* (*x* = 0.00, 0.02, 0.04, 0.06 and 0.08), and the enlarged diffraction peaks in the 2*θ* range of 31-33o and 38–40° can be seen in Fig. 1(b) and 1(c), respectively. Although a stable ABO3 perovskite structure has been formed in all samples, a small amount of Bi-rich impurities (Bi24Fe2O39) can be also detected. Obviously, both (110) and (111) diffraction peaks shift to lower angles as *x* increases, indicating an increase in the lattice parameters when the larger radii ions Mg2+ (0.72Å, CN = 12) and Zr4+ (0.72Å, CN = 12) enter the lattice and substitute the Fe3+ (0.64Å, CN = 12) and Ti4+(0.604Å, CN = 12) sites. In order to further identify the phase structure of BF-BT-*x*BMZ ceramics, well Rietveld refinement with low Sig and Rw values has been conducted based on rhombohedral (*R*) and pseudocubic (*P*C) phases, as shown in Fig. 1(d). And the *P*C phase gradually dominant with the increase of BMZ components while *R* phase decreases [Fig. 1(e)].

Figure 2(a) presents the surface SEM images of the well-polished and thermally etched ceramics. All samples exhibit irregular lumps with clear grain boundaries. Compared to *x* = 0, samples modified with BMZ display much lower porosity and enhanced density. In addition, the grain size distribution is presented in the Fig. 2(b). The average grain size of *x* = 0 is 7.41 *µ*m. And the grain size first decreases to 2.74 *µ*m and then increases to 4.71 *µ*m with increasing BMZ content. The decrease in the grain size may be related to the solute drag mechanism. It’s well known that the difference in ionic radii between dopant and host ions can induce lattice strain energy (∆*G*strain), which can be described according to the following equations[18]:

where *r*0, *r*d, *N*A, and *M* represent the optimal radius of the lattice site, ionic radius of dopant, Avogadro’s number, and Young’s modulus, respectively. In this study, the ionic radius of Mg2+ (0.72Å, CN = 12) and Zr4+ (0.72Å, CN = 12) are larger than those of Fe3+ (0.64Å, CN = 12) and Ti4+(0.604Å, CN = 12). Therefore, as the concentration of Fe3+ and Ti4+ increase, ∆*G*strain becomes larger, which inhibits the grain boundary mobility and suppresses the grain growth. When *x* > 0.04, the excessive Mg2+ into the *B*-sites could results in the generation oxygen vacancies which could help the mass transport during sintering process and then promote the grain growth.[19] Therefore, appropriate introduction of BMZ can promote the grain growth and benefit the enhancement of density. As presented in Fig. 2(c), all elements (Bi, Fe, Ba, Ti, Mg, Zr, O) are distributed with good chemical uniformity for *x* = 0.04.

To establish the relationship between the crystal structure and dielectric properties, the temperature-dependent dielectric constant (*ε*r-*T*) are displayed for BF-BT-*x*BMZ ceramics in Figure S1, measured at 25 to 600 ℃, and 1, 10, 100, and 1000 kHz. Strong frequency dispersion and broaden peaks are obtained around *T*m with introduced BMZ. To characterize the diffuseness/relaxation degree, four different relaxation calculation methods were adopted, as shown in Fig. 3. According to the Curie-Weiss law for the dielectric response above *T*C, the temperature-dependent *ε*r is expressed in different forms around *T*m or *T*C intercepted from the *ε*r-*T* curves measured at 300–600°C to characterize the relaxation degree. Firstly, it is noted that *ε*r follows the Curie-Weiss law above 550°C, as shown in the Fig. 3(a). The relaxor *γ* is the slope of the fitting lines for the logarithmic plot of (1/*ε*r-1/*ε*m) versus (*T*-*T*m), according to the equation:[20]

where *ε*r, *ε*m and *C* represent dielectric constant, maximum dielectric constant and the Curie constant, respectively. In general, *γ =* 1 and *γ* = 2 represent normal ferroelectrics and ideal relaxors, respectively, while the values between 1 and 2 denote relaxor ferroelectrics. Obviously, the *γ* increases as BMZ content elevates, as shown in Fig. 3(b). And, the normalized *ε*r/*ε*m curves are obtained as a function of *T*-*T*m based on *ε*m and *T*m. Then *W*2/3-H and *W*2/3-L are gained as the difference between the *T*m and the temperature where *ε*r reaches 2/3 of the *ε*m from the high-temperature and low-temperature sides, respectively, as shown in Fig. 3(c).[21–22] Because of the Curie-Weiss law is applied above *T*C or *T*m, the diffuseness temperature (*T*diffuseness) equals *W*2/3-H instead of *W*2/3-L. It is shown that the *T*diffuseness gradually increases with increasing BMZ, as shown in Fig. 3(e). One more method is to calculated *T*diffuseness according to a Gaussian distribution of *T*C. *ε*r as a function of the temperature is expressed as

After expanding and simplifying with the limitation of 1 < *ε*m/*ε*r < 1.5, where *δ*g is the *T*diffuseness.[22–24] According to the fitting results for *T* ≥ *T*m(*T*C) [Figure 3(d)], increasing *δ*g values could be gained, exhibiting the same change tendency to the value of *W*1/2-H, as presented in Fig. 3(e). Thus, it is confirmed that the relaxation degree increases with elevated BMZ component. According to previous studies, relaxation is mainly induced by structural disorder, heterogeneous microscopic distribution and stoichiometric fluctuation, et al.

Obviously, the more BMZ, the stronger dielectric relaxation. The addition of BMZ can break the long-range ferroelectric orders, leading to the domination of pseudo-cubic phase and strengthened dielectric relaxation. With strengthened random field, the long-range ordering is gradually destroyed, contributing structure/composition heterogeneity. And refined domain structure and enlarged polar nanoregions are formed. Therefore, the relaxation degree enhances with increases components as well as the phase transition from ferroelectric to relaxor states.

The structure evolution can greatly influence electrical properties of BF-BT-BMZ ceramics. The *P*-*E* loops and *S*-*E* curves are shown in Fig. 4(a-c), measured at 50 kV/cm, 1 Hz and 100 ℃. Figure 4(d-e) present the electrical parameters as a function of BMZ content, including remanent polarization (*P*r), coercive field (*E*C) and strain (*S*pos and *S*neg derived from bipolar *S*-*E* curves, *S*uni derived from unipolar *S*-*E* curves). The *P*r decreases as the *x* the increases while the *P*-*E* loops become slim, indicating the increased relaxor ferroelectric states with lowered polarization and refined domain.[25] Significantly, the strain is optimized with modified structure. Butterfly *S*-*E* curves gradually change into sprout *S*-*E* curves when *x* increases, and typical unipolar strain curves are realized for the ceramics. Both *S*pos and *S*uni increase first and then decrease with increasing *x*, as well as the maximum values (*S*pos = 0.24% and *S*uni = 0.25%) for *x* = 0.04. The enhanced strain is due to the coexistence of long-range order matrix and PNRs which can promote the domain switching and reversion of poled domains. Meanwhile, the poor strain is attributed to excess PNRs and little long-range order ferroelectric matrix from further elevated BMZ content. Obviously, long-range order ferroelectric matrix with appropriate PNRs can strengthen the strain effectively with large polarization and easy domain evolution. Notably, the negative strain (*S*neg) gradually decreases with elevating *x* because of the reduced irreversible non-180o domain from lowered ferroelectric *R* phase and increased *P*C phase.

It is well known that good temperature stability is necessary in the actuator applications. Therefore, the temperature stability of electrical properties is further explored for the ceramics. The temperature-dependent strain properties are measured at *E* = 60 kV/cm, 30 to 120 ℃, as shown in Figure S2 and Fig. 5. The strain values are presented in Figure S2 (e), (f) and (g). As typical bipolar strain curves maintain with increasing temperature [Figure 5(a-c) and Fig. S2(a-b)], both *S*pos gradually increases with elevating temperature, which is due to promoted piezoelectric effect and thermally activated domain motion/switching at thermal field for ferroelectric state, as shown in Fig. S2(f). *S*neg exhibits a decreasing tendency when *x* = 0.02–0.08 which is attributed to the weakened non-180o domain switching.[26–27], while an enhancement is found when *T* ≥ 100℃ for *x* = 0 [Fig. S2(e)]. And typical unipolar strain curves are obtained with increasing temperature, as shown in Fig. 5(d-f). Obviously, *S*uni increases with elevating temperature due to promoted domain switching under thermal field [Fig. S2(g)], which is similar to the change tendency of *S*pos. In order to better compare the thermal stability of different components, parameter *S*/*S**T*=30 ℃ are presented with increasing temperature, as shown in Fig. 5(g-i). Obviously, weakened enhance of *S*pos and *S*uni as well as decline of *S*neg is found with increasing temperature for *x* = 0.02–0.08 with respect to that for *x* = 0, implying improved temperature stability. This phenomenon is due to the weakened polarization enhancement at increased temperature from the refined domain size and strengthen relaxor behavior as well as induced *P*C phase.

Figure 6. (a-c) and Figure S3 exhibit the *P*-*E* loops measured at 30 to 120 ℃. *P*max and *P*max/*P*max *T*=30℃ increase with increasing temperature [Fig. 6(d) and Fig. S3(c)], which is due to the thermal field induced promoting domain switching.[28–29] For ceramics dominated with *R* phase (*x* = 0-0.02), the *P*r and *P*r/*P*r *T*=30℃ increase with the increasing temperature [Fig. S3(d)], indicating that the thermal field is beneficial to promote non-180o domain switching, which is a typical characteristic for BF-based ferroelectrics. For the ceramics with *x* = 0.04–0.06, the *P*r and *P*r/*P*r *T*=30℃ slightly fluctuate with increasing temperature, indicating the weaking of non-180o domain switching and the increasing of nanodomain. And a decline is observed for *x* = 0.08 at high temperature due to the strong relaxor behavior. For all ceramics, the *E*C and *E*C/*E*C *T*=30℃ decrease with increasing temperature [Fig. 6(f) and Fig. S3(e)], indicating the promoted domain switching. Thus, improved temperature stability could be gained with the coexistence of macrodomains and polar nanoregions, yielding to appropriate realxor behavior.

Figure 7 illustrates the amplitude images of the ceramics. Notably, the domain morphology depends on grain with obvious grain boundary, which conforms to the result of reported BF-based ceramics.[30] Combing with the phase images [Fig. S4], the initial state with dominated randomly distributed long-range ordered domains is formed for *x* = 0 [Fig. 7(a1) and (a2)]. As relaxor behavior is strengthened with increasing *P*C phase (*x* = 0.04), lowered ferroelectric domain and increased PNRs is produced, which could induce the reduced domain size [Fig. 7(b1) and (b2)]. When *P*C phase further increases (*x* = 0.08), relaxor ferroelectric state with dominated PNRs is revealed, exhibiting seriously deteriorative ferroelectric domain [Fig. 7(c)]. And enhanced piezoelectric response is achieved for the ceramics when *x* = 0.04, manifesting the promoted domain switching and lowered energy barrier from the ferroelectric-relaxor state [Fig. 7(d)~(f)].