To the theory of relaxation phenomena in quasi-one-dimensional ferroelectrics with hydrogen bonds . Nonequilibrium statistical operator approach

Kinetic equations for the dielectric properties of CsH2PO4-type ferroelectrics are derived within the nonequilibrium statistical operator method proposed by D.N.Zubarev. These equations are solved with two different approaches. In the first approach the closure relation in the spirit of a two-site cluster approximation is used. In the second approach the exact relation between the static correlation functions is extrapolated to a dynamic case and used as a closure relation. The both approaches give similar results. The obtained dynamic characteristics are close to those of the Glauber model and are shown to describe the dynamic dielectric permittivity of CsH2PO4 and CsD2PO4 crystals.


Introduction
The Ising model (IM) appears to be very effective for the formulation and solution of problems in the theories for electrically or magnetically ordered compounds, binary alloys, etc. [1][2][3][4][5][6][7][8].In the case of low-dimensional structures, the model allows an accurate calculation of the thermodynamic quantities (order parameters, static susceptibility, specific heat) and the spin correlation functions [7][8][9].Accurate solutions of the two-dimensional IM [2,9] and the one-dimensional Kac model make a substantial contribution to the understanding of the nature of the second order phase transitions.
Great interest is given to the dynamics of the IM.The study of Ising systems' kinetics was initiated by Glauber's paper [10], where a one-dimensional (1D) model was treated.The papers [11][12][13] extend Glauber's approach to 2D and 3D Ising lattices.As in these cases the problem becomes rather intractable, the mean field approximation [11] and the cluster approximation [12,13] had to be used.Paper [13] differs from the preceding work [12] due to the account of a long-range interaction.Glauber's idea also inspired studies of the model with several modifications [14,15].Within Glauber's approach one can derive [32] a chain of coupled equations for time-dependent distribution functions (DF).Their solution can be found only for 1D model.For pseudo 1D model Zumer [20] suggested the simplest decoupling of the chain, which leads to some interpolative equation for one-site DF that yields the generally known exact result for static susceptibility.But Zinenko [30] discovered that approximations of papers [21,22] lead to incorrect results in the lowtemperature region.Paper [32] suggests a detailed analysis of the coupled equations for time-dependent DF obtained from the Glauber master equation.Several decoupling schemes are suggested.One of them yields thermodynamic and dynamic quantities of the model similar to those obtained within a cluster approach developed in papers [28,29].These theoretical results were extensively used for the discussion and analysis of experimental data on ferroelectric of CsH 2 PO 4 type [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32].
The crystal structure of CsH 2 PO 4 can be described as "PO 4 tetrahedrons connected by hydrogen bonds of two types".The shorter H-bonds connect PO 4 groups forming zig-zag chains along the b-axis.In the paraelectric phase the protons of these bonds are distributed symmetrically with respect to the centre of an H-bond, whereas in the ferroelectric phase there appears a spontaneous asymmetry.The protons of the longer H-bonds are ordered in both phases.Deuteration significantly changes the temperature of the phase transition, the dielectric and thermal properties of CsH 2 PO 4 crystal [16][17][18][19][20][21][22][23][24][25][26][27][28][29].It is the evidence of the important role of H-bonds in the ferroelectricity of these crystals.
Ferroelectric crystals PbHPO 4 contain only short H-bonds connecting PO 4 groups along the c-axis [17,18,19,33,34,35].The H-bonds in these crystals are almost parallel to the direction of the spontaneous polarization.The phase transition in these crystals is also closely connected with the proton ordering.However, the effects of deuteration in PbHPO 4 crystals are much less pronounced in comparison with CsH 2 PO 4 crystals.
It is nowadays commonly accepted that the ordering of protons (deuterons) in the ferroelectrics of CsH 2 PO 4 type is of a quasi-one-dimensional nature implying strong intrachain and weak interchain interactions.This results in some specific properties of such compounds which should be explained on a microscopic level.
Papers [16][17][18][19][20][21][22][23][24][25][26][27] are focused on the discussion of the selected experiment and the suggested model parameters which provide the best fit to those data.Such a one-sided approach sometimes leads to incorrect conclusions.Thus, the authors of papers [22][23][24][25], comparing predictions of the model with such parameters with other measured quantities, found a marked disagreement and concluded that the model cannot describe the compounds under investigation.However, papers [28,29,32] show that the theory [28,32] which allows for short-and long-range interactions, describes thermodynamics and dynamics of these compounds with a satisfactory accuracy.
Glauber's approach, however, has several deficiencies.This is, first of all, lack of substantiation for the master equation, which contains phenomenologically chosen spin-flip probabilities and unknown kinetic parameters that should be provided by a separate ab-initio theory.Such drawbacks are not intrinsic to the kinetic equations method [33][34][35][36][37][38][39][40][41][42][43].There exists a number of approaches in this field.In this paper we shall use a nonequilibrium statistical operator (NSO) method suggested by D.N.Zubarev.In the presence of a small parameter, the NSO method allows one to obtain kinetic equations for the model [41][42][43].We intend to use the method in order to describe the dynamics of quasi-1D ferroelectrics with hydrogen bonds.In section 2 a kinetic equation is derived.It is solved in section 3 within a two-site cluster approximation.In section 4 a chain of coupled equations for time-dependent DFs is derived and a decoupling scheme is proposed, which supports the results of the cluster approximation.If the kinetic parameters are independent of frequency, the NSO method yields the results that coincide with those [28,32] obtained within the Glauber master equation approach.

Nonequilibrium statistical operator. Kinetic equation
In this section we shall consider relaxation phenomena in systems with the Hamiltonian where Here Ĥ1 (s) denotes the Hamiltonian of the quasispin subsystem, Ĥ2 (t) is a Hamiltonian of the dissipative subsystem, H ′ (s, t) describes their interaction.Non-spin degrees of freedom constitute a heat bath, i.e. their heat capacity is so large that any perturbation caused by the quasispin subsystem cannot push the heat bath out of the equilibrium.The nonequilibrium statistical operator of the entire system satisfies the Liouville equation where iL is a Liouville operator defined by In order to calculate NSO for the system with a small parameter, the following form of (2.3) is convenient [43]: Solution of the Liouville equation (2.5) (or (2.3)) requires boundary conditions.Following the Zubarev NSO method [40] we shall look for such solutions of (2.5) which depend on time only via several observable variables.These variables of a reduced description set the nonequilibrium state of the system.The desired solutions can be selected, including an infinitely small source in the right-hand side of equation (2.5) [40]: Here ρq (t) is a quasiequilibrium statistical operator that has to minimize the entropy of the system at fixed variables of a reduced description and the normalization Spρ q (t) = 1.The source in the right-hand side of (2.6) violates the timereversal symmetry of the equation and selects retarded solutions corresponding to a reduced description.After the thermodynamic limit one has to set ε = 0 and to obtain the solution of the equation (2.5).
The described above approach is based on Bogolubov's idea that the nonequilibrium dynamics of the macroscopic physical system is characterized by a hierarchy of relaxation times [33]: parts of the system reach the equilibrium long before the entire system does.Therefore, these parts can be described by a small number of quasiintegrals of motion which slowly change under the action of weaker interactions that equilibrate all the subsystems.Thus, for times that are longer than some relaxation times, the nonequilibrium dynamics of the system can be described by a small number of parameters.In such a way a reduced description becomes sufficient [33].Let Pm be a set of reduced description variables.The quasiequilibrium statistical operator takes the form [40]: where ensures the normalization of ρq (t), F m (t) being the corresponding Lagrange parameters that have to be found from self-consistency conditions where These conditions lead to It means that F m (t) are conjugated to Pm .An expression for entropy clarifies the physical meaning of the parameters F m (t): and use the identity then integral equation (2.16) takes the form: This expression for NSO coincides with that obtained by Pokrovskii [41], if the following relation is taken into account: The evolution of the free system is determined by quantities a ml Transfer equations can be written using these parameters: If one calculates the last term in (2.23) using a linear in interaction V approximation for NSO (for the details of that calculation the reader is referred to [41]), the following kinetic equations for Pm are obtained (to the second order approximation with respect to V ): where L( 1) The kinetic equation (2.24) can be rewritten in the form: (2.28) Now we must be more explicit and write the quasispin-heat bath interaction term H ′ (s, t) of the Hamiltonian.Let us consider only one-particle interactions where U α j (f ) are operators related to the heat bath, and Let us introduce (the equilibrium) statistical operator of the heat bath and operator Ĥ′ (s, f ) defined by whose evolution is of the form: Using (2.29) one finds that With the help of the time Fourier transformation Hereafter we assume that fluctuations of heat bath variables at different sites do not correlate.Using also a translational invariance of the model we obtain: If there exists a diagonal basis for Ĥ2 (f ), where matrix elements of operator U α (f ) are real, the following relations hold [44]: (2.40) In this case and (2.37) can be written in the form where the following notations are used: Using a spectral representation for time-dependent correlation functions in (2.43) If correlations in the dissipative subsystem have a Markovian character and M α ′ α (t) takes the form [7] where τ f α are correlation times, then coefficients N and K can be related using (2.43), (2.46) and (2.47).We obtain: Finally, The terms with Ω α µ = Ω α ′ ν give contributions of order K α ′ α (Ω)/Ω to the solution of the kinetic equation.Such contributions are of no importance in many cases [7].Condition Ω α µ = Ω α ′ ν , in its turn, allows one to keep only the terms with α = −α ′ .The last condition can be proven in a very general form [7,42,43].The relaxation term (2.50) takes the simplest form if the set { Pm } contains only diagonal operators.In this case imaginary terms in (2.50) and R 2 ( Pm ) vanish, and the kinetic equation (2.27) takes the form [7,42,43]: Here the following notations are used: It is equation (2.51) that will be used in the following section in order to describe the relaxation dynamics of deuterated quasi-1D ferroelectrics with hydrogen bonds.

Dynamics and thermodynamics of deuterated quasi-1D ferroelectrics with hydrogen bonds in the two-site cluster approximation
We shall consider a system of protons located at O-H. ..O bonds constituting zig-zag chains.Elementary cells of such compounds contain two PO 4 tetrahedrons with two short H-bonds related to one of them (tetrahedron of type A).H-bonds related to the other (B-type) tetrahedron belong to two nearest structural elements surrounding it [28].The Hamiltonian of the proton subsystem of the considered ferroelectrics has the form [28]: Here the first term describes short-range configurational interactions of protons in chains near the A-type tetrahedron (the first Kronecker symbol) and the B-type one (the second Kronecker); the second term is an effective long-range interaction between the protons including interactions mediated by lattice vibrations, the last term includes an interaction of protons with a longitudinal external field, σ z R i f = ±1 is an operator of the z-component of the quasispin, which we juxtapose to the proton located in elementary cell R i at the f th H-bond (f = 1, 2); r 2 is a vector of a relative position of H-bonds in the cell; µ denotes an effective dipole moment of the elementary cell along the polar axis.We shall study thermodynamics and dynamics of the systems with the Hamiltonian (3.1) in the two-site cluster approximation [12,13,28].The approximation is based on the one-site and two-site Hamiltonians [28]: Ĥ(1) Here ∆ R i f denotes an effective field generated by those neighbouring H-bonds which are not included in the cluster, ν 1 2 η (1) is a molecular field generated by remote H-bonds via a long-range interaction.The fields ∆ R i f are variational parameters that have to be found via self-consistency equations.Following the Hamiltonian H (2) the evolution of quasispin operators is given by or in a frequency representation where the eigenfrequencies are and On the other hand, the Hamiltonian H (1) suggests ) )e −itαΩ (1) . (3.9) In TCA, operators Pm coming into kinetic equation (2.51) are Now one can calculate Q ± R j µα (P m ) and, taking into account the boundary conditions, obtain a system of equations for one-site and two-site DFs within the "H (2) approach": and an equation for one-site DF within the "H (1) approach" d dt η (1) = −2K (1) η (1) + 2K (1) tanh βΩ Here the following notations are used: 1 z (2) 1 z (2) (2) −1 ), where Thus, the TCA algorithm leads to a system of nonlinear differential equations (3.12), (3.13) for variables η (1) , η (2) , ∆.This system seems to be quite intractable.Therefore, we consider only the case of small deviations from equilibrium values η(1) , η(2) , ∆, Ẽ and neglect nonlinear in η (1) t , ∆ t , E t terms in the kinetic equations.This leads to separate equations for static quantities and dynamic deviations t + a 2 η (2) where Equations (3.20) coincide with the results of the conventional static theory.The general solution of (3.21) is of the form: with the relaxation times and (3.30) The quantity C(iω) determines the dynamic susceptibility of the model where Consequently, the dynamic permittivity is of the form: where ε ∞ is a high-frequency contribution to permittivity.Let us note that, if , the obtained in this section kinetic equations (3.12) coincide with that used in the stochastic Glauber model in [28].Therefore, Glauberian equations describe a physical system where Fourier-transforms of the heat bath correlation functions do not depend on frequency [7,42,43].In this case α is a constant of time dimensionality which characterizes an interaction of the spin subsystem with the heat bath.Now let us note that the cluster approach yields exact results (obtained in [5,[7][8][9][10]) for the static properties (free energy, specific heat, static permittivity ε(0, T ) (see (3.33)), order parameter (3.20)), equilibrium distribution functions [32,37] of quasi-one-dimensional systems.For the results of the cluster approach for the systems of higher dimensions the reader is referred to [13].In the case of a zero long-range interaction, formula (3.33) yields the result of [12].
To describe the property of some compound on the basis of the obtained above results one has to consider a mechanism of relaxation of that compound and calculate the kinetic parameters K (2)   µ and K (1) .We shall not discuss this problem here.Figures 1 and 2 show that dynamic permittivity of CsH 2 PO 4 and CsD 2 PO 4 crystals in the paraelectric phase can be described within the Glauberian approach under certain choice of model parameters (see table 1).To describe the dynamic permittivity in the ferroelectric phase one should suppose that bare relaxation time α depends on temperature.
In the theory of the CsH 2 PO 4 -type ferroelectrics it is important to describe dramatic isotope effects in these crystals.Partial substitution of hydrogens by deuterons can be described within the disordered Ising model.In such a way a problem of the calculation of dynamic properties of the disordered Ising model arises.Papers [49][50][51][52] treat the Glauberian dynamics of the model within the twosite cluster approximation.NSO dynamics of the disordered Ising model is an open problem.
In the next section we shall consider alternative methods for the solution of the kinetic equation in order to verify the cluster approach developed here.

Relaxation dynamics of the quasi-one-dimensional Ising model
For the sake of convenience we write the Hamiltonian of quasi-1D ferroelectrics of the order-disorder type in the following form [20,32] The first term describes the short-range intrachain interactions between the hydrogen bonds, the second one accounts for their effective long-range interactions, the last term is an interaction of hydrogen bonds with an external longitudinal electric field; S z ij is a quasispin operator corresponding to the hydrogen bond The spatial position of each quasispin is described by two indices (ij).The first one denotes the sequential number of the quasispin in a chain, the second is a chain number.
We shall treat the long-range interactions in the mean field approximation.In  1, the symbols present the experimental data [48,25,47].
this case the effective Hamiltonian takes the form: Here the notations are used.The time evolution of the quasispin operator is defined by (2.33): where  1), the symbols show experimental data [25].
The choice of operators P k , which come into the kinetic equation, is determined by the system's properties and the form of its Hamiltonian.In our case the following choice of P k is convenient: The kinetic equation (2.51) takes the form: where Now, on the basis of (4.9) one can obtain a system of equations for time-dependent spin correlation functions where the following notations are introduced: η m = σ z ij σ z i+m,j σ z i+(m+1),j = η m,m+1 , η 1,2 = η 3 , (4.14) η q = σ z ij q , η iq = ... q , i = 1, 2, 3.The system of equations (4.11) is not closed.Unfortunately, we do not know other relations between the nonequilibrium quasispin distribution functions of the model.The equilibrium spin distribution functions for the quasi-one-dimensional Ising model (4.1) were studied in papers [2,5,[7][8][9].The exact results for several correlations functions are known [7,8]: ) One can easily verify that all the correlation functions in (4.16) and (4.17) can be expressed in terms of the equilibrium averages η 0 and η 1 (see [7,8]): In view of relations (4.18) we can suggest closure relations for the kinetic equations (4.11): we derive a closed system of kinetic equations The system of equations (4.20) is nonlinear and quite complex.Therefore, we restrict further consideration to small deviations from the equilibrium and represent the distribution functions as sums of their equilibrium values and (small) timedependent deviations: In this case the following relations hold: Z where where Equations (4.25) yield the following expression for dynamic susceptibility of the system: In the paraelectric phase our results have a particularly simple form: .
Like in the preceding section, in the case K 0 = K 1 = K −1 , expressions (4.26) and (4.27) become simplified and we obtain the results of the Glauber approach derived in [32].

Conclusions
In this paper we study relaxational phenomena in quasi-one-dimensional ferroelectrics of the order-disorder type.We use Zubarev's NSO method in order to describe the kinetics of these compounds.Two closure relations are used to deal with the obtained kinetic equations.The first one uses the two-site cluster approximation in order to calculate the dynamic and thermodynamic characteristics of the system.The second one uses suggested in [32] decoupling for a chain of equations for quasispin distribution functions.The both approaches give numerically close results.For the description of experimental data within the NSO method one must calculate kinetic parameters K (2)  µ and K (1) .In the course of such a calculation one must consider mechanisms of relaxation in given crystals.Such problems will be considered in our future works.Here we only note that if the kinetic parameters are independent of frequency, our results coincide with those obtained in [28,32] within the Glauber model and provide a good description of the experimental data on quasi-one-dimensional ferroelectric crystals CsH

Figure 2 .
Figure 2. The dynamic dielectric permittivity of the crystal CsH 2 PO 4 .The lines present theoretical results (parameters in table1), the symbols show experimental data[25].

Table 1 .
The model parameters that describe experimental data on CsH 2 PO 4 and CsD 2 PO 4 : the nearest neighbour coupling K; the long-range interaction J; the effective dipole moment of a proton (deuteron) on a hydrogen bond µ; the parameter of the kinetic equation -relaxation time of the noninteracting proton (deuteron) α.
19))One can see that(4.19)isconsistent with the exact relations (4.18), though there is no guarantee that such a dependence among η, η 1 , η 2 and η 3 is preserved even in weakly nonequilibrium states.Using (4.19) and the condition 24))Substitution of (4.21) and (4.22) into (4.20) and linearization of the obtained system of equations with respect to η i , η it and E t allow one to obtain two systems of equations.The first one contains equilibrium correlation functions, and its solution yield relations (4.15) -(4.17).The second one contains a dynamic part of the correlation functions and has the form: