Monte Carlo simulation of the dielectric susceptibility of Ginzburg-Landau mode relaxors

Monte Carlo simulation of the dielectric susceptibility of Ginzburg-Landau mode relaxors

J.-M. Liu*and X. Wang

Laboratory of Solid State Microstructures, Nanjing University, Nanjing 210093, China and International Center for Materials Physics, Chinese Academy of Sciences, Shenyang, China

H. L. W. Chan and C. L. Choy

Centre for Smart Materials, Hong Kong Polytechnic University, Kowloon, Hong Kong 共Received 11 August 2003; revised manuscript received 16 October 2003; published 24 March 2004兲 The electric dipole configuration and dielectric susceptibility of a Ginzburg-Landau model ferroelectric lattice with randomly distributed defects are simulated using the Monte Carlo method. The simulated charac- teristics of the lattice configuration and dielectric susceptibility indicate that the model lattice evolves from a normal ferroelectric state to a typical relaxor state with increasing defect concentration. Consequently, the energy and dielectric susceptibility characteristics associated with the ferroelectric phase transitions become smeared. The simulated results approve the applicability of the Ginzburg-Landau model in approaching relaxor ferroelectrics.

DOI: 10.1103/PhysRevB.69.094114 PACS number共s兲: 77.80.Dj, 77.84.Jd, 77.22.Gm

Ferroelectric relaxors have been receiving attention from physicists and materials researchers for over 40 years, mainly because of their fascinating electric-dipole ordering and phase transition phenomena and the excellent dielectric and electromechanic properties for practical applications.^{1– 6} The well-established features for the dielectric response of relaxors include the diffusive ferroelectric phase transition, strong frequency dispersion, and sensitivity to external elec- tric bias. Over the past 40 years, a number of theoretical models were proposed to explain these features that are quite different from normal ferroelectrics, and those well- documented models include the compositional inhomoge- neous model,¹spin共dipole兲-glass-like model,³superparaelec- tric model,² and defect model,^2,4 etc. One essence of these models is the existence of internal random field as induced by either compositional inhomogeneity or defects. This con- cept is popular and well accepted.⁵The major effect of the internal random field is characterized by the coexistence of nanosized electric dipole clusters embedded in the matrix of paraelectric phase. These nanoclusters hold their stability over a wide range of temperature T, leading to a weak hys- teresis of polarization above the phase transition point and the frozen microregions of electric dipoles at low T. In fact, the system at low T may exhibit the same behaviors as the normal ferroelectrics.

While some of the models mentioned above consider re- laxors as spinlike systems,⁷the Hamiltonian includes a term accounting the randomly distributed internal field. In addi- tion to the statistical mechanics approach, the Monte Carlo 共MC兲 method represents a popular technique employed to investigate the static and dynamic dielectric behaviors of relaxors.⁸ The predicted behaviors are similar somehow to those identified for spin-glass systems. The simplest and rep- resentative dynamics is the multistate Potts model with a random Potts field,^9,10in which the random field is imposed by assigning a variable spin-interaction factor or a variable internal electrostatic energy obeying some assumed distribu- tion function, such as the Gaussian distribution. In these models, the dielectric relaxation is dynamically modulated

by the randomly distributed field. Consequently, the dielec- tric behaviors under different internal and external conditions have been investigated.^11–13However, for an electric-dipole ordered lattice, the involved interactions may be more com- plicated. One needs to consider the Landau potential f_L, the dipole-dipole interaction f_dip, the gradient energy f_Gassoci- ated with the spatial distribution of dipoles and long-range elastic interaction f_E, in addition to the electric field induced electrostatic energy f_S.^{14 –16} It can be argued that in these interactions, the latter four terms are the resultant terms upon the dipole moment and alignment in the lattice, and only the Landau free energy is the origin to generate an electric di- pole. Therefore, one may argue that introduction of a defect into the lattice will mainly influence f_L, and consequently modulate the other free energy terms. In fact, this is the main argument of the Ginzburg-Landau 共GL兲 model for relaxors recently developed by Semonovskaya et al.,¹⁷ which allows us to access the evolution of the dipole configuration and dielectric property in a ferroelectric lattice as function of the induced defects.

In this report, we would like to perform a Monte Carlo study on the dielectric relaxation behaviors of a GL-model relaxor. The MC simulation is performed on a two- dimensional 共2D兲 L⫻L lattice with periodic boundary con- ditions applied. Four equivalent orientation states ⫾关0,1兴,

⫾关1,0兴 for each dipole are allowed. We note that a three- dimensional 共3D兲 simulation would produce more reliable simulations than the 2D simulation. In the earlier simulations based on the random-field models, 3D 16⫻16⫻16 lattice was often employed, however, the one-dimensional size 共16 lattice units兲is too small for the present simulation based on the GL model, since the dipole-cluster statistics seems to be quite bad if one notes that the intercluster separation is⬃10 lattice units. A 3D lattice larger than 40⫻40⫻40 makes the computational capability unavailable to us because we would simulate the lattice evolution over a wide temperature range.

As for the 2D simulations, Semonovskaya et al.¹⁷ employed a lattice as large as 128⫻128. In our presimulation, we did not find substantial difference when L is reduced to 64 in

terms of the system energy and dielectric susceptibility evaluation. Therefore, we take L⫽64 in our simulations.

The GL theory on normal ferroelectrics was proposed earlier.^15–17 For each lattice site, a dipole vector P

⫽关P_x(r), P_y(r)兴 with its moment and orientation defined by the system energy minimization, where P_x and P_y are the two components along the x axis and y axis, respectively, is imposed. It is known that for ferroelectric crystals the fer- roelastic effect cannot be ignored. However, this effect would be weak for relaxors where no long-range ordering structure exists. The Landau double-well potential f_Lis writ- ten as¹⁶

f_L共P_i兲⫽A₁共P_x²⫹P_y²兲⫹A₁₁共P_x⁴⫹P_y⁴兲⫹A₁₂P_x²P_y²

⫹A₁₁₁共P_x⁶⫹P_y⁶兲, 共1兲 where subscript i refers to lattice site i and A₁, A₁₁, A₁₂, and A₁₁₁are the energy coefficients, respectively. For normal ferroelectrics, the first-order phase transitions occur if A₁

⬍0. When a spatial distribution of the dipoles exists, the as-induced gradient energy is f_G:^16,17

f_G共P_{i, j}兲⫽1

2G₁₁共P_x,x² ⫹P_{y ,y}² 兲⫹G₁₂P_x,xP_{y ,y}

⫹1

2G₄₄共P_x,y⫹P_{y ,x}兲²⫹1

2G₄₄

⬘

共P_x,y⫺P_{y ,x}兲², 共2兲 where P_{i, j}⫽⳵^Pi/⳵^xj. Since parameters G₁₁, G₁₂, G₄₄, and G₄₄

⬘

are all positive, f_G⬎0 in general, which favors the ho- mogeneous dipole alignment in the lattice. The dipole-dipole interaction f_dipis long-ranged. In the SI unit, this term at site i is written as¹⁶

f_dip共P_i兲⫽ 1

8␲␧0␹

兺

_具j典

冋

^{P兩共rrii⫺兲Pr共j兩r3j兲}

⫺3关P共r_i兲共r_i⫺r_j兲兴关P共r_j兲共r_i⫺r_j兲兴

兩r_i⫺r_j兩⁵

册

^{, 共3兲}

where具j典represents a summation over all sites within a cycle of infinite radius (R⇒⬁) centered at site i, parameters r_i, r_j, P(r_i), and P(r_j) here should be vectors, r_iand r_jare the coordinates of sites i and j, respectively. However, in simu- lation, a finite cutoff is needed and we take R⫽8 in our simulation. A difference of ⬃2% is estimated as compared with the value obtained at R⫽30 because f_dip decays very rapidly with vector (r_i⫺r_j). Therefore, taking R⫽8 will not produce substantial error in the simulation. A minimizing of f_dipover the lattice favors a head-to-tail alignment of dipoles.

Finally, the external electric field E introduces the electro- static energy

f_S共P_i兲⫽⫺P_i•E, 共4兲 where P_i and E should be vectors. The total interaction en- ergy counting these interactions is

H⫽

兺

_具i典 ^{关fL⫹fG⫹fdip⫹fS兴, 共5兲}

which serves as the Hamiltonian for the system used in our simulation, where具i典represents a summation over the whole lattice.

As mentioned above, in the GL model,¹⁷introduction of a defect into the lattice is thought to change the stability of a local dipole. For example, in Pb(Mg_1/3Nb_2/3)O₃共PMN兲with acceptor dopants, where O^⫺2 vacancies are introduced. Each vacancy may be combined with a metal impurity ion to form a defect dipole which may be imposed onto the local lattice dipole. Therefore, the moment of dipoles at these sites at- tached with defects can be enhanced or suppressed. In other words, the stability of these sites for a local dipole becomes site dependent. Because the stability of a dipole共its moment兲 is mainly determined by coefficient A₁of the Landau poten- tial Eq. 共1兲, the randomly distributed defects in fact impose the spatial-dependent coefficients A₁, A₁₁, A₁₂, and A₁₁₁in Eq.共1兲. Following the argument of Semenovskaya et al., it is assumed that only A₁ is affected by the defects¹⁷

A₁共r_i兲⫽A₁₀⫹b_mc,

A₁₀⫽␣共T⫺T⁰兲, ␣⬎0, 共6兲 where␣⬎0 is a materials constant, T⁰ is a critical tempera- ture for a normal ferroelectric crystal of the first-order phase transitions, parameter c(⫽0,1) labels the defect state of a site. c⫽1 means that site i is imposed with a defect and it remains perfect if c⫽0. Parameter b_mis the coefficient char- acterizing the influence of defects on the dipole stability. In this model, coefficient b_mcan be positive or negative in our simulation in order to model a suppression or enhancement of the dipole stability.

The MC simulation is performed via the following proce- dure. For a lattice, each site is imposed a dipole with its moment P and orientation randomly chosen within 共0–1.0兲 and the four equivalent states. Also, on each site is probably imposed a defect, with the probability determined by C₀, the defect concentration of the lattice. A random number R₁ is generated and this site is attached with a defect (c⫽1) if R₁⬍C₀, and not (c⫽0) otherwise. The defect type to this site is measured by a random number R₂ within共⫺0.5, 0.5兲 associated with a choice of b_mwithin (⫺b_M,b_M), where b_M is the maximal of b_m. The simulation begins at an extremely high temperature T⫽8.0 at which no freezing effect remains (T⁰⫽3.0). The simulation follows the standard kinetic ME- TROPOLIS algorithm for dipole flip among the all candidate states. For a site chosen randomly, all those energy terms defined in Eqs.共1兲–共4兲are calculated for the whole lattice to obtain H⫽H₀from Eq.共5兲and then this site is imposed with a new dipole chosen randomly to simulate the dipole flip.

Those energy terms are calculated again and we obtain H

⫽H₁. Consequently, a third random number R₃is generated to compare with probability p defined below:

p⫽exp关⫺共H₂⫺H₁兲/T兴 if H₂⬎H₁,

p⫽1 if H₂⭐H₁, 共7兲

where T is the temperature scaled by energy 共kT兲 and the Boltzmann constant k is omitted. If R₃⬍p, the dipole flip is approved and rejected otherwise. Then one cycle is com- pleted. This cycle is repeated until a given number of cycles is reached.

The time of simulation is scaled by the Monte Carlo step 共MCS兲and one MCS represents L⫻L chosen flip events. In our simulation, at each temperature, the initial 600 MCS runs are discarded and then the configuration averaging is per- formed over the subsequent 2500 MCS. Here it should be pointed out that such a short time for configuration averaging is not long enough for simulating many statistical mechani- cal phenomena. However, it was verified that for spin-glass- like systems such as relaxors studied here, the short-time Monte Carlo simulation gave excellent agreement with experiments.¹⁸In fact, we performed one averaging sampling over a time series as long as 40 000 MCS and no substantial difference of the data from our short-time data was found. In addition, the data presented below represent an averaging over four runs with different seeds for random number gen- erator of both the initial lattice and defect distribution. Under an ac-electric field of frequency␻, the lattice dielectric sus- ceptibility␹is written as¹⁹

␹

⬘

^⫽K_N

冓 ^兺

^{Ni 1⫹共1␻␶兲2}

冔

^,

␹

⬙

⫽K

N

冓 ^兺

^{Ni 1⫹␻␶共␻␶兲2}

冔

^{, 共8兲}

where具 典 represents the configuration averaging, ␹

⬘

^and␹

⬙

are the real and imaginary parts of ␹^, ␶ is the time for the dipole at site i flipping from one state to another, N⫽L² and K is a temperature-dependent constant. Because the lattice is inhomogeneous once C₀⬎0, time␶becomes site dependent and it can be expressed in the Arrhenius form¹⁹

␶⫽␶0exp共⌬H/T兲⫽␶00exp共⫺f₀/T兲exp共⌬H/T兲, 共9兲 where␶00is the preexponential factor scaling the character- istic time for lattice vibration,␶0is the typical flip time for a noninteracting system which should be dependent of the de- fect concentration C₀, f₀ is the Landau potential for a de- fective system of no dipole-interaction, and⌬H is the differ- ence in Hamiltonian after and before the dipole flip. Under the assumption of nondipole interaction, f₀ is actually the Landau potential f_L in the mean-field approximation. Ex- cluding the higher-order terms共fourth order and above兲, one has f₀⬃A₁P₀²⬃关␣^(T⫺T⁰)⫹b_mC₀兴P₀², where P₀ is the av- eraged dipole moment 共magnitude兲. Here, ␶00⫽1.0 and ␶

⫽␶0/n is obtained from the statistics of the MC sequences where n is the number of dipole flips at site i per MCS.

In the simulation, b_M is given and C₀ is treated as vari- ables to investigate. The other lattice parameters are chosen and the dimensionless normalization of them is done follow- ing the works by Hu et al. on the dynamics of domain switching in BaTiO₃.¹⁶These parameters are given in Table I. In addition, the external electric field takes the form E

⫽E₀⫹E_msin(2␲␻t), where E₀ is the dc bias, E_m is the ac- signal amplitude, and t is time.

Figure 1 presents the simulated dipole configurations at T/T⁰⫽0.3 for three lattices of various defect concentration C₀. As C₀⫽0.0 关Fig. 1共a兲兴, a normal ferroelectric configu- ration with multidomained structure is shown. The parallel dipole alignment within each domain and the 90° head-to-tail domain walls can be identified. All dipoles within the do- mains have similar moment while those at the walls are smaller in moment. As C₀⬎0 关C₀⫽0.4, Fig. 1共b兲兴, one sees clearly the lattice inhomogeneity and the moment of those dipoles near the walls begins to shrink, while it is interesting to note that the defects are randomly distributed in the lattice.

TABLE I. System parameters used in the simulation.

Parameter Value Parameter Value Parameter Value

T⁰ 3.0 ␣ 1.0 A11 ⫺0.5

A₁₂ 9.0 A₁₁₁ 0.8 G₁₁ 1.0

G₁₄ 0.2 G₄₄ 1.0 L 64

b_M 10 C₀ 0–1.0

FIG. 1. Snapshot dipole configuration of mode ferroelectric lat- tice with different defects concentration C₀ at T/T⁰⫽0.3. E₀

⫽0.0, Em⫽0.2, and␻⫽0.01.

As C₀⫽0.8关Fig. 1共c兲兴, the two-phase coexistence picture in the lattice becomes quite clear. The lattice consists of local ferroelectric regions embedded in the matrix of paraelectric phase, a typical pattern for relaxor ferroelectrics. We also observe the temperature dependence of this two-phase coex- isted structure, as shown in Fig. 2, where C₀⫽0.8. While the ferroelectric phase dominates over the lattice as T/T⁰⫽0.2, at T/T⁰⫽0.7 the ferroelectric phase becomes minor and its size is much smaller than that at T/T⁰⫽0.2. As T/T⁰⫽1.4, the ordered dipole clusters become too small to be easily identified. Anyhow, one is shown that some small-sized di- pole clusters remain stable at a temperature much higher than the stability point T⁰. Therefore, the present model produces a lattice configuration of dipoles consistent with our compre- hensive understanding of relaxors.

The evaluated dielectric susceptibility as a function of T for several lattices of different defect concentrations as indi- cated is presented in Fig. 3共a兲for real part and Fig. 3共b兲for imaginary part. These curves are shifted for a clear illustra- tion, and in fact the value of these curves at T⫽6.0 are very close to each other. The bigger C₀ is, the slightly higher the value is. As C₀⫽0, a typical ferroelectric phase transition

with sharp dielectric peak ␹m

⬘

^{at T}⬃T_m is observed. With increasing defect concentration, the temperature dependence of both␹

⬘

^and␹

⬙

becomes diffusive and the transition peak shifts upward and T_m is slightly up to a higher value, indi- cating the typical dielectric characteristics for relaxors. As C₀⬎0.6, the temperature range covered by the ferroelectric transition is already very broad. It is well known that for a relaxor the dielectric susceptibility above T_m can be de- scribed by the following equation:^20,21

1

␹

⬘

^⫺ 1

␹m

⬘

^{⫽C⫺1共T⫺Tm兲␥, 共10兲} where C is a material constant similar to the Curie-Weiss constant and ␥ is the transition exponent characterizing the diffusivity of the phase transition. The bigger␥is, the more diffusive the transition is. A best fitting of the data shown in Fig. 3共a兲by Eq.共10兲produces the fitted parameters␹m

⬘

^{, C,} T_m, and ␥ as a function of defect concentration C₀, as shown in Figs. 4共a兲and 4共b兲, respectively. While␹m

⬘

^remains roughly unchanged, constant C drops down slowly and␥^and T_m increases with C₀, until ␥⬃1.8 at C₀⫽0.6 and above.

This exponent is quite close to the experimentally evaluated values for several relaxors, such as ␥⫽1.64 for PMN and 1.76 for PZN 关Pb(Zr_1/3Nb_2/3)O₃兴.²⁰ This is a positive evi- dence to support the present GL model.

Our simulation reveals a stronger frequency dispersion of the dielectric susceptibility for a lattice of higher C₀. A lower ␹m

⬘

value and a higher T_m value are observed when frequency␻increases. We also simulate the effect of ac-field FIG. 2. Snapshot dipole configuration of mode ferroelectric lat-

tice at different normalized temperatures T/T⁰ as C0⫽0.8. E0

⫽0.0, E_m⫽0.2, and␻⫽0.01. The cycles label the dipole clusters.

FIG. 3. Simulated dielectric susceptibilities␹⬘共a兲and␹⬙共b兲as a function of temperature kT for different defect concentration C₀ 共from bottom to top, C0⫽0.0, 0.1, 0.2, 0.4, 0.6, 0.8, and 1.0兲. E0

⫽0.0, E_m⫽0.2, and␻⫽0.01. The curves are shifted for eye guid- ance共see text兲.

amplitude and dc-electric bias on the dielectric susceptibility.

It is revealed that the value of␹m

⬘

drops down and T_mshifts toward the higher value with increasing dc-electric bias. The effect of the ac-field amplitude is opposite to that of the dc-electric bias. With enhancing ac-field amplitude, the tran- sition peak shifts upward and T_mapproaches a lower value in the same time. All of these features are reflected in typical relaxors.

To understand this broadening behavior of the phase tran- sition as reflected by the dielectric susceptibility as a func- tion of T in a defective lattice, one may look at the evolution of the Landau potential and interaction terms f_L, f_dip, and f_Gwith defect concentration. In Fig. 5 are shown these terms as a function of T, respectively. As C₀⫽0, each of these terms as a function of T can be divided into two temperature regions: the paraelectric region and ferroelectric region, separated by a clean boundary around which the phase tran- sition occurs. A prominent feature as reflected upon increas- ing C₀ is the obvious smearing of the boundary region. For f_L, the slope jump of the linear relation is disappearing with increasing C₀. The rapid change of both f_dipand f_Gover the low-T range is weakened as the lattice contains more defects.

It is noted that in the low-T range, the dipole-dipole interac- tion and the gradient energy are lifted, and as a compensation a big drop of the Landau free energy is observed.

For the perfect lattice of no any defect, the ferroelectric phase transitions are mainly determined by coefficient A₁ in Eq. 共1兲. As T⬍T⁰, the paraelectric phase loses its stability.

The introduction of defects into the lattice generates an in- homogeneity over the lattice where the stability for paraelec- tric phase varies from site to site. Some sites favor ferroelec- tric phase as TⰇT⁰, while some others favor paraelectric phase as TⰆT⁰. Therefore, the essence of the GL model for relaxors is to broaden the temperature range at which the

ferroelectric phase transitions proceed and the recorded tran- sition region extends towards both the low-T range and high-T range in the same time.

Unfortunately, to the best knowledge of the authors, there has never been an experimental system reported in which the defect concentration can be modulated to cover the whole composition range, so that a direct checking of the present model becomes possible. It has been reported recently that an irradiation of some ferroelectric copolymers by electrons, ions, or protons introduces defects into the systems and re- sults in an evolution of the dielectric behaviors from a nor- mal ferroelectric state to a relaxorlike state.²²However, such an irradiation basically suppresses the ferroelectric phase transitions, probably by amorphorization, while no dipole cluster can be stably retained at a temperature above the Curie point for the nonirradiated sample. Thus, such a defec- tive system seems not compatible with the present model.

In conclusion, we have performed a Monte Carlo simula- tion on the electric-dipole configuration and dielectric behav- FIG. 4. 共a兲Evaluated maximal value of␹⬘^,␹mand constant C,

and共b兲evaluated temperature T_mfor␹mand transition exponent␥ as a function of defect concentration C0. E0⫽0.0, Em⫽0.2, and

␻⫽0.01.

FIG. 5. Simulated Landau potential f_L, dipole-dipole interac- tion f_dip, and gradient energy f_G 共per site兲 as a function of tem- perature kT, respectively, for different defect concentration C₀ 共from bottom to top, C0⫽0.0, 0.1, 0.2, 0.4, 0.6, 0.8, and 1.0兲. E0

⫽0.0, E_m⫽0.2, and␻⫽0.01.

ior of a Ginzburg-Landau ferroelectric lattice with randomly distributed defects. It has been revealed that introduction of the defects results in a gradual evolution of the system from a normal ferroelectric state to a typical relaxor state, charac- terized by the diffusive phase transitions, strong frequency dispersion, and enhancement of the dielectric susceptibility.

A smearing effect of the Landau potential, the dipole-dipole interaction, and the gradient energy over the phase transition region has been observed. It is suggested that the present

Ginzburg-Landau model represents a realistic approach to the phase transition and dielectric property of relaxor ferro- electrics.

The authors would like to acknowledge the financial sup- port from the Natural Science Foundation of China through the innovative group project and Project No. 50332020, the National Key Project for Basic Research of China共Grant No.

2002CB613303兲, and LSSMS of Nanjing University.

*Author to whom correspondence should be addressed. Electronic address: [email protected]

1G. Smolenski and A. Agranovska, Sov. Phys. Solid State 1, 1429 共1960兲.

2L. E. Cross, Ferroelectrics 76, 241共1987兲.

3D. Viehland, S. J. Jang, L. E. Cross, and M. Wutting, J. Appl.

Phys. 68, 2916共1990兲; Phys. Rev. B 43, 8316共1991兲.

4A. K. Tagantsev, Phys. Rev. Lett. 72, 1100共1994兲.

5B. E. Vugmeister and M. D. Glinchuk, Rev. Mod. Phys. 62, 993 共1990兲.

6A. E. Glazounov, A. K. Tagantsev, and A. J. Bell, Phys. Rev. B 53, 11 281共1996兲.

7Y. Imry and S. K. Ma, Phys. Rev. Lett. 35, 1399共1975兲.

8R. Fisch, Phys. Rev. B 67, 094110共2003兲.

9K. Eichhorn and K. Binder, J. Phys.: Condens. Matter 8, 5209 共1996兲.

10H. Qian and L. A. Bursill, Int. J. Mod. Phys. B 10, 2027共1996兲.

11Z. Y. Cheng, R. S. Katiyar, X. Yao, and A. S. Bhalla, Phys. Rev.

B 57, 8166共1998兲.

12Q. Tan and D. Viehland, Phys. Rev. B 53, 14 103共1996兲.

13Z. R. Liu, B. L. Gu, and X. W. Zhang, Phys. Rev. B 62, 1共2000兲.

14W. Cao and L. E. Cross, Phys. Rev. B 44, 5共1995兲.

15S. Numbu and D. A. Sagala, Phys. Rev. B 50, 5838共1994兲.

16H. L. Hu and L. Q. Chen, Mater. Sci. Eng., A 238, 182共1997兲; J.

Am. Ceram. Soc. 81, 492共1998兲.

17S. Semenovskaya and A. G. Khachaturyan, J. Appl. Phys. 83, 5125共1998兲.

18K. Binder and A. P. Young, Rev. Mod. Phys. 58, 801共1986兲.

19C. C. Su, B. Vugmeister, and A. G. Khachaturyan, J. Appl. Phys.

90, 6345共2001兲.

20G. A. Smolensky, J. Phys. Soc. Jpn. 28, 26共1997兲.

21K. Uchino and S. Nomura, Ferroelectr., Lett. Sect. 44, 55共1982兲.

22Q. M. Zhang, V. Bharti, and X. Zhao, Science 280, 2101共1998兲.