- Email: [email protected]
Relaxation dynamics of metastable systems: application to polar medium
Physica A 340 (2004) 196 – 200
www.elsevier.com/locate/physa
Relaxation dynamics of metastable systems: application to polar medium E. Klotins∗ Institute of Solid State Physics, 8 Kengaraga Str., Riga LV 1063, Latvia
Abstract Ginzburg–Landau theory for ferroelectric phase instability is connected with Fokker–Planck equation techniques to model impact of thermal noise in the kinetics of ferroelectric polarization. Distinguishing features of this approach accounting for spatially inhomogeneous external 0eld, 0nite size e2ects and boundary conditions are based on an auxiliary relationship addressing the problem to imaginary time Schr4odinger equation with well-developed numerical techniques for its solution. c 2004 Elsevier B.V. All rights reserved. PACS: 05:45: − a; 05:70: − a; 05.70.L; 77.80.D; 77.22.G Keywords: Nonequilibrium thermodynamics; Ferroelectrics; Polarization reversal
The e2ect of thermal noise on a variety of nonequilibrium phenomena has attracted much attention in a number of areas in physical science raising, in particular, conceptual and technical questions regarding the application of the thermodynamic approach to metastability of small systems. Time-dependent Ginzburg–Landau approach [1] and generalizations of Boltzmann–Gibbs statistical mechanics [2] are the concerns in this context. Crucial questions related to this type of analysis address the rate of polarization switching, dynamics of domain walls, the role of 0nite size, and dimensionality. The conventional mathematical technique is based on a variation equation for the Ginzburg–Landau energy functional and yields spatio-temporal relaxation of polarization in zero temperature limit whereas a more advanced approach accounting for 0nite temperatures is based on transformation of the problem to the Fokker–Planck equation for density distribution of polarization. Although the Fokker–Planck equation ∗
Tel./fax: +371-7112583. E-mail address: [email protected] (E. Klotins).
c 2004 Elsevier B.V. All rights reserved. 0378-4371/$ - see front matter doi:10.1016/j.physa.2004.04.007
E. Klotins / Physica A 340 (2004) 196 – 200
197
may be derived for the arbitrary energy functional, the methods for its integration are problem-speci0c with transformation to the imaginary time Schr4odinger equation as the most suEcient technique for nonconservative energy functional and 0nite size systems. For quartic energy functional and in case of spatial homogeneity this technique is well-developed [3–5] and yields thermal noise controlled dielectric relaxation, dielectric response on time-periodic driving voltage, and nonlinear susceptibility [6]. Extension of this imaginary time Schr4odinger equation approach to more realistic systems, for example ferroelectric materials, comprising electric domains and surfaces is a challenge and the key problem faced in this work. Formally, the spatial extension is introduced by a gradient term in the energy functional contributing in an auxiliary function (ansatz) coupling the multivariate density distribution of polarization with an analogue of wave function found by symplectic integration of the imaginary time Schr4odinger equation. The mathematical tools are real-space/real-time matrix recurrence relations well-developed for simulating quantum electron dynamics [7] by large-scale numerical integration. This work is restricted to relaxation rate of polarization found as gradient-speci0c and, consequently, giving a possible answer to the question of polarization switching that occurs at a much lower 0eld than the conventional approach [1] suggests. The approach starts with the probability density of electric polarization (P; x) derived from Ginzburg–Landau energy functional for uniaxial ferroelectrics (in dimensionless units) P(x)2 P(x)4 P (x)2 dx (1) H= −P(x) − + + 2 4 2 resulting in variational derivative of the energy functional Eq. (1) as H = − − P(x) + P(x)3 − P (x) : P(x)
(2)
In the spatially homogeneous case illustrated at the beginning, it is convenient to introduce volume density of energy U = −P −
P2 P4 + 2 4
(3)
and write the Fokker–Planck equation for density distribution of the polarization (P; t) as (P; ˙ t) =
#2 # (P; t) : [(P; t) U (P)] + #P #P 2
(4)
Here the di2usion coeEcient condenses most features of the system, namely, the coeEcients in Eq. (1), averaged time of barrier crossing, and the temperature. The ansatz for spatially homogeneous energy functional Eq. (1) is de0ned as (P; t) = eF(P) G(P; t) :
(5)
Substitution of the ansatz Eq. (5) in Eq. (4) constitutes di2erential equation for the F-function as 2 F (P) = U (P) and its solution is given by F(P) = U (P)=2 . As a
198
E. Klotins / Physica A 340 (2004) 196 – 200
result the ansatz eliminates the 0rst derivative of G(P; t) and generates the imaginary time Schr4odinger equation for the G(P; t)-function as 2 #2 U U ˙ G(P; t) = 2 − + G(P; t) : (6) #P 4 2 In case of nonconservative energy functional Eq. (1) given by = (t) an extra U˙ (P; t) term appears in Eq. (6) [8], estimated as negligible in semiadiabatic limit and omitted 2 here for brevity. Considering V (P) = −U =4 + U =2 as a local potential [3,4] the problem, Eq. (6), becomes similar to its quantum mechanical analogue for the wave function of electrons with well-developed symplectic integration methods for its solution [7]. The relevant recurrence relations for G(P; t)-function starts with the initial condition found by considering thermal equilibrium under driving voltage 1 as the initial state. Then, after alteration the driving voltage from 1 to 2 at the time t = 0, the initial condition for subsequent relaxation yields [9] −U (P; 1 ) U (P; 2 ) −U (P; 1 ) exp + dP = G(P; 0) : (7) 2 This mapping between the Fokker–Planck and the imaginary time Schr4odinger equations [3,4,10] is presented here as a background for extension of the problem toward multivariate probability density of polarization = (P1 ; : : : ; PM ; t) in spatial mesh x = xmin + Jx(m − 1) = xm , m = 1; 2; : : : ; M , Jx = (xmax − xmin )=M . In these terms the ansatz Eq. (5) transforms in M −U (Pm ; 2 ) + Pm P m (P1 ; P2 ; : : : ; PM ; t) = exp G(P1 ; : : : ; PM ; t) (8) 2 m=1
and Eq. (6) is replaced by the spatially extended form M M #2 ˙ 1 ; P2 ; : : : ; PM ; t) = G(P + Vm G(P1 ; P2 ; : : : ; PM ; t) ; #Pm2 m=1
(9)
m=1
here the V -potential is nonlocal and reads as Vm = −
[U (Pm )]2 U (Pm ) [Pm ]2 U (Pm )Pm + − + : 4 2 4 2
(10)
Going back to the initial condition, Eq. (7), one can 0nd for spatial extension M −2U (Pm ; 1 ) + U (Pm ; 2 ) + Pm Pm G(P1 ; P2 ; : : : ; PM ; 0) = C exp : 2
(11)
m=1
Here the normalization constant reads as M −U (Pm ; 1 ) + Pm P m dP1 ; dP2 ; : : : ; dPM : C(1 ) = exp m=1
(12)
a.u.
E. Klotins / Physica A 340 (2004) 196 – 200
199
1
V POTENTIAL
0
-1
-2
-3 -1
-0.5
0
POLARIZATION
0.5
1
a.u.
Fig. 1. V -potential landscape for variable driving voltage speci0ed by an external 0eld term in Eq. (1). The thin, medium and bold curves match = 0, = −1=20 and = −0 , correspondingly (here 0 is the static coercive 0eld). The e2ect of polarization gradient is proved similarly to enhancement of the applied 0eld.
Accounting for thermal noise in Eqs. (9)–(12) is a principal departure from the conventional Landau–Khalatnikov approach [1]. While in Ref. [1] the system obeys a macroscopic damped equation exhibiting cooperative phenomena, the Fokker–Planck equation, Eq. (4), for energy density, Eq. (3), represent the dynamics toward uniquely stationary state limit and model both the relaxation and the dynamic hysteresis [6,9]. Also the puzzling feature of 0eld-induced polarization switching that occurs at much lower values than the conventional approach may be interpreted by accounting for the thermal noise term. A technical development is the transformation of Eqs. (9)–(12) in recurrence matrix representation. In spatially homogeneous case, Eq. (3), N grid points are introduced in the P space with increment JP=(Pmax −Pmin )=N [9]. In these terms Pn =Pmin +JP(n− 1), n = 1; 2; : : : ; N + 1, the G-function transforms in a N vector, and the operator part in Eq. (9) in N × N matrices as the key entities for large-scale modelling of relaxation dynamics either in spatially homogeneous case, Eq. (6), or in case of weak nonlocality, Eqs. (9) and (10), similarly to this for symplectic integration of quantum Schr4odinger equations [7]. It is now possible to estimate the e2ect of gradient term on the relaxation rate of polarization as the second result of this work. Considering a model initial state speci0ed by zero driving 0eld and zero polarization one can 0nd the V -potential symmetric as shown by the thin plot in Fig. 1. More realistic is an inhomogeneous polarization 0eld for which both polarization P and spatial derivative P are zero at the boundaries between regions with opposite polarization. In vicinity of the boundaries the product PP may be estimated as a small quantity in Eq. (8) and the major contribution in probability density are given by gradient terms in the potential Eq. (10), namely, by the last U (Pm )Pm =2 term contributing in V -potentional symmetry violation. Without lost of generality, we identify a boundary at x = xm with Pm−1 ¡ 0, Pm+1 = −Pm−1 (1) and, obviously, Pm+1 ¡ 0. Then a negative driving 0eld m+1 ¡ 0 distorts the symmetry of V -potential as shown by the medium plot in Fig. 1 and the polarization Pm+1 is (1) (2) | ¡ |m+1 | the relaxed toward a negative value. By enlarging the driving 0eld |m+1
200
E. Klotins / Physica A 340 (2004) 196 – 200
V -potential transforms into a strongly asymmetric one as shown by the bold plot in Fig. 1, and the relaxation rate of Pm+1 toward a negative value grows if compared with (1) this for m+1 driving 0eld. Accounting for U (Pm )Pm =2 term in Eq. (10) one can
2
(Pm )] 0nd the impact of the [U 4 term enhanced in the vicinity of boundaries at which P exhibited a maximum. Consequently, the relaxation of polarization is initialized in the vicinity of domain walls that agree well with nucleation and domain wall motion mechanism in ferroelectrics studied in particular cases by other techniques [11,12] and give clues to the problem of polarization switching that occurs at much lower 0eld than the conventional approach suggests [1]. To summarize, in this paper we have considered the application of Fokker–Planck imaginary time Schr4odinger equation approach to the relaxation dynamics of spatially extended systems exampli0ed by Ginzburg–Landau-type energy functional in the limit of weak nonlocality. Key constituents of this approach includes probability density of electric polarization and an auxiliary relationship (ansatz) accounting, primarily following our knowledge, for spatial extension. The relaxation rate of polarization is found to be gradient speci0c, namely, the polarization switching to the opposite direction is favored in the vicinity of boundaries at which the variational derivative of polarization reaches its maximal value. It gives a possible answer to the question of polarization switching occurring at much lower 0eld than the conventional approach suggests. Another question concerned is the linkage of noise and cooperative phenomena formally given by Eq. (4)-type relation with a feedback term [13] its integration at arbitrary driving captured by routine symplectic integration. We therefore look forward to future applications of this approach for a variety of problems associated with thermally activated and spatially extended structural instability.
This work has been supported by the Contract No. ICA1-CT-2000-70007 of European Excellence Center of Advanced Material Research and Technology (Riga), J. Hlinka and the Institute of Physics ASCR (Praha), and A. Kholkin and University of Aveiro (Portugal) are greatly appreciated for hospitality, discussions and comments. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
R. Ahluwalia, T. Lookman, A. Saxena, W. Cao, cond-mat/0308232 v1, 2003. F.D. Nobre, C. Tsallis, Physica A 213 (1995) 337. B. Caroli, C. Caroli, B. Roulet, J. Statist. Phys. 21 (4) (1979) 415–437. A. Bunde, J.-F. Gouyet, Z. Phys. B: Conden. Matter 42 (1981) 169–181. P. Talkner, New J. Phys. 1 (1999) 4.1–4.25. E. Klotins, Ferroelectrics 295 (2003) 577–588. N. Watanabe, M. Tsukada, Phys. Rev. E 62 (2000) 2914. J. Hlinka, private communication. E. Klotins, J. Hlinka, J. Kaupuzs, Ferroelectrics, in press. G. Parisi, Statistical Field Theory, Frontiers in Physics, Lecture Notes Series, Addison-Wesley, Reading, MA, 1988. [11] M. Iwata, H. Orihara, Y. Ishibashi, Jpn. J. Appl. Phys. 40 (2001) 708–712. [12] A.M. Bratkovsky, A.P. Levanyuk, Phys. Rev. B 66 (2002) 184109. [13] M. Shiino, J. Math. Phys. 43 (2002) 2654–2669.
www.elsevier.com/locate/physa
Relaxation dynamics of metastable systems: application to polar medium E. Klotins∗ Institute of Solid State Physics, 8 Kengaraga Str., Riga LV 1063, Latvia
Abstract Ginzburg–Landau theory for ferroelectric phase instability is connected with Fokker–Planck equation techniques to model impact of thermal noise in the kinetics of ferroelectric polarization. Distinguishing features of this approach accounting for spatially inhomogeneous external 0eld, 0nite size e2ects and boundary conditions are based on an auxiliary relationship addressing the problem to imaginary time Schr4odinger equation with well-developed numerical techniques for its solution. c 2004 Elsevier B.V. All rights reserved. PACS: 05:45: − a; 05:70: − a; 05.70.L; 77.80.D; 77.22.G Keywords: Nonequilibrium thermodynamics; Ferroelectrics; Polarization reversal
The e2ect of thermal noise on a variety of nonequilibrium phenomena has attracted much attention in a number of areas in physical science raising, in particular, conceptual and technical questions regarding the application of the thermodynamic approach to metastability of small systems. Time-dependent Ginzburg–Landau approach [1] and generalizations of Boltzmann–Gibbs statistical mechanics [2] are the concerns in this context. Crucial questions related to this type of analysis address the rate of polarization switching, dynamics of domain walls, the role of 0nite size, and dimensionality. The conventional mathematical technique is based on a variation equation for the Ginzburg–Landau energy functional and yields spatio-temporal relaxation of polarization in zero temperature limit whereas a more advanced approach accounting for 0nite temperatures is based on transformation of the problem to the Fokker–Planck equation for density distribution of polarization. Although the Fokker–Planck equation ∗
Tel./fax: +371-7112583. E-mail address: [email protected] (E. Klotins).
c 2004 Elsevier B.V. All rights reserved. 0378-4371/$ - see front matter doi:10.1016/j.physa.2004.04.007
E. Klotins / Physica A 340 (2004) 196 – 200
197
may be derived for the arbitrary energy functional, the methods for its integration are problem-speci0c with transformation to the imaginary time Schr4odinger equation as the most suEcient technique for nonconservative energy functional and 0nite size systems. For quartic energy functional and in case of spatial homogeneity this technique is well-developed [3–5] and yields thermal noise controlled dielectric relaxation, dielectric response on time-periodic driving voltage, and nonlinear susceptibility [6]. Extension of this imaginary time Schr4odinger equation approach to more realistic systems, for example ferroelectric materials, comprising electric domains and surfaces is a challenge and the key problem faced in this work. Formally, the spatial extension is introduced by a gradient term in the energy functional contributing in an auxiliary function (ansatz) coupling the multivariate density distribution of polarization with an analogue of wave function found by symplectic integration of the imaginary time Schr4odinger equation. The mathematical tools are real-space/real-time matrix recurrence relations well-developed for simulating quantum electron dynamics [7] by large-scale numerical integration. This work is restricted to relaxation rate of polarization found as gradient-speci0c and, consequently, giving a possible answer to the question of polarization switching that occurs at a much lower 0eld than the conventional approach [1] suggests. The approach starts with the probability density of electric polarization (P; x) derived from Ginzburg–Landau energy functional for uniaxial ferroelectrics (in dimensionless units) P(x)2 P(x)4 P (x)2 dx (1) H= −P(x) − + + 2 4 2 resulting in variational derivative of the energy functional Eq. (1) as H = − − P(x) + P(x)3 − P (x) : P(x)
(2)
In the spatially homogeneous case illustrated at the beginning, it is convenient to introduce volume density of energy U = −P −
P2 P4 + 2 4
(3)
and write the Fokker–Planck equation for density distribution of the polarization (P; t) as (P; ˙ t) =
#2 # (P; t) : [(P; t) U (P)] + #P #P 2
(4)
Here the di2usion coeEcient condenses most features of the system, namely, the coeEcients in Eq. (1), averaged time of barrier crossing, and the temperature. The ansatz for spatially homogeneous energy functional Eq. (1) is de0ned as (P; t) = eF(P) G(P; t) :
(5)
Substitution of the ansatz Eq. (5) in Eq. (4) constitutes di2erential equation for the F-function as 2 F (P) = U (P) and its solution is given by F(P) = U (P)=2 . As a
198
E. Klotins / Physica A 340 (2004) 196 – 200
result the ansatz eliminates the 0rst derivative of G(P; t) and generates the imaginary time Schr4odinger equation for the G(P; t)-function as 2 #2 U U ˙ G(P; t) = 2 − + G(P; t) : (6) #P 4 2 In case of nonconservative energy functional Eq. (1) given by = (t) an extra U˙ (P; t) term appears in Eq. (6) [8], estimated as negligible in semiadiabatic limit and omitted 2 here for brevity. Considering V (P) = −U =4 + U =2 as a local potential [3,4] the problem, Eq. (6), becomes similar to its quantum mechanical analogue for the wave function of electrons with well-developed symplectic integration methods for its solution [7]. The relevant recurrence relations for G(P; t)-function starts with the initial condition found by considering thermal equilibrium under driving voltage 1 as the initial state. Then, after alteration the driving voltage from 1 to 2 at the time t = 0, the initial condition for subsequent relaxation yields [9] −U (P; 1 ) U (P; 2 ) −U (P; 1 ) exp + dP = G(P; 0) : (7) 2 This mapping between the Fokker–Planck and the imaginary time Schr4odinger equations [3,4,10] is presented here as a background for extension of the problem toward multivariate probability density of polarization = (P1 ; : : : ; PM ; t) in spatial mesh x = xmin + Jx(m − 1) = xm , m = 1; 2; : : : ; M , Jx = (xmax − xmin )=M . In these terms the ansatz Eq. (5) transforms in M −U (Pm ; 2 ) + Pm P m (P1 ; P2 ; : : : ; PM ; t) = exp G(P1 ; : : : ; PM ; t) (8) 2 m=1
and Eq. (6) is replaced by the spatially extended form M M #2 ˙ 1 ; P2 ; : : : ; PM ; t) = G(P + Vm G(P1 ; P2 ; : : : ; PM ; t) ; #Pm2 m=1
(9)
m=1
here the V -potential is nonlocal and reads as Vm = −
[U (Pm )]2 U (Pm ) [Pm ]2 U (Pm )Pm + − + : 4 2 4 2
(10)
Going back to the initial condition, Eq. (7), one can 0nd for spatial extension M −2U (Pm ; 1 ) + U (Pm ; 2 ) + Pm Pm G(P1 ; P2 ; : : : ; PM ; 0) = C exp : 2
(11)
m=1
Here the normalization constant reads as M −U (Pm ; 1 ) + Pm P m dP1 ; dP2 ; : : : ; dPM : C(1 ) = exp m=1
(12)
a.u.
E. Klotins / Physica A 340 (2004) 196 – 200
199
1
V POTENTIAL
0
-1
-2
-3 -1
-0.5
0
POLARIZATION
0.5
1
a.u.
Fig. 1. V -potential landscape for variable driving voltage speci0ed by an external 0eld term in Eq. (1). The thin, medium and bold curves match = 0, = −1=20 and = −0 , correspondingly (here 0 is the static coercive 0eld). The e2ect of polarization gradient is proved similarly to enhancement of the applied 0eld.
Accounting for thermal noise in Eqs. (9)–(12) is a principal departure from the conventional Landau–Khalatnikov approach [1]. While in Ref. [1] the system obeys a macroscopic damped equation exhibiting cooperative phenomena, the Fokker–Planck equation, Eq. (4), for energy density, Eq. (3), represent the dynamics toward uniquely stationary state limit and model both the relaxation and the dynamic hysteresis [6,9]. Also the puzzling feature of 0eld-induced polarization switching that occurs at much lower values than the conventional approach may be interpreted by accounting for the thermal noise term. A technical development is the transformation of Eqs. (9)–(12) in recurrence matrix representation. In spatially homogeneous case, Eq. (3), N grid points are introduced in the P space with increment JP=(Pmax −Pmin )=N [9]. In these terms Pn =Pmin +JP(n− 1), n = 1; 2; : : : ; N + 1, the G-function transforms in a N vector, and the operator part in Eq. (9) in N × N matrices as the key entities for large-scale modelling of relaxation dynamics either in spatially homogeneous case, Eq. (6), or in case of weak nonlocality, Eqs. (9) and (10), similarly to this for symplectic integration of quantum Schr4odinger equations [7]. It is now possible to estimate the e2ect of gradient term on the relaxation rate of polarization as the second result of this work. Considering a model initial state speci0ed by zero driving 0eld and zero polarization one can 0nd the V -potential symmetric as shown by the thin plot in Fig. 1. More realistic is an inhomogeneous polarization 0eld for which both polarization P and spatial derivative P are zero at the boundaries between regions with opposite polarization. In vicinity of the boundaries the product PP may be estimated as a small quantity in Eq. (8) and the major contribution in probability density are given by gradient terms in the potential Eq. (10), namely, by the last U (Pm )Pm =2 term contributing in V -potentional symmetry violation. Without lost of generality, we identify a boundary at x = xm with Pm−1 ¡ 0, Pm+1 = −Pm−1 (1) and, obviously, Pm+1 ¡ 0. Then a negative driving 0eld m+1 ¡ 0 distorts the symmetry of V -potential as shown by the medium plot in Fig. 1 and the polarization Pm+1 is (1) (2) | ¡ |m+1 | the relaxed toward a negative value. By enlarging the driving 0eld |m+1
200
E. Klotins / Physica A 340 (2004) 196 – 200
V -potential transforms into a strongly asymmetric one as shown by the bold plot in Fig. 1, and the relaxation rate of Pm+1 toward a negative value grows if compared with (1) this for m+1 driving 0eld. Accounting for U (Pm )Pm =2 term in Eq. (10) one can
2
(Pm )] 0nd the impact of the [U 4 term enhanced in the vicinity of boundaries at which P exhibited a maximum. Consequently, the relaxation of polarization is initialized in the vicinity of domain walls that agree well with nucleation and domain wall motion mechanism in ferroelectrics studied in particular cases by other techniques [11,12] and give clues to the problem of polarization switching that occurs at much lower 0eld than the conventional approach suggests [1]. To summarize, in this paper we have considered the application of Fokker–Planck imaginary time Schr4odinger equation approach to the relaxation dynamics of spatially extended systems exampli0ed by Ginzburg–Landau-type energy functional in the limit of weak nonlocality. Key constituents of this approach includes probability density of electric polarization and an auxiliary relationship (ansatz) accounting, primarily following our knowledge, for spatial extension. The relaxation rate of polarization is found to be gradient speci0c, namely, the polarization switching to the opposite direction is favored in the vicinity of boundaries at which the variational derivative of polarization reaches its maximal value. It gives a possible answer to the question of polarization switching occurring at much lower 0eld than the conventional approach suggests. Another question concerned is the linkage of noise and cooperative phenomena formally given by Eq. (4)-type relation with a feedback term [13] its integration at arbitrary driving captured by routine symplectic integration. We therefore look forward to future applications of this approach for a variety of problems associated with thermally activated and spatially extended structural instability.
This work has been supported by the Contract No. ICA1-CT-2000-70007 of European Excellence Center of Advanced Material Research and Technology (Riga), J. Hlinka and the Institute of Physics ASCR (Praha), and A. Kholkin and University of Aveiro (Portugal) are greatly appreciated for hospitality, discussions and comments. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
R. Ahluwalia, T. Lookman, A. Saxena, W. Cao, cond-mat/0308232 v1, 2003. F.D. Nobre, C. Tsallis, Physica A 213 (1995) 337. B. Caroli, C. Caroli, B. Roulet, J. Statist. Phys. 21 (4) (1979) 415–437. A. Bunde, J.-F. Gouyet, Z. Phys. B: Conden. Matter 42 (1981) 169–181. P. Talkner, New J. Phys. 1 (1999) 4.1–4.25. E. Klotins, Ferroelectrics 295 (2003) 577–588. N. Watanabe, M. Tsukada, Phys. Rev. E 62 (2000) 2914. J. Hlinka, private communication. E. Klotins, J. Hlinka, J. Kaupuzs, Ferroelectrics, in press. G. Parisi, Statistical Field Theory, Frontiers in Physics, Lecture Notes Series, Addison-Wesley, Reading, MA, 1988. [11] M. Iwata, H. Orihara, Y. Ishibashi, Jpn. J. Appl. Phys. 40 (2001) 708–712. [12] A.M. Bratkovsky, A.P. Levanyuk, Phys. Rev. B 66 (2002) 184109. [13] M. Shiino, J. Math. Phys. 43 (2002) 2654–2669.