Phase field simulation of de-aging process in acceptor-doped ferroelectrics

Journal Pre-proof Phase field simulation of de-aging process in acceptor-doped ferroelectrics Chao Yang, Enwei Sun, Zhen Liu, Yunfei Chang, Bin Yang, Wenwu Cao PII:

S0925-8388(19)33749-1

DOI:

https://doi.org/10.1016/j.jallcom.2019.152503

Reference:

JALCOM 152503

To appear in:

Journal of Alloys and Compounds

Received Date: 1 July 2019 Revised Date:

23 September 2019

Accepted Date: 28 September 2019

Please cite this article as: C. Yang, E. Sun, Z. Liu, Y. Chang, B. Yang, W. Cao, Phase field simulation of de-aging process in acceptor-doped ferroelectrics, Journal of Alloys and Compounds (2019), doi: https:// doi.org/10.1016/j.jallcom.2019.152503. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Phase field simulation of de-aging process in acceptor-doped ferroelectrics Chao Yang a, Enwei Sun a,b,*, Zhen Liu a, Yunfei Chang a,b, Bin Yang a,b, Wenwu Cao a,b,c,*

,

a

Condensed Matter Science and Technology Institute and Department of Physics, Harbin Institute

of Technology, Harbin 150080, China b

School of Instrumentation Science and Engineering, Harbin Institute of Technology, Harbin

150080, China c

Department of Mathematics and Materials Research Institute, The Pennsylvania State

University, University Park, Pennsylvania 16802, USA *

Corresponding authors: [email protected] (E. Sun), [email protected] (W. Cao)

Abstract: A phenomenological de-aging model with uniformly distributed defect dipoles for acceptor-doped ferroelectrics is proposed with the dynamic switching process of defect dipoles being simulated by the kinetic Monte Carlo method. The shifted P-E loops at different frequencies and temperatures were theoretically investigated. The results showed that the internal bias field increases with frequency and decreases with temperature. The scaling relation of the internal bias field with frequency was calculated at 298 K, 318 K and 338 K. The number of defect dipoles along different directions was calculated during the de-aging process, which is a key factor affecting the shape of P-E hysteresis loop. At 0.1 Hz, the defect dipoles can align with the polarization direction, but when the frequency is more than 10 Hz, they could not follow the field change and the internal bias field will saturate to a fixed 1

value as the field frequency is beyond 100 Hz. Keywords: acceptor-doped ferroelectrics, de-aging process, defect dipoles, phase field simulation, Monte Carlo method

1. Introduction Acceptor-doped ferroelectric materials have a wide range of applications in high power electromechanical devices due to their large coercive field, low dielectric loss and high mechanical quality factors [1-3]. In general, the domain wall stability in acceptor-doped ferroelectric materials is higher than undoped ferroelectrics due to the existence of an internal bias electric field [4-7]. This is why the acceptor-doped ferroelectric materials generally show more recoverable electro-strain [8, 9]. However, for acceptor-doped ferroelectric materials, the occurrence of de-aging process is inevitable since the aligned defect dipoles may become randomized. The de-aging process is accompanied by the decrease in coercive field and domain stability, which will reduce the electric field stability of the material and make the materials gradually fails to meet the demand for high-power devices [10]. Therefore, establishing a dynamic model for the de-aging process is practically useful for many applications of acceptor-doped ferroelectric materials. In general, there are many defect dipoles distributed inside the acceptor-doped ferroelectrics. These defect dipoles cause the acceptor-doped ferroelectric to exhibit “hard”

characteristics

as

mentioned

above 2

[7,

11-14].

Ren

proposed

a

symmetry-conforming short-range ordering (SC-SRO) theory to explain the alignment of defects [15]. After aging, defect dipoles in the acceptor-doped ferroelectric materials are always aligned with the spontaneous polarizations [16-18]. The electric fields caused by these aligned defect dipoles are superimposed to form a strong internal bias electric field. The internal electric field will make it harder for the external stimuli to disturb the established domain configuration [19, 20]. While in the de-aging process, defect dipoles gradually become randomly oriented, i.e., to drift away from the aligned state. The randomization of defect dipole orientations will cause the internal electric field to become smaller. The weakened internal bias electric field will make the domain walls less stable, hence, increases the mechanical loss (or reducing the mechanical Qm). Therefore, the dynamic process of de-aging is directly associated with the dynamic behavior of defect dipoles. However, previous models about defect dipoles fixed the orientation of defect dipoles in order to study the steady state of acceptor-doped ferroelectrics, for instance, the study of P-E loops after sufficient aging [21-24]. In the de-aging process, it is necessary to know the time evolution of defect dipoles with applied external electric field [25, 26]. In fact, in acceptor-doped ferroelectrics, the rotation of defect dipoles is attributed to the hopping of oxygen vacancies. Defect dipoles could form by the combination of oxygen vacancies and doped ions. For a tetragonal phase ferroelectric material, oxygen vacancies can only stay at the six vertices of the oxygen octahedron, so there are six orientations for the defect dipoles [27-29]. In terms of energy, these six directions are not equivalent under field. By knowing the energy states of oxygen 3

vacancies, their hopping probabilities at different positions can be calculated. The dynamic model of defect dipoles can be established using this strategy. In the process of de-aging, the dynamic change of defect dipoles can be affected by external factors [1]. In order to quantitatively study the influence of these factors on de-aging process, it is necessary to know the transition rate of oxygen vacancies and their effect on the de-aging rate. For the hopping of oxygen vacancies, the time spends from 10 -5 s to 10 6 s at 380 K. When the temperature is reduced, the needed time will increase significantly. For such a time scale, the kinetic Monte Carlo (kMC) method is an effective tool to simulate the process [30, 31]. In addition, the influence of domain structure on the hopping of oxygen vacancies should also be taken into account. Therefore, in this work, by coupling the evolution of the defect dipole to the spontaneous polarization, we developed a two-dimensional model consisting of domain evolution with rotatable defect dipoles to simulate the de-aging process of acceptor-doped ferroelectric materials. Phase field method has been used for the study of

domain

structure

evolution

in

ferroelectrics

employing

the

Landau-Ginzburg-Devonshire (LGD) theory, and the dynamic behavior of defect dipoles is modeled using the kinetic Monte Carlo method. The model here simultaneously takes into account the long-range electric and elastic interactions. 2. Theoretical model The phase field and kMC method are employed here to form our simulation

4

model. The spontaneous polarization vector Pi ( i = 1,2,3) is chosen as the primary order parameter. The time-dependent Ginzburg-Landau equation [32] is given by

∂Pi δF = −L δ Pi ∂t

(1)

where L is a kinetic coefficient associated with the system relaxation rate, F is the total free-energy and is written in the form of polarization components, the partial derivative of polarization Pi, j = ∂Pi / ∂x j and elastic strains ε ij =

1 ( ui, j + u j ,i ) , 2

i, j = 1, 2,3 , which can be written as:

F = ∫ dr fGLD { Pi , Pi , j , ε ij } + Fdip { Pi } + Fdef { Pi }

(2)

where the first term in the bracket is the bulk free-energy density, which takes the following form [33-35]: f GLD {Pi , Pi , j , ε ij } = α1 ∑ Pi 2 + α11′ ∑ Pi 4 + α12′ ∑ Pi 2 Pj2 + α111 ∑ Pi 6 i

i> j

i

i

(3)

+ α P P P − PE i i + 1 2 cijkl ε ij ε kl − qijkl ε ij Pk Pl 2 2 2 123 1 2 3

+ 1 2 Gijkl Pi , j Pk ,l

The coefficients α1 , αij′ , and α ijk are Landau coefficient, cijkl , qijkl and Gijkl stand for elastic stiffness constants, electrostriction coefficients and gradient energy coefficients, respectively. In the simulation, the Voigt notation is used [36]: q11=q1111, q12=q1122, q44=2q1122, c11=c1111, c12=c1122, c44=c1212, g11=G1111, g12=G1122, and g44=G1212.

Ei stands for the component of the external electric field. The zero-strain coefficients

αij′ can be expressed as follows[36]: 1  qˆ 2

2qˆ 2 

α11′ = α11 +  11 + 22  6  Cˆ11 Cˆ 22 

5

(4)

2 2  1  2qˆ112 2qˆ22 3q44 α12′ = α12 +  − +  6  Cˆ11 Cˆ 22 C44 

(5)

where cˆ11 = c11 + 2c12

cˆ22 = c11 − c12 qˆ11 = q11 + 2q12

qˆ22 = q11 − q12 The electrostatic energy due to spontaneous polarization is written as [37]: Fdip { Pi } = − 1 2 ∫ dr  E dip (r ) ⋅ P (r ) 

(6)

where E dip ( ri ) =

 P ( r ) 3 ( ri − r j )  P ( r j ) ⋅ ( ri − r j )    j   3 d − r  j  3 5 ∫ 4πε 0 ri − r j  ri − r j  1

(7)

denotes the inhomogeneous depolarization field due to dipoles. The last term in Eq. (2) is the free energy due to defect dipoles, which take the form of: Fdef (Pi ) = − ∫ dr (E def ⋅ P (r ))

(8)

where Edef = −∇Vdef is the internal field associated with defect dipoles [38]. nd  q0 (r j ) q0 (r j )  Vdef (ri ) = ∑  −  ( ) r − r − δ j  ri − (r j + δ j )  i j j 

(9)

The related parameters in Eq. (9) are explained in detail in Ref. [21]. The parameters of Pb(Zr0.2Ti0.8)O3 are used in this work and we used normalized coefficients for the convenience of calculations. All coefficients are listed in Table 1. P0 = 0.8 Cm −2 is spontaneous polarization at room temperature.

6

Table 1 Values of the normalized coefficients used in the simulation.

α11( e)′

α12( e)′

1.3

2.1

′ α112 ′ α111

′ α123

g14

g 44

q11′

q12′

′ q44

c11′

c12′

′ c44

0.40

-8.0

0.0

1.0

70

1.7

48

2400

1090

1527

1.6

Previous theoretical models in the literature assumed that defect dipoles absorbed energy from the electric field [39]. Once defect dipoles got enough energy they could jump from one local minimum to another. But this theory cannot provide quantitative results. In fact, defect dipoles are different from the spontaneous polarization in ferroelectrics because they don’t directly absorbed energy from the electric field [27]. According to first principle calculations, the energy barriers of oxygen vacancy hopping are deeply affected by the polarization direction, while the influence of the external electric field on the energy barrier is rather small [27]. So the “direct” influence of electric field on the defect dipoles could be ignored, while the polarization could be changed by an external electric field, which would then change the energy state of defect dipoles. So the electric field has an “indirect” influence on defect dipoles. This mechanism is not popular in the literature but it catches the essence of defect dipoles [30, 39]. According to the thermodynamic theory, the transition rate of state i jumping to state j is expressed by the equation:  -∆E (mi , j )  k (i , j ) = k0( i , j ) exp   k T   B 

(10)

where k 0( i , j ) is the prefactor, which is assumed as 2 THz in our calculation [27]. The value represents the characteristic frequency for atomic vibrations, k B is the 7

Boltzmann constant, ∆Em(i , j ) is the energy barrier associated with the transition. In this paper, the kMC method is conducted as follows: At each kMC step, three possible transition ways are prepared. The three possible ways of transition correspond to the other three local minima on the energy profile for the 2D model. The transition rates are controlled by Eq. (10). Then, one transition way is randomly chosen

within

the

three

possible

ways,

meeting

the

condition

of

n=m

K (m−1) ≤ r1K ( N ) ≤ K (m) , where K (m ) = ∑ k (n ) is the spatial sums of the rates and n =1

K(

N)

is the sum of all rates. The time is updated by a step ∆t =

1 1   r ( N ) ln  . 1 and K  r2 

r2 are two random numbers distributing within [0,1] with a uniform probability. Finally, put the oxygen vacancy associated with this transition probability to its new site and go to the next kMC step. In fact, even a well-polarized single crystal ferroelectric could not maintain single domain in nature [40]. So a ferroelectric system of multidomain should be set as initial state before an alternating electric field is applied. The initial state of the domain is set as shown in Fig. 1(a). After the de-aging process, the direction of the defect dipoles will be randomized as shown in Fig. 1(b). The domain size is 45 nm in our simulation. In Fig. 1, the big arrows indicate the directions of spontaneous polarizations and the small arrows indicate defect dipoles in the [01], [10], [01] and [10] directions. Under the influence of AC electric field along [11] and domain patterns, defect dipoles along [01] and [10] are equivalent, while those along [01] and [10] are equivalent under a field along [11], so we will not distinguish [01] and [10] 8

(or [01] and [10] ) in our calculations. Those defect dipoles are randomly located within domains, and are aligned with the spontaneous polarizations according to the SC-SRO theory [11]. There are only defect dipoles along [01] and [10] in the initial state, as shown in Fig. 1(a). The simulated concentration of defects is chosen as 5%, so there are total 819 defect dipoles in the area of 128×128 grid points.

Fig. 1. Schematic diagram of the domain structures and defect dipoles, (a) the initial domain state, (b) domain state after de-aging process, (c) the external AC electric field for two circles, and (d) defect dipole directions in two dimensions.

Periodic boundary conditions are applied here to account for long range interactions. Alternating sinusoidal electric field of the form E = E0 sin(2π ft − π / 2) is applied along [11] direction for two circles, as shown in Fig. l(c), where E0 is set at 13.2 kV/cm. For the convenience of computer simulation, we discretize the electric field by m steps, so the time interval of per step of the electric field is tw=T/m, where T is the total applying time of electric field. The principle of calculation is shown in Fig. 2. First, a step of electric field is applied, then the spontaneous polarization dipoles 9

evolve rapidly according to the time-dependent Ginzburg-Landau equation, i.e. Eq. (1). Subsequently, the defect dipoles evolve according to the kMC method mentioned above. Since the time required for the evolution of defect dipoles is much larger than that of spontaneous polarization dipoles, it can be considered that the evolution of the spontaneous polarization dipoles is completed immediately while the defect dipoles are constantly evolving after each step of the electric field. After each step of electric field, the total time cost by the evolution of defect dipoles is compared with tw, if the total time is larger than tw, the next step of electric is applied, otherwise, defect dipoles continues to evolve.

Fig. 2. Flow chart of electric field applied method and defect dipole evolution.

3. Simulation results and discussion The shifted P-E loops at room temperature (298 K) were simulated under an AC 10

electric field with frequencies ranging from 0.1 Hz to 100 Hz, as shown in Fig. 3. When the frequency increases from 0.1 Hz to 100 Hz, Ec+ increase from 0.42 to 0.69 in normalized scale, and Ec- decreases from 0.34 to 0.20 as shown in Fig. 3, which indicates the intense offset of P-E loops at high frequencies. The frequency affects the defect dipole switching ratio, while the switching ratio has a major impact on the shape of the P-E loops.

Fig. 3. The shifted P-E loops at 298 K under AC electric field with various frequencies for acceptor-doped PZT system by phase field simulation.

In order to verify the effect of the defect dipole switching ratio on the P-E loops, we simulated the time evolution of defect dipole percentage along different directions. Every curve in Fig. 4 consists of thousands of points and each point indicates percentage of defect dipoles at that moment. The simulations were performed at frequencies of 0.1, 1, 10, and 100 Hz at 298 K. The red line denotes defect dipoles at [01] and [10], which is zero at the beginning of the simulation. When the electric field is changed from [11] to [11], the polarizations are switched to [01] and [10] with certain degree of phase lagging, then, defect dipoles oriented along [01] and [10] will 11

try to follow. Because their response time is longer, only part of them may be switched to [01] and [10]. At 0.1 Hz, the field change is slow enough to allow defect dipoles to follow, so that nearly all of them switched to [01] and [10], as shown in Fig. 4(a). The two lines cross with each other at the end of simulation in Fig. 4(a), which indicates that the hysteresis loop is almost symmetric after the second electric field circle. At 1 Hz and 10 Hz (shown in Fig. 4(b) and Fig. 4(c)) not all defect dipoles can follow the electric field switching. The oscillation of the curve in Fig. 4(b) is more intense than that in Fig. 4(c), indicating the more number of switched defect dipoles at 1 Hz than at 10 Hz. Those switched defect dipoles will contribute to the de-aging effect. At 100 Hz (shown in Fig. 4(d)), the two curves are nearly unchanged with the electric field. This means that the defect dipoles were unable to follow the field change, so their effect would become be to provide an internal bias field, causing the shift of the hysteresis loop to one direction. The simulation results clearly demonstrated the graduate formation process of the internal bias with the increase of frequency.

12

Fig. 4. The time evolution of defect dipole percentage at various directions and frequencies, (a) 0. 1 Hz, (b) 1 Hz, (c) 10 Hz, and (d) 100 Hz.

The domain structure configurations with rotatable defect dipoles at low frequency (0.1 Hz) and high frequency (100 Hz) are simulated to reveal the correlation between the defect dipole switching degree and the microstructure of ferroelectrics as shown in Fig. 5. Fig. 5(a) is the initial state, which is consistent with that in Fig. 1. When the electric field is at the maximum along [11], intrinsic dipoles in all domains are either along [01] or [10] directions and all defect dipoles are along the spontaneous polarization directions in each domain. Fig. 5(b) and Fig. 5(d) denotes the domain configurations when the electric field is at the maximum along [11]. The frequency is 0.1 Hz in Fig. 5(b). The engineered domain has switched to [01] and [10] and most of the defect dipoles are aligned with domain directions as seen in 13

the enlarged area in Fig. 5(c), so it is inferred that those defect dipoles has experienced an 180o switching at this frequency. In Fig. 5(d), the frequency is 100 Hz and the domain state is quite different from that in Fig. 5(b). Through the enlarged area in Fig. 5(e), one can see that the direction of polarization in each domain is no longer unidirectional. The local polarization deviation is caused by the electric field generated by the un-switched defect dipoles. In fact, most defect dipoles maintain the initial setting because they could not follow up with the field change. So, the ferroelectric system has been completely switched at 0.1 Hz the switching could not be completed at 100 Hz, the unswitched defect dipoles produced an internal bias field causing the hysteresis loop become asymmetric.

Fig. 5. Domain configurations at different frequencies of electric field, (a) initial domain state, (b) 14

domain state at the largest electric field along [11] at 0.1 Hz, (c) domain state at the largest electric field along [11] direction at 100 Hz, (d) an enlarged area in (b), and (e) an enlarged area in (c).

Besides the frequency of electric field, another important factor affecting the shape of the de-aging process is the temperature [41, 42]. P-E loops at 3 different temperatures were simulated at 100 Hz as shown in Fig. 6. As can be seen from Fig. 6, the effect of lowering temperature is equivalent to increasing frequency of electric field. As the temperature increases from 298 K to 378 K, the coercive fields gradually become smaller and the internal bias field decrease from 0.24 to 0.19 in normalized scale. So the de-aging process is easier to achieve at higher temperatures.

Fig. 6. The P-E loops at 100 Hz at different temperatures for the acceptor-doped PZT system by phase field simulation.

This phenomenon is caused by the increasing mobility of the oxygen vacancies at higher temperatures [1, 43]. According to the transition rate Eq. (10), the “speed” of defect dipole rotation would increase drastically with increasing temperature, which means that under the same field and a given time period, more defect dipoles could be rotated to the new orientation if the temperature is higher. Fig. 7 shows the time 15

evolution of defect dipole percentage of different directions at different temperatures. At room temperature (298 K), most defect dipoles could not follow the electric field change. As the temperature increased to 318 K, 338 K and 358 K, more defect dipoles could catch up with the electric field, so the oscillation of the curves become more intense, as shown in Fig. 7(b), 7(c) and 7(d).

Fig. 7. The time evolution of defect dipole percentage of different directions under 100 Hz AC field at different temperatures (a) 298 K, (b) 318 K, (c) 338 K, and (d) 358 K.

The domain configurations in the de-aging process are simulated at different temperatures and different frequencies in Fig. 8, where the instantaneous electric field is along [11]. At 0.1 Hz, the domain configurations at four temperatures are roughly the same, indicating that at this frequency, all defect dipoles are aligned with the spontaneous polarization even when the electric field has been changed. As the frequency increases, the switching of the defect dipole gradually fails to keep up with 16

the electric field at low temperatures, but this effect could be countered by increasing temperature. For example, at 10 Hz, there are still a lot of un-aligned defect dipoles at 338 K, while only a few defect dipoles are un-aligned at 358 K, indicating that more complete de-aging occurred at 358 K. Hence, both low frequency and high temperature will help the de-aging process.

Fig. 8. Domain configuration for the [11]-poled acceptor-doped ferroelectric system at various temperatures and frequencies under AC electric field.

It is known that the internal bias field is an indicator to the de-aging degree. If more defect dipoles point to one direction, the absolute value of the internal field would become larger. The relationship between the calculated normalized internal bias field and frequency at three temperatures are shown in Fig. 9. At 298 K, the normalized internal bias field increases sharply from 0.04 to 0.22 in normalized scale as the frequency changing from 0.1 to 10 Hz and then grows steadily from 0.22 to 0.24 as the frequency changes from 10 to 100 Hz. As the temperature increases, the 17

value of the internal bias field gradually becomes smaller, which means a more symmetrical P-E loop. It can be seen from Fig. 9 that the internal bias field does not increase significantly when the electric field frequency is larger than 10 Hz at various temperatures, indicating that the switching behavior of defect dipoles under the electric field frequency larger than 10 Hz is almost the same. The theoretical relation between the value of internal bias field and field frequency is similar to experiment results reported in relaxor-PbTiO3 system [1, 44].

Fig. 9. The relation between normalized internal bias field and frequency at various temperatures.

4. Conclusion In summary, a phase field model combined with kMC method is proposed to study the de-aging process of acceptor-doped ferroelectrics. The shifted P-E loops at different frequencies were simulated at room temperature (298 K). The P-E loops shift as the frequency increases and the shift amount soon saturates with frequency. The percentages of defect dipoles along different orientations were calculated at different 18

frequencies. At low frequency (0.1 Hz), switching of defect dipoles could catch up with the AC electric field, the P-E loop is symmetric. As the frequency increases, the oscillation of the percentage curve gradually becomes smaller. High temperature effect is similar to low frequency effect in terms of influencing the P-E loops due to the high mobility of oxygen vacancies at higher temperatures. Domain configurations before and after de-aging process were simulated at low frequency (0.1Hz) and high frequency (100 Hz) and we show that the directions of defect dipoles have a significant influence on the final domain state. Domain structures were simulated in the de-aging process at different temperatures and frequencies to study the effect of defect dipoles on the formation of domain structures. In addition, the relation between the normalized internal bias field and frequency at various temperatures was calculated. It was found that the normalized internal bias field gradually becomes larger as the frequency increases and saturate quickly beyond 100 Hz. Acknowledgments This research was supported by the NSFC under Grant No. 11304061, 51572056, 11572103, and 51502055.

19

References [1] X. Qi, E. Sun, S. Li, W. Lü, R. Zhang, B. Yang, W. Cao, Dynamic scaling of internal bias field in Mn-doped 0.24Pb(In1/2Nb1/2)O3–0.42Pb(Mg1/3Nb2/3)O3–0.34PbTiO3 ferroelectric ceramic, J. Mater. Sci., 53 (2018) 12762-12769. [2] P. Tung, J. Daniels, M. Major, D. Schneider, R. Chen, H. Luo, T. Granzow, Achieving large electric-field-induced strain in lead-free piezoelectrics, Mater. Res. Lett., 7 (2019) 173-179. [3] H. Qi, R. Zuo, X. Zhou, D. Zhang, Phase structure dependence of acceptor doping effects in (Bi0.5Na0.5)TiO3-BaTiO3 lead-free ceramics, J. Alloy. Compd., 802 (2019) 6-12. [4] L. Wang, Z. Zhou, X. Zhao, Z. Liu, R. Liang, X. Dong, Enhanced strain effect of aged acceptor-doped BaTiO3 ceramics with clamping domain structures, Appl. Phys. Lett., 110 (2017) 102904. [5] L. Zhang, X. Ren, Aging behavior in single-domain Mn-doped BaTiO3 crystals: implication for a unified microscopic explanation of ferroelectric aging, Phys. Rev. B, 73 (2006) 094121. [6] D. Xue, J. Gao, L. Zhang, H. Bao, W. Liu, C. Zhou, X. Ren, Aging effect in paraelectric state of ferroelectrics: implication for a microscopic explanation of ferroelectric deaging, Appl. Phys. Lett., 94 (2009) 082902. [7] D. Lee, B. Jeon, S. Baek, S. Yang, Y. Shin, T. Kim, Y. Kim, J. Yoon, C. Eom, T. Noh, Active control of ferroelectric switching using defect-dipole engineering, Adv. Mater., 24 (2012) 6490-6495. [8] E. Sapper, R. Dittmer, D. Damjanovic, E. Erdem, D. Keeble, W. Jo, T. Granzow, J. Rödel, Aging in the relaxor and ferroelectric state of Fe-doped (1-x)(Bi1/2Na1/2)TiO3-xBaTiO3 piezoelectric ceramics, J. Appl. Phys., 116 (2014) 104102. [9] C. Chen, Z. Zhou, Z. Liu, R. Liang, W. Zhang, X. Dong, Enhanced and stable strain memory in Mn-doped Pb(Mn1/3Sb2/3)O3-Pb(Zr,Ti)O3 ceramics realized by sesquipolar loading, J. Am. Ceram. Soc., 101 (2018) 1020-1024. [10] J. Glaum, Y. Genenko, H. Kungl, L. Schmitt, T. Granzow, De-aging of Fe-doped lead-zirconate-titanate ceramics by electric field cycling: 180°-vs. non-180° domain wall processes, J. Appl. Phys., 112 (2012) 034103. [11] Y. Gao, K. Uchino, D. Viehland, Effects of thermal and electrical histories on hard piezoelectrics: a comparison of internal dipolar fields and external dc bias, J. Appl. Phys., 101 (2007) 054109. [12] R. Eichel, P. Erhart, P. Traskelin, K. Albe, H. Kungl, M. Hoffmann, Defect-dipole formation in copper-doped PbTiO3 ferroelectrics, Phys. Rev. Lett., 100 (2008) 095504. [13] G. Du, R. Liang, L. Wang, K. Li, W. Zhang, G. Wang, X. Dong, Linear temperature scaling of ferroelectric hysteresis in Mn-doped Pb(Mn1/3Sb2/3)O3-Pb(Zr,Ti)O3 ceramic with internal bias field, Appl. Phys. Lett., 102 20

(2013) 142903. [14] C. Long, H. Fan, W. Ren, J. Zhao, Double polarization hysteresis and dramatic influence of small compositional variations on the electrical properties in Bi4Ti3O12 ceramics, J. Eur. Ceram. Soc, 39 (2019) 4103-4112. [15] X. Ren, Large electric-field-induced strain in ferroelectric crystals by point-defect-mediated reversible domain switching, Nat. Mater., 3 (2004) 91-94. [16] D. Lupascu, Y. Genenko, N. Balke, Aging in Ferroelectrics, J. Am. Ceram. Soc., 89 (2006) 224-229. [17] C. Guvenc, U. Adem, Influence of aging on electrocaloric effect in Li+ doped BaTiO3 ceramics, J. Alloy. Compd., 791 (2019) 674-680. [18] B. Akkopru‐Akgun, W. Zhu, M. Lanagan, S. Trolier-McKinstry, The effect of imprint on remanent piezoelectric properties and ferroelectric aging of PbZr0.52Ti0.48O3 thin films, J. Am. Ceram. Soc., 102 (2019) 5328-5341. [19] D. Xue, J. Gao, H. Bao, Y. Zhou, L. Zhang, X. Ren, In situ observation of thermally activated domain memory and polarization memory in an aged K+-doped (Ba, Sr)TiO3 single crystal, J. Phys. Condens. Mat., 23 (2011) 275902. [20] M. Zhu, L. Liu, Y. Hou, H. Wang, H. Yan, Microstructure and electrical properties of MnO-doped (Na0.5Bi0.5)0.92Ba0.08TiO3 lead-free piezoceramics, J. Am. Ceram. Soc., 90 (2007) 120-124. [21] C. Yang, E. Sun, B. Yang, W. Cao, Theoretical study on local domain pinning effect due to defect dipole alignment, J. Phys. D, 51 (2018) 415303. [22] Y. Zhang, J. Li, D. Fang, Oxygen-vacancy-induced memory effect and large recoverable strain in a barium titanate single crystal, Phys. Rev. B, 82 (2010) 064103. [23] Y. Ma, B. Xu, K. Albe, A. Grünebohm, Tailoring the electrocaloric effect by internal bias fields and field protocols, Phys. Rev. Appl., 10 (2018) 024048. [24] Y. Zhou, D. Xue, X. Ding, L. Zhang, J. Sun, X. Ren, Modeling the paraelectric aging effect in the acceptor doped perovskite ferroelectrics: role of oxygen vacancy, J. Phys. Condens. Mat., 25 (2013) 435901. [25] Y. Genenko, Space-charge mechanism of aging in ferroelectrics: an analytically solvable two-dimensional model, Phys. Rev. B, 78 (2008) 214103. [26] X. Cao, H. Tian, C. Hu, F. Huang, Y. Wang, X. Sun, Z. Zhou, Defect dipole evolution and its impact on the ferroelectric properties of Fe-doped KTN single crystals, J. Am. Ceram. Soc., 102 (2019) 3117-3122. [27] P. Erhart, P. Träskelin, K. Albe, Formation and switching of defect dipoles in acceptor-doped lead titanate: a kinetic model based on first-principles calculations, Phys. Rev. B, 88 (2013) 024107. [28] Y. Sheng, Y. Huang, C. Chen, M. Zhang, N. Deng, L. Ma, Effect of oriented defect-dipoles on the ferroelectric and piezoelectric properties of CuO-doped (K0.48Na0.52)0.96Li0.04Nb0.805Ta0.075Sb0.12O3 ceramics, Ceram. Int., 44 (2018) 10141-10146. [29] T. Wang, Y. Liao, D. Wang, Q. Zheng, J. Liao, F. Xie, W. Jie, D. Lin, Cyclingand heating-induced evolution of piezoelectric and ferroelectric properties of CuO-doped K0.5Na0.5NbO3 ceramic, J. Am. Ceram. Soc. , 102 (2019) 351-361. 21

[30] G. Geneste, J. Kiat, H. Yokota, Y. Uesu, Dielectric relaxation in Li-doped KTaO3 studied by kinetic Monte Carlo, Phys. Rev. B, 83 (2011) 184202. [31] C. Yang, E. Sun, B. Yang, W. Cao, Modeling dynamic rotation of defect dipoles and poling time dependence of piezoelectric effect in ferroelectrics, Appl. Phys. Lett., 114 (2019) 102902. [32] S. Nambu, D. Sagala, Domain formation and elastic long-range interaction in ferroelectric perovskites, Phys. Rev. B, 50 (1994) 5838. [33] P. Ondrejkovic, P. Marton, M. Guennou, N. Setter, J. Hlinka, Piezoelectric properties of twinned ferroelectric perovskites with head-to-head and tail-to-tail domain walls, Phys. Rev. B, 88 (2013) 024114. [34] R. Ahluwalia, T. Lookman, A. Saxena, W. Cao, Domain-size dependence of piezoelectric properties of ferroelectrics, Phys. Rev. B, 72 (2005) 014112. [35] D. Wang, X. Ke, Y. Wang, J. Gao, Y. Wang, L. Zhang, S. Yang, X. Ren, Phase diagram of polar states in doped ferroelectric systems, Phys. Rev. B, 86 (2012) 054120. [36] P. Marton, I. Rychetsky, J. Hlinka, Domain walls of ferroelectric BaTiO3 within the Ginzburg-Landau-Devonshire phenomenological model, Phys. Rev. B, 81 (2010) 144125. [37] H. Hu, L. Chen, Three-dimensional computer simulation of ferroelectric domain formation, J. Am. Ceram. Soc., 81 (1998) 492-500. [38] R. Ahluwalia, W. Cao, Influence of dipolar defects on switching behavior in ferroelectrics, Phys. Rev. B, 63 (2000) 012103. [39] S.M. Yang, J.Y. Jo, T.H. Kim, J.G. Yoon, T.K. Song, H.N. Lee, Z. Marton, S. Park, Y. Jo, T.W. Noh, Ac dynamics of ferroelectric domains from an investigation of the frequency dependence of hysteresis loops, Phys. Rev. B, 82 (2010) 174125. [40] T. Rojac, S. Drnovsek, A. Bencan, B. Malic, D. Damjanovic, Role of charged defects on the electrical and electromechanical properties of rhombohedral Pb(Zr,Ti)O3 with oxygen octahedra tilts, Phys. Rev. B, 93 (2016) 014102. [41] S. Upadhyay, V. Reddy, K. Sharma, A. Gome, A. Gupta, Study of aging & de-aging behaviour in un-doped polycrystalline ferroelectric BaTiO3, Ferroelectrics, 437 (2012) 171-180. [42] S. Yun, X. Wang, J. Shi, J. Zhu, D. Xu, Ferroelectric properties of barium calcium titanate ceramics doped with bismuth oxide, Mater. Lett., 63 (2009) 1883-1885. [43] K. Yan, F. Wang, D. Wu, X. Ren, K. Zhu, Ferroelectric aging effects and large recoverable electrostrain in ceria-doped BaTiO3 ceramics, J. Am. Ceram. Soc., 102 (2018) 2611-2618. [44] S. Upadhyay, V. Reddy, Study of ferroelectric hysteresis scaling exponents in aged polycrystalline BaTiO3, Ferroelectrics, 445 (2013) 147-151.

22

A phenomenological de-aging model for acceptor-doped ferroelectrics is proposed. The dynamic switching process of defect dipoles was simulated by Monte Carlo method. The number of defect dipoles along different directions is a key factor affecting the P-E loop. The internal bias field gradually becomes larger as the frequency increases and saturate beyond 100 Hz.