Antiferroelectric polarization switching and dynamic scaling of energy storage: A Monte Carlo simulation
B. Y.Huang,1Z. X.Lu,1,2Y.Zhang,1,a)Y. L.Xie,1M.Zeng,3Z. B.Yan,1and J.-M.Liu1,3
1Laboratory of Solid State Microstructures and Innovation Center of Advanced Microstructures, Nanjing University, Nanjing 210093, China
2College of Mechanical Engineering, Linyi University, Linyi 276000, China
3Institute for Advanced Materials and Laboratory of Quantum Engineering and Materials, South China Normal University, Guangzhou 510006, China
(Received 14 January 2016; accepted 20 April 2016; published online 5 May 2016)
The polarization-electric field hysteresis loops and the dynamics of polarization switching in a two- dimensional antiferroelectric (AFE) lattice submitted to a time-oscillating electric field E(t) of frequency f and amplitude E0, is investigated using Monte Carlo simulation based on the Landau–Devonshire phenomenological theory on antiferroelectrics. It is revealed that the AFE double-loop hysteresis areaA, i.e., the energy loss in one cycle of polarization switching, exhibits the single-peak frequency dispersionA(f), suggesting the unique characteristic time for polarization switching, which is independent ofE0as long asE0is larger than the quasi-static coercive field for the antiferroelectric–ferroelectric transitions. However, the dependence of recoverable stored energyW on amplitudeE0seems to be complicated depending on temperatureTand frequencyf.
A dynamic scaling behavior of the energy loss dispersionA(f) over a wide range ofE0is obtained, confirming the unique characteristic time for polarization switching of an AFE lattice. The present simulation may shed light on the dynamics of energy storage and release in AFE thin films.
Published by AIP Publishing.[http://dx.doi.org/10.1063/1.4948476]
I. INTRODUCTION
It was not until 1951 when Kittel investigated and put for- ward related phenomenological theory on antiferroelectric (AFE) materials as a subclass of ferroelectric materials.1In an AFE lattice, unlike a ferroelectric (FE) one, adjacent electric dipoles are oriented in opposite directions. Thus, its total spon- taneous polarizationPis zero under zero electric field, since the adjacent dipoles cancel out each other. Accompanying with the double-loop hysteresis of electric polarization P against electric fieldE, an AFE system shows quite different characteristics from an FE system. The phase transitions, do- main structure, and polarization switching of an AFE lattice, exhibiting some analogies with those of an FE lattice, also depend strongly on temperatureT, electric fieldE, microstruc- tural details, and so on.2,3 It is noted that these issues have been much less concerned due to the limited practical applica- tions of AFE materials in comparison with FE materials.
Nevertheless, one exception is the application potentials of AFE materials in high-energy-storage capacitors, which was suggested as early as 1961 by Jaffe and has been contin- uously investigated mainly from experimental aspects.4–8 When an electric fieldEis applied to an AFE material in a capacitor structure, the induced AFE-FE phase transitions take place, which thus stores the input energy. After the field is removed, the reversal FE-AFE transitions occur rapidly,9 leaving the input energy to be stored inside the capacitor at the open-circuit state and allowing its rapid release upon the connection to external loads. The advantage of AFE energy
storage and release lies in the fact that an AFE material usu- ally has relatively high dielectric constant but also high dielectric breakdown field on one hand, and on the other hand the field-induced double-loop hysteresis has a very small remnant polarization, giving rise to high energy release efficiency. In other words, except for a small heat loss due to the polarization switching, most of the input energy can be released in this cycle. Additionally, AFE materials possess good fatigue endurance.10,11
So far more than 100 AFE materials have been synthe- sized, among which PbZrO3 (PZO) represents one of the most popular ones.3,12–14 In spite of the attractive high- energy-storage application potentials of these materials, sub- stantial dependences of the energy storage performance on the AFE-FE-AFE transition cycle and microstructural details have been repeatedly investigated,6,11,15–17 raising quite a few of issues on the complexity in optimizing the energy storage performance. The three most concerned issues are:
(1) the stored energy densityW0; (2) the energy release effi- ciency gwhich is defined as the ratio of releasable energy densityW overW0(g¼W/W0) in one hysteresis cycle; and (3) the thermal and frequency stability and the energy stor- age/release cycling lifetime. Physically, the three issues are relevant with the so-called dynamic hysteresis in an AFE system. This dynamic hysteresis deals with the polarization switching behaviors in response to temperature T and time- varying electric field E(t) by generating a hysteresis from which the two parametersWandW0can be evaluated, where A¼W0W is the hysteresis loop area A¼Ð
PdE in one cycle and efficiencyg¼1A/W0. In particular and as a con- ventional scenario, the dynamic hysteresis treats the fre- quency dispersion A(f) and amplitude dependenceA(E0), or
a)Author to whom correspondence should be addressed. Electronic mail:
0021-8979/2016/119(17)/174103/9/$30.00 119, 174103-1 Published by AIP Publishing.
simplyA(f,E0), at a givenT, where parametersfandE0are defined by time-varying electric fieldE(t)¼E0sin(xt) with x¼2pf and t the time. This physics can be schematically illustrated in Fig.1where a model double-loop hysteresis for a typical AFE material is sketched.
Surely, a thorough understanding of the dynamic hystere- sis for an AFE lattice and associated polarization switching becomes imperative for efficient energy-storage device design. Our target is to attain the largest stored energy density W and highest energy release efficiency g simultaneously.
At the same time, the thermal stability (i.e.,T-dependence of A(f,E0)) is also important, which will be discussed elsewhere.
However, very differently, the dynamic hysteresis for conven- tional FE materials has been extensively investigated since 1990s.18–20 These studies revealed the intrinsic relationship between the dynamic hysteresis and FE domain switching.
For an FE lattice, the domain switching can be kinetically described by a mono-dispersive domain size distribution and domain wall motion speed, resulting in the well-known sin- gle-peak frequency dispersion A(f), while the domain wall motion speed is dependent of field amplitude. In details, for a time-varyingE(t), the frequency dispersionA(f) of an FE lat- tice can be scaled by
A fð Þ ¼A0þEa0fbg f
Ec0 ; (1)
whereA0is the area in the f )0 limit counting the effect from non-dynamic origins,a,b, and care the scaling expo- nents, andgis a non-monotonic function which meetsg(x) )0 asx)0 or1.19,21This scaling reflects the nature of a unique characteristic timesfor domain switching at a given E0. In the simplest cases, this scaling behavior can be approximately described in the two extreme ends of fre- quencyf18,22
Aðf;E0Þ /E2=30 f1=3 asf )0 Aðf;E0Þ /E20f1 as f ) 1;
(
(2) noting that the scaling exponents toE0andfmay be different from the stated values here, as discussed for various theoreti- cal models19,21,23and experiments on different materials.24–26
It is clear that the scaling behavior grasps the core nature of domain switching for an FE system.
The dynamic hysteresis and scaling behavior discussed earlier may not be applicable directly to an AFE lattice. For two stubbornly coupled antiparallel FE sublattices, no clear domain concept can be defined in the AFE state unless one sublattice is at least partially switched by electric field. The AFE-FE switching sequence may be kinetically scaled by one characteristic time (sc1) too. This switching sequence completes the AFE-FE transitions. Subsequently, the FE domains in the whole lattice may experience growth and coarsening, which is scaled by another characteristic time (sc2). Therefore, it seems that the dynamic hysteresis for an AFE lattice may not be scalable by one characteristic time but two times if any.
In fact, some experiments on the scaling behavior of AFE systems were reported. Kim and Kim proposed two scaling relations of the AFE hysteresis areaA(f,E0) in a mixed AFE betaine phosphate-arsenate (BP0.9A0.1) crystal at a frequency below 200 Hz.27 Chenet al. established the scaling behavior of energy storage densityWfor the saturated double-loop hys- teresis of Pb0.99Nb0.02[(Zr0.60Sn0.40)0.95Ti0.05]O3.28 However, these experiments covered relatively narrowf-range, which is insufficient for evaluating the dynamic hysteresis. More ex- perimental data covering broaderfandE0ranges are needed, which will be helpful for not only checking theoretical model if any but also revealing the underlying physical scenario for the dynamic hysteresis. On the other hand, there have been only a few simulation reports on the polarization switching of AFE materials.29,30In addition, no simulation on the scaling studies of AFE materials has been available to the best of our knowledge. We are only aware of the work of Misirlioglu et al.which deals with a two-dimensional Ising model to sim- ulate the FE-AFE hysteresis.31
In this work, we address the polarization switching and its scaling behavior in a model AFE lattice. For ferroelectrics, one may adopt phase-field model to obtain the domain struc- tures and to track its evolution.32However, it is rather difficult to simulate the situation of AFE.33 As discussed earlier, we can simply divide an AFE lattice into two coupled FE sublattices and consider a coupling energy between them to favor anti-parallel dipole alignment. We start from the Landau–Devonshire phenomenological theory on a tetragonal AFE lattice and present domain structures at different states associated with a double-loop P-E hysteresis. Subsequently, we pay attention to the energy loss, i.e., hysteresis area A(f, E0)¼Ð
PdE, and the recoverable energy W¼Ð EdP (upon an electric discharging), which are both sketched in Fig.1. It will be shown that the frequency dispersion can still be scaled by one characteristic time instead of two times if field amplitude E0 is sufficiently large over the quasi-static field for AFE-FE transition, EAF1, as shown in Fig.1. Based on the simulation results, one may derive some conclusion on the energy storage/release as a function offandE0.
Prior to the model description and results presentation, the difficulty in experimental checking of our simulation results should be mentioned. As subsequently shown, in spite of some qualitative consistencies between our proposed scal- ing behaviors and existing experimental data, to be discussed
FIG. 1. A sketch of energy storage for an AFE system as shown in a double- loopP-Ehysteresis, where the coercive fieldEAF1for the AFE-FE transition and fieldEAF2for the FE-AFE transition are marked.
later, a quantitative comparison becomes impossible since so far no any physical parameters in the framework of the Landau–Devonshire phenomenological theory on any AFE system are available. The available experimental data are insufficient for a convincing comparison. Therefore, our sim- ulation conclusion would appeal for further experimental verification, especially with respect to AFE domain structure and dynamic hysteresis.
II. MODEL AND SIMULATIONS A. Model description
We utilize the Monte Carlo simulation to track the dynamic hysteresis driven by field E(t) for a two- dimensionalLLAFE lattice.34–36The lattice consisting of two interpenetrated FE sublattices labeleda andb, respec- tively, is applied with the periodic boundary conditions.
Each site in sublattice a is surrounded with four nearest neighbor sites in sublatticeband four next nearest neighbor sites inaand vice versa. On each site are imposed two order parameters: electric dipole vectorP(r)¼(Px, Py), normalized byP0, and displacement vectoru(r)¼(ux, uy), normalized by lattice constanta0. This lattice is employed as an approach to an AFE thin film lattice where all electric dipoles are aligned in the x-y (a-b) plane. This approach seems to be over- simplified but our main motivation is to make a qualitative insight into the energy storage/release associated with an AFE thin film capacitor.
Based on the Landau phenomenological description of a tetragonal FE lattice,37 our simulation starts from the total free energy on this AFE lattice
Ftotal¼FabþX
a;b
FldþFgrþFelaþFele; (3) whereFabis an exclusive term associated with the exchange interactions between sublatticesaandb, whileFld,Fgr,Fela, and Fele are the Landau potential, gradient energy, elastic energy, and electric energy of sublattice a or b, respec- tively.37For simplicity,Fgr,Fela, andFeleconsider only the interaction between electric dipoles from the identical sublattice.
For accounting the AFE interaction, we follow the standard treatment of termFabwhich favors the anti-parallel alignment of the nearest dipole pairs1
Fab¼JX
hi;ji
PiPj; (4)
where hi,ji means a summation over all the nearest- neighboring dipole pairs and J>0 is the antiferroelectric coupling constant normalized by Landau coefficient ja0j.
Thus, all free energy terms are normalized by the energy unitja0jP20.
Similar to previous simulations on FE lattice,37for each FE sublattice, the Landau free energy polynomial function is expanded up to the sixth-order
Fld¼A1ðP2xþP2yÞ þA11ðP4xþP4yÞ þA111ðP6xþP6yÞ þA12ðP2xP2yÞ þA112ðP4xP2yþP4yP2xÞ; (5)
where the Landau coefficients are represented in coefficients A with various subscripts and A1¼A10(TT0), whereT0is the critical temperature for the paraelectric to FE transition of the corresponding FE sublattice.
The lowest-order expression of the gradient energy is Fgr ¼1
2
G11P2x;xþP2y;y
þG12Px;xPy;y
þG44ðPx;yPy;xÞ2þG044ðPx;yþPy;xÞ2
; (6)
where Pi,j¼@Pi/@rj, and coefficients Gij>0. Subsequently, termFelais defined as
Fela¼1
2C11u2x;xþu2y;y
þC12ux;xuy;yþ1 2C44u2x;y ðux;xgxxþuy;ygyyþux;ygxyÞ; (7) where ui,i¼@ui/@ri, ui,j¼@ui/@rjþ@uj/@ri with strain gxx
¼Q11Px2
þQ12Py2
,gyy¼Q11Py2
þQ12Px2
, andgxy¼Q44PxPy. Additionally, Cij are the elastic coefficients and Qij are the electrostrictive coefficients. It should be mentioned that the total elastic energy Felagiven in Eq. (7)looks quite different from that reported in literature.38,39 In fact, the total elastic energy Fela is usually written as Fela¼1/2Cijkl(eije0ij) (ekle0kl) whereCijklis the elastic stiffness tensor, andeijand e0ijare the total strains and stress-free strain, respectively, with the tensor subscripts (i,j,k,l¼1, 2, 3).39This energy term can be re-formulated into three parts. The first part is related only with polarizationPwhich can be merged into the Landau free energy term. The second part counts the conventional elastic energy and the third one counts the electrostriction energy, i.e., the second part of Eq. (7).39 While the details of this re- formulation are omitted here, Eq. (7) is actually equivalent with the original formulation of Fela for the latter two parts, while the first part is included into the Landau term Eq.(5). In fact, Eq.(7)has also been often employed.40
Finally, termFeleconsisting of dipole–dipole interaction and electrostatic energy can be written as
Fele¼ 1 4pe
X
hi;ji
PiPj
jrirjj33PiðrirjÞPiðrirjÞ jrirjj5
" #
ðExPxþEyPyÞ; (8) where the pre-factor (1/4pe)¼1 is taken for normalization purpose and fieldE(t) is applied along thex-axis. Instead of using the Ewald summation scheme for this long-range term, we adopt a more tractable calculation done by finite trunca- tion treatment. It is noted that for a two-dimensional lattice, this treatment is accurate as long as the truncating distance Ris big (R¼8 in our simulation).41,42
As mentioned earlier, one key problem with the present simulation is that so far neither phenomenological theory other than the above scheme on AFE lattice is available nor physical parameters for those AFE materials have been reported. As a qualitative approximation, all the physical coefficients except the AFE coupling constant J in Eq. (4) are chosen following the reference system BaTiO3, which
are listed in Table I.37 Unfortunately, the choice of all the coefficients regarding BaTiO3instead of any real AFE sys- tem makes a quantitative comparison with experiments impossible. ParameterJis treated as a variable in the simula- tion so that a representative double-loop hysteresis can be generated in a quasi-static condition at low T. In our case, it is shown thatJ¼1.00 is appropriate in terms of qualitative agreement with the major features of an AFE system. In fact, the most concerned AFE material Pb(Zr1–xTix)O3 with x0.5 has similar perovskite structure as BaTiO3.
B. Procedure of simulations
The MC simulation is performed in two steps. First, we simulate an AFE lattice using the gradual annealing method.
A lattice with randomly distributedPanduat a highTT0
is submitted to the thermal evolution toward the equilibrium, step by step down to a sufficiently lowT. In order to relax the electric dipoles and elastic strains in the lattices, we update the values ofPanduby the standard Metropolis pro- cedure, in which the trial dipole orientation is chosen at ran- dom. Considering that the response time of elastic strain is far shorter than that for dipole relaxation, we give ten chan- ces to its strain relaxation every time we relax one electric dipole.41,42 As a result, the initial dipole configuration of sublattice ais a ferroelectrically ordered state with all the dipoles aligned along thex-axis, while those from sublattice bare oriented in the opposite directions. The annealing pro- cess continues for 106mcs with 1 mcs standing forL2dipole flip attempts.
Second, the AFE lattice is then submitted to field E(t) below T0. The instant polarization P(t) averaged over the whole lattice configuration is plotted againstE(t) to produce the P-E hysteresis. For the high-f cases, the P-E hysteresis may be time-dependent in the initial periods of cycling and the presented hysteresis is different from the quasi-stable one in the very late stage. In general, the first 10 loops were dis- carded and the subsequent 20 loops were taken for the data averaging ofP(t). Finally, the loop areaA(f,E0) data package are evaluated. In our simulations, we pay our attention to the cases ofE0>EAF1, while the other cases will be addressed in the future due to the different physics.
III. RESULTS AND DISCUSSION A. Domain structures and evolution
We have performed extensive simulations on the evolu- tion of dipole configuration and dynamic hysteresis over broad ranges of variablesf, E0, and T. For each sublattice,
one can define the so-called FE domains, and such a domain refers to a region in which most dipoles align along one direction although minor dipoles may align differently.
We first look at the domain structure in various states (A, B, C, and D) of a typical double-loop hysteresis, as shown in Fig.2(a)wheref¼4.88104mcs1andE0¼8.0 atT¼0.2. The domain patterns of the two sublatticesaand b are plotted in Fig. 2(b). The domain orientations are marked by colors and red arrows, and anglehis defined with respect to the x-axis. Two major features deserve mention here. First, at state A, adjacent dipoles from different sublat- tices (aandb) are oppositely aligned, constituting an AFE state. However, each sublattice shows its own FE domain structure. For example, at state A, sublatticeadoes not ex- hibit a mono-domain but multi-domain structure, similarly does sublattice b. Second, at state B, all dipoles antiparallel with external field are completely switched, and an FE state over the whole lattice is reached. State C allows some local dipole fluctuations in the two sublattices, i.e., small cluster- like FE domains of opposite polarization embedded in the FE matrix. This evolution becomes more remarkable and those cluster FE domains grow and coarsen in size at state D in comparison with them at state C.
The above highlighted evolution of the AFE configura- tion in association with the double-loop hysteresis seems to hint two characteristic sequences. One is the dipole reversal of the two sublattices along the path from state A to state B.
The other is the appearance of cluster-like domains in the oppositely aligned matrix along the path from state B to states C and D. Indeed, the FE domain patterns and sizes of the two sublattices following the two paths are quite differ- ent. However, the two sequences are actually miscible, sug- gesting that the dipole switching sequences in the two sublattices are coherent and synchronous. It is thus implied that the polarization switching over the whole lattice may be possibly described by one characteristic time. This conjec- ture will be confirmed by the hysteresis dispersion and scal- ing analysis.
B. Dynamic hysteresis and energy storage
Now we investigate the dynamic hysteresis along two lines. One is the frequency dispersion A(f) at given E0and the other is theE0-dependence at givenf. The representative examples at T¼0.2 are plotted in Figs. 3(a)–3(d). For the cases of E0(5.0) only slightly larger thanEAF1, as shown in Fig.3(a), the double-loop feature can be observed only in the low-f range and the hysteresis evolves into the inclined fat single loop shape with increasing f. The loop area A(f)
TABLE I. Physical coefficients chosen for simulation (1/4pe¼1, Boltzmann constantkB¼1).
Parameter (unit) Value Parameter (unit) Value Parameter (unit) Value
A10(a10) 0.10 A11(a11P02/ja0j) 0.24 A12(a12P02/ja0j) 4.50 A111(a111P04/ja0j) 0.49 A112(a112P04/ja0j) 1.20 G11(g11/a20ja0j) 1.60 G12(g12/a20ja0j) 0.00 G44(g44/a20ja0j) 0.80 G’44(g044/a20ja0j) 0.80 C11(c11/ja0jP20) 2.75 C12(c12/ja0jP02) 1.79 C44(c44/ja0jP02) 0.54
Q11(q11/ja0j) 0.142 Q12(q12/ja0j) 0.0074 Q44(q44/ja0j) 0.0157
L(a0) 64 T0(Tc/(TcT1)) 5.00 J(ja0j) 1.00
first increases and then decreases. For the high-E0 cases (E0¼8.0>EAF1) shown in Fig. 3(b), qualitatively similar features to the low-E0 cases can be identified while the double-loop feature may sustain until a much higherf0. On the other hand, the E0-dependence exhibits different
characters, as shown in Figs.3(c)and3(d), respectively, for the low-f and high-f cases (f¼2.4107mcs1 and 2.4 104mcs1). It is seen that the double-loop feature remains nearly unchanged with increasing E0 in the low-f cases. However, the hysteresis can become fat upon increas- ing E0in the high-fcases, where the double-loop shape dis- appears at large E0, suggesting more complicated evolution of the domain structure in the high-frange.
The dynamic hysteresis discussed earlier suggests the significant dependences of the energy storage/release on f, E0, and T. To evaluate these dependences, we present the releasable energyWas a function ofE0at different frequen- cies in Fig.4(a) forT¼0.2 and Fig.4(b)for T¼0.01. At a relatively highT(T¼0.2),Wis enhanced with increasingE0
in the low-frange and then tends to be saturated in the large E0 range. However, a monotonous decreasing of W with increasingE0is identified in the high-frange due to the grad- ual disappearance of the double-loop feature. The hysteresis becomes expanded and fat with increasing E0, leading to reducedW. These results are qualitatively consistent with the experimental observations.28,43,44
Physically, it is reasonable to argue that a larger E0
always brings about largerdE/dtand shorter time for domain
FIG. 2. (a) A simulatedP-Ehysteresis atf¼4.88104mcs1andT¼0.2, givenE0¼8. Various states in the hysteresis are labeled with A, B, C, and D, respectively. (b) The snapshots of the domain structures for sublatticea (left column) and sublatticeb(right column) at states A, B, C, and D. The domain orientations are marked by colors and red arrows, and anglehis defined with respect to thex-axis. Dark yellow represents the polarization along thex-axis, and blue stands for the polarization along thex-axis. The light blue areas between the dark yellow and blue domains mark the domain walls. The magnitude of electric dipole vectors in white areas is zero.
FIG. 3. SimulatedP-Eloops for the AFE lattice atT¼0.2: (a)E0¼5 at var- iousf, (b)E0¼8 at variousf, (c)f¼2.4107mcs1with differentE0, and (d)f¼2.4104mcs1with differentE0.
FIG. 4. The recoverable energy densityW(E0) calculated from theP-Eloops simulated for the low-fand high-fcases with (a)T¼0.2 and (b)T¼0.01. The value offfor each curve is labeled numerically with the unit 106mcs1.
switching during field-driven transitions at any fixedf, result- ing in more remarkable hysteresis effect. Given a relatively highTand in the low-fcases where characteristic relaxation timesf1, the enhancedW withE0in the low-E0range is mainly ascribed to the field-driven AFE-FE transitions, which lead to much larger polarization P, but the double- loop hysteresis remains roughly unchanged in shape. In the high-f cases where sf 1, the decayed W withE0comes from the increased remnant polarization due to the intensi- fied hysteresis effect with gradual disappearance of the double-loop feature, as shown in Fig. 3(d). If temperature Tis reduced down so thatsf1becomes available even at extremely low f(107mcs1), it is inevitable that energy Was a function ofE0always falls down in thef-range cov- ered here above 107mcs1. This is the reason for the observations shown in Fig.4(b).
C. Hysteresis dispersion scaling
To evaluate the characteristic relaxation time as a func- tion ofE0at differentT, we plot the evaluated dispersionA(f) for variousE0given T¼0.1 and 0.4 in Figs.5(a) and5(b), respectively, where frequencyfcovering from 107mcs1to 101mcs1 is presented in the logarithmic scale. It is noted that the smallestE0is still larger thanEAF1. Referring to ear- lier works on the frequency dispersion on FE lattice, several points deserve mention here. First, similar to the FE lattice, all the dispersion curves at the two temperatures exhibit the single-peak pattern, with the peak coordinate designated as (fm, Am), suggesting a resonant state of the polarization switching at ffm. In other words, the characteristic time sfor polarization switching is sfm1. Second, it is found that peaked area Am increases with E0, but the peaked fre- quency fm does not show remarkable shifting toward the high-fside withE0, very different from the FE lattice.20The identifiableE0-dependence offm is very weak indeed. Third, the T-dependences of both fm and Am are also very weak
although a slight shifting of fm toward the high-f side with increasingTcan be detected.
Given the single-peak dispersion behavior, one can rescale these dispersions at various E0following the simple exponential scaling relationships:
A0ðf0Þ ¼AðfÞ=Ea0
f0¼f=Ec0; (9) whereaandcare the scaling exponents that can be obtained from the best fitting of the simulated data so that the whole set of dispersion curves at a given Tcan fall onto a master curve.
The best fitted results give the scaling exponents a¼1.0060.05 andc¼0.0560.05 forT<0.3, anda¼1.05 60.05 andc¼0.0560.05 forT>0.3 in the present simula- tions. The rescaled dispersion curves forT¼0.1 and 0.4 are plotted in Figs. 5(c) and 5(d), respectively, demonstrating that these curves do fall onto one master curve within the simulation uncertainties. The scaling of the data implies at least three facts: (1) the frequency dispersion A(f) can be indeed scaled by simple power-law transforms based on a single parameter E0 over a broad range of E0; (2) time s shows very weak dependence of E0; and (3) time s shows very weak dependence ofTas long asTT0.
To understand the reason for such weakE0-dependence of characteristic times, we come back to the hysteresis itself and correlate it with the kinetics of polarization switching in the characteristic time language. As discussed earlier, one refers to Fig.2(a)and starts from state A where the AFE order is the ground state. Given fixedf,E0, andT, fieldE(t) drives the evolution of the lattice from state A to state B via the AFE-FE transition. The corresponding characteristic time for such a transition is denoted as sAFE-FE. Subsequently, the dynamic hysteresis proceeds from state B to states C and D in sequence via the FE-AFE transition. The corresponding characteristic time is denoted assFE-AFE. Clearly, one expects that a well-expanded and fat hysteresis will be developed if sFE-AFEsAFE-FE; otherwise, a slim and thin hysteresis is in- evitable ifsFE-AFEsAFE-FE.
For comparing the time scales sFE-AFE and sAFE-FE, one discusses the single dipole switching sequence as a zero-order approximation of the polarization switching, considering the domain growth is much faster than the do- main nucleation for typical FE systems. The dipole switching for an FE lattice and an AFE lattice is shown in Fig. 6(a).
We consider the first quadrant of the hysteresis. For the FE lattice where all dipoles spontaneously align rightward and the AFE lattice where all dipoles are driven to align right- ward, a decreasing electric field E (its positive direction aligns rightward) may trigger one rightward dipole (cycled) to flip leftward (red arrow). The energy change associated with such a flip is
DH¼DHFE2pE; for FE lattice
DH¼DHAFE2pE; for AFE lattice; (10) whereDHFE>0 andDHAFE<0 are, respectively, the changes of interactions between this dipole and its neighbors upon
FIG. 5. Simulated frequency dispersionsA(f) at variousE0for (a)T¼0.1 and (b)T¼0.4, and re-scaled frequency dispersionsA0(f0) for (c)T¼0.1 and (d)T¼0.4. It is seen that all theA0(f0) data at differentE0collapse onto one single curve for both cases.
such a flip for the FE lattice and AFE lattice. The characteris- tic dipole switching timescan be roughly expressed as
s¼s0ðTÞexpðDH=kBTÞ; (11) wherekBis the Boltzmann constant ands0is the characteris- tic flip time for a free dipole at a givenT. For the AFE lat- tice, it is clear that the relationsFE-AFEsAFE-FE is always true, qualitatively illustrating why the peak frequencyfmis only weakly E0-dependent although the peak height Am is stronglyE0-dependent. The weakT-dependence offmcan be understood in a similar way.
For a better illustration of the above scenario, we present in Fig.6(b)theP-Ehysteresis loops corresponding tofm6.25 102mcs1at five different values ofE0>EAF1. Here, the electric field is re-scaled byE0and it is seen that the five loops are roughly overlapped in spite of delicate differences among them. While the delicate differences reflect the weak E0- dependence, the overlapping reflects the E0-independence.
This property is favored for AFE storage energy applications and one need not take care of the frequency sensitivity.
D. Energy storage density
Finally, we discuss the energy storage density W in response to variations off,E0, andTfor the present system.
The simulated results are summarized in Figs.7(a)and7(b) for two differentT. Based on the dynamic hysteresis analy- sis, it is obvious that the storage densityW always decays monotonously with increasing f over the whole f-range.
Nevertheless, such decay shows distinct features in different f-subranges. The rapid decay occurs in the intermediate f-subrange, while a plateau is observed in the low-fsubrange and the storage density becomes nearly zero in the high-f subrange. An unexpected fact is that the frequency threshold beyond which the decay becomes rapid is lower when E0is larger, suggesting that the frequency stability of the energy storage is worse at largerE0.
The continuous decay ofWoutside the high-fsubrange can be well understood considering the relationship W(f)
¼W0(f)A(f), while the total energy densityW0(f) is rela- tively insensitive to f as f<fm. As f>fm, the hysteresis shrinks along the polarization axis rapidly, leading to serious suppression ofW0and thus less and less energy can be stored and released. Besides, given the scaling behavior ofA(f,E0), it would be interesting to check the possible scaling over W(f) at differentE0. Similarly, we propose
W0ðf0Þ ¼WðfÞ=Em0
f0¼f=En0; (12) where exponentsmandncorrespond to exponentsaandcin Eq.(9)on scaling of dispersion A(f). The best fitting of all the W(f) data gives rise tom¼0.0560.05 andn¼ 1.0060.05 for T¼0.1, and m¼0.2060.05 and n¼ 0.9060.05 for T¼0.4. The re-scaled data are plotted in Figs.7(c)and 7(d) for the two temperatures and it is seen that the scaling on data atT¼0.1 is slightly better than that atT¼0.4. These scaling behaviors confirm again the single characteristic time for polarization switching in the present AFE system.
What should be noted is that the scaling behaviors as revealed in the present model rely on the assumption that the exchange interaction is completely attributed to the adjacent dipoles from different sublattices. Surely, for realistic AFE materials, the antiferroelectric coupling interactions mecha- nism is much more complicated, rendering a simulation tougher. Besides, we only consider the 2D thin film lattice with in-plane AFE domains here, which also deserves
FIG. 6. (a) A proposedP-switching model in an FE lattice and an AFE lat- tice, respectively. Each arrow represents an electric dipole. In the AFE lat- tice, blue and green arrows represent different sublattices. (b) Simulated P-E/E0loops for the AFE lattice atT¼0.2 andf¼6.25102mcs1with differentE0.
FIG. 7. Simulated energy storage densityWas a functionfat variousE0for (a) T¼0.1 and (b)T¼0.4, and re-scaled W0(f0) data for (c)T¼0.1 and (d)T¼0.4.
caution for guiding practical applications. Even though the proposed model is far from perfect and the basic assumption needs further verification, the current simulations qualita- tively agree well with some experiments and are expected to serve as a guide for predictive modeling of enhanced energy storage properties of AFE materials. The present simulation scenario can also be easily modified to study other AFE sys- tems. In this sense, demonstrating these scaling behaviors will allow further exploration of strategies to reliably boost the energy efficiency of materials in proper conditions.
Finally, as repeated several times already, we unfortu- nately have a brick-wall to be overcome in checking our simulated results with experimental data in quantitative sense. First, AFE materials as a kind of energy storage sup- porters have not been highly concerned until recently and thus no sufficient experimental data on the domain dynamics of AFE-FE inter-transitions driven by electric field have been available. This issue seems to be receiving attention recently. Second, the present simulations rely majorly on the availability of the free energy parameters on an AFE system.
Unfortunately, this availability is not true for any of so far concerned AFE materials. Third, a qualitative rather than quantitative comparison seems to be difficult either if not impossible due to the lack of experimental data of the dynamic hysteresis over broad ranges off, E0, and T. This difficulty, on the contrary, raises strong appealing for future experimental investigation of the dynamics of field-driven AFE-FE inter-transitions.
Even then, we would like to mention some recent experi- mental data for a crude comparison with our simulated results.
First, our simulations suggest that the energy storage density W is only weakly f-dependent in the low-f range, while the measured data on AFE Pb0.97La0.02(Zr0.95Ti0.05)O3thin films do show nearly f-independent energy density over f¼2 Hz to 10 kHz at room temperature.44 On the other hand, Chen et al. proposed a scaling relation Wf0.02(E0EAF1)0.08 on AFE Pb0.99Nb0.02[(Zr0.60Sn0.40)0.95Ti0.05]O3bulk ceramics, based on their measured data in the low-f range (1 Hz to 100 Hz) and low-E0 range (38 kV/cm to 56 kV/cm).28 The covered f-range and E0-range are narrow for checking our simulated data over broad f- and E0-ranges. However, from our simulated results, we do see the weakly enhanced energy density with increasingE0, as shown in Fig.4(a), which seems crudely consistent with this experiment.28
IV. CONCLUSION
In summary, we have simulated the electric field driven AFE-FE transitions and domain structural evolution in association with the double-loop hysteresis for a two- dimensional model AFE system, based on the Landau phe- nomenological theory by using Monte Carlo method. The dynamic hysteresis over broad ranges of parametersf,E0, and T for the cases of E0>EAF1 and TT0 have been carefully investigated. It is revealed that the energy storage densityWis strongly dependent onf,E0, andT, consistent with some experimental observations. Both the frequency dispersion A(f) and energy storage density W(f) can be scaled by simple power laws based on a single parameter
E0, which confirms that one characteristic timesis enough for the description of polarization switching over the whole lattice. The weak E0-dependences of characteristic time s and monotonous decaying behaviors of energy storage W in differentf-subranges are discussed and explained quali- tatively. The present work raises strong appealing for ex- perimental investigation of the dynamics of field-driven AFE-FE inter-transitions in AFE materials.
ACKNOWLEDGMENTS
This work was supported by the National 973 Projects of China (Grant No. 2015CB654602), the Natural Science Foundation of China (Grant Nos. 51431006 and 51301084), and the Jiangsu Provincial Natural Science Foundation of China (Grant No. BK20130576).
1C. Kittel,Phys. Rev.82, 729 (1951).
2G. Shirane, E. Sawaguchi, and Y. Takagi,Phys. Rev.84, 476 (1951).
3X. Hao, J. Zhai, L. B. Kong, and Z. Xu,Prog. Mater Sci.63, 1 (2014).
4B. Jaffe,Proc. IRE49, 1264 (1961).
5J. Parui and S. Krupanidhi,Appl. Phys. Lett.92, 192901 (2008).
6I. Burn and D. M. Smyth,J. Mater. Sci.7, 339 (1972).
7B. Xu, Y. Ye, Q. M. Wang, N. G. Pai, and L. E. Cross,J. Mater. Sci.35, 6027 (2000).
8X. Hao, J. Zhou, and S. An,J. Am. Ceram. Soc.94, 1647 (2011).
9W. Pan, C. Dam, Q. Zhang, and L. Cross,J. Appl. Phys.66, 6014 (1989).
10J. H. Jang, K. H. Yoon, and H. J. Shin, Appl. Phys. Lett. 73, 1823 (1998).
11X. Hao, J. Zhai, and X. Yao,J. Am. Ceram. Soc.92, 1133 (2009).
12H. Fujishita and S. Hoshino,J. Phys. Soc. Jpn.53, 226 (1984).
13D. J. Singh,Phys. Rev. B52, 12559 (1995).
14X. Tan, C. Ma, J. Frederick, S. Beckman, and K. G. Webber, J. Am.
Ceram. Soc.94, 4091 (2011).
15J. F. Li, D. D. Viehland, T. Tani, C. D. E. Lakeman, and D. A. Payne, J. Appl. Phys.75, 442 (1994).
16J. Ge, D. Remiens, J. Costecalde, Y. Chen, X. Dong, and G. Wang,Appl.
Phys. Lett.103, 162903 (2013).
17X. Hao, Y. Wang, L. Zhang, L. Zhang, and S. An,Appl. Phys. Lett.102, 163903 (2013).
18M. Rao, H. R. Krishnamurthy, and R. Pandit,Phys. Rev. B42, 856 (1990).
19F. Zhong, J. Dong, and D. Xing,Phys. Rev. Lett.80, 1118 (1998).
20J. M. Liu, W. M. Wang, Z. G. Liu, H. L. Chan, and C. L. Choy,Appl.
Phys. A75, 507 (2002).
21G. F. Mazenko and M. Zannetti,Phys. Rev. B32, 4565 (1985).
22M. Rao and R. Pandit,Phys. Rev. B43, 3373 (1991).
23J. M. Liu and Z. G. Liu,Appl. Phys. A70, 113 (2000).
24J. F. Scott, F. M. Ross, C. A. P. De Araujo, M. C. Scott, and M. Huffman, MRS Bull.21, 33 (1996).
25J. M. Liu, H. P. Li, C. K. Ong, and L. C. Lim,J. Appl. Phys.86, 5198 (1999).
26B. Pan, H. Yu, D. Wu, X. H. Zhou, and J. M. Liu,Appl. Phys. Lett.83, 1406 (2003).
27Y.-H. Kim and J.-J. Kim,Phys. Rev. B55, R11933 (1997).
28X. Chen, F. Cao, H. Zhang, G. Yu, G. Wang, X. Dong, Y. Gu, H. He, and Y. Liu,J. Am. Ceram. Soc.95, 1163 (2012).
29O. Maksimova, T. Petrova, and A. Maksimov, Bull. Russ. Acad. Sci., Phys.77, 1068 (2013).
30J. Fornier and B. Verweire,Ferroelectrics213, 159 (1998).
31I. Misirlioglu, L. Pintilie, K. Boldyreva, M. Alexe, and D. Hesse,Appl.
Phys. Lett.91, 202905 (2007).
32Y. Li, S. Hu, Z. Liu, and L. Chen,Appl. Phys. Lett.78, 3878 (2001).
33F. Xue, L. Liang, Y. Gu, I. Takeuchi, S. V. Kalinin, and L. Q. Chen,Appl.
Phys. Lett.106, 012903 (2015).
34Y. Wang, W. Zhong, and P. Zhang,Phys. Rev. B51, 5311 (1995).
35S. Li, J. Eastman, Z. Li, C. Foster, R. Newnham, and L. Cross,Phys. Lett.
A212, 341 (1996).
36J. Wang, M. Kamlah, T. Y. Zhang, Y. Li, and L. Q. Chen,Appl. Phys.
Lett.92, 162905 (2008).
37F. Xue, X. Gao, and J. M. Liu,J. Appl. Phys.106, 114103 (2009).
38L. Q. Chen,J. Am. Ceram. Soc.91, 1835 (2008).
39Y. Li, S. Hu, Z. Liu, and L. Chen,Acta Mater.50, 395 (2002).
40H. L. Hu and L. Q. Chen,Mater. Sci. Eng. A238, 182 (1997).
41B. Li, X. Liu, F. Fang, J. Zhu, and J. M. Liu,Phys. Rev. B73, 014107 (2006).
42X. Wang, J. M. Liu, H. L. Chan, and C. Choy,J. Appl. Phys.95, 4282 (2004).
43Y. Wang, Z. Lv, H. Xie, and J. Cao,Ceram. Int.40, 4323 (2014).
44C. Liu, S. Lin, M. Qin, X. Lu, X. Gao, M. Zeng, Q. Li, and J. M. Liu, Appl. Phys. Lett.108, 112903 (2016).
