Proton ordering model of phase transitions in hydrogen bonded ferrielectric type systems : the GPI crystal ∗

A microscopic model based on the consideration of the proton ordering is proposed for describing the H-bonded ferroelectric crystalline systems with a complex structure of the hydrogen bond network. The model has been used for the investigation of thermodynamics and dielectric properties of the GPI crystal. The symmetry analysis of the order parameters responsible for the mixed (ferroand antiferroelectric) nature of ordering is performed within the model. The phase transition into the ferroelectric state is described. Changes in the dielectric susceptibility of the crystal are studied in the presence of the transverse external electric field acting along the c axis. The results of measurements of temperature and field dependences of dielectric permittivity ε c in the paraelectric phase are presented. The microscopic mechanism of the observed effects is discussed based on the comparison of theoretical results and experimental data. A conclusion is made about the significant role of the ionic groups connected by hydrogen bonds in the charge transfer. So they make an important contribution into the polarizability of the GPI crystal along the direction of H-bonded chains.


Introduction
Ferroelectric crystals with hydrogen bonds have been the subject of investigations for many years.The phase transitions in such systems take place at the decrease of temperature mainly to the ferroelectric or antiferroelectric phase and are connected with the ordering of protons in hydrogen bonds (the proton on a bond is in a double potential well).Due to the features of the lattice structure, spontaneous polarization can appear in the direction which not necessarily coincides with the orientation of hydrogen bonds.Such an effect can be caused by participation of the other elements of structure (ionic groups) in the creation of spontaneous dipole moments (as it takes place, for example, in the KDP crystal).At the same time, the situation is possible when the hydrogen bond network is complex and the proton ordering produces the dipole moments diversely oriented in various sublattices.It corresponds to the ferrielectric state, where the total dipole moment is different from zero and arises due to partial decompensation of certain components of separate dipole moments.Dielectric properties of such objects are of interest because of the existing internal correlation between contributions of protons into the crystal polarization in different directions.One should expect the manifestation of peculiar effects connected with the external electric field action.Similar to the field effect in magnetic systems, the electric field can lead here to the reorientation of dipole moments and can induce the corresponding phase transitions.
Among the H-bonded ferroelectric type systems, in which the above mentioned situation takes place, one can consider the glycine phosphite, NH 3 CH 2 COOH 3 PO 3 (abbreviated as GPI) crystal belonging to the group of recently synthesized crystals with the hydrogen bond network possessing a rather complicated structure.
Their special feature is the presence of a chain-like arrangement of hydrogen bonds sequentially linking the phosphite H(HPO 3 ) − anions.Besides, the glycinium cations [H 3 NCH 2 COOH] + alternately interconnect the chains by the hydrogen bonds.At room temperature the crystal is monoclinic (the space group P2 1 /a).The unit cell with dimensions a = 9.792 Å, b = 8.487 Å, c = 7.411 Å, β = 100.43• contains four molecular units [1].The hydrogen bonds linking the phosphite groups are symmetric (with hydrogens equally disordered over two sites).Investigation of symmetry and structure of crystal in the ferroelectric phase (below temperature of the second order phase transition T c = 224.7 K) was recently performed by means of X-ray structure analysis [2] and neutron diffraction study [3].The space group of this phase is P2 1 .Spontaneous polarization appears along the b-axis and is perpendicular to the direction of the chains of hydrogen bonds (the c-axis).Protons are ordered at positions approximately equal to one of the two previous sites.Phosphite tetrahedra are distorted with respect to their shape in the paraelectric phase.
Physical properties of GPI crystal were investigated experimentally in a number of papers.The effect of deuteration on T c was studied.The effect is significant: the phase transition point shifts considerably to higher temperatures (T D c − T H c = 97 K [4]).The dependence of T c against deuterium contents is linear [5].The measurements of dielectric constant along the b-axis were performed.The temperature behavior of the ε b component as well as the spontaneous polarization P s were studied; the conclusion was made that the phase transition can be regarded as the one close to a tricritical [6].The absolute value of P s in the state of saturation (∼ 5 × 10 −3 C/m 2 [7]) is considerably smaller in comparison to other crystals of hydrogen bonded fam-ily.A large anomaly of permittivity in the c-direction has also been detected; at room temperature ε c is higher than ε b [7].
Apart from the structure investigations, the measurements of frequency dependence of dielectric constant also lead to the conclusion that the phase transition in GPI is of the order-disorder type [8,9].However, the observed dynamics of the soft relaxational mode connected with the proton flipping in the double potential well is more complicated in GPI than in the other similar H-bonded ferroelectrics (for example, BPI).The estimates of the corresponding activation energy and relaxation time imply that para-ferroelectric phase transition in GPI does not seem to be strictly connected only with the movement of protons [9,10].A reorientation and deformation of the ionic group (phosphite anions) could also play an important role (the corresponding structure changes can be extracted from the neutron scattering data [3]).
On the other hand, the results of experiments with deuterated crystals (DGPI) give a strong support to the proton ordering mechanism of the phase transition and point to the important role of hydrogen bonds.Besides the large shift of T c , one can mention the "geometric" isotope effect (a linear change of the R 0−0 distance in hydrogen bond at the deuteration of crystal).There is a direct correlation between ∆R 0−0 and δT c for partially deuterated GPI crystals [11], which is very similar to the one for a variety of H-bonded ferroelectrics with phase transitions caused by the ordering of protons.
One should also mention the high hydrostatic pressure studies of the transition into ferroelectric phase in GPI.The effect of pressure is found to be significant.With increase of pressure, T c decreases with dT c /dp = −11.0K/kbar [12].Such a decrease is comparable with the value dT c /dp = −11.8K/kbar reported for the BPI crystal [13], but is larger than for KDP-type ferroelectrics.The hydrostatic pressure effects also the saturation value of spontaneous polarization.Sensitivity of GPI crystal to the external stress was also studied by means of ultrasonic investigations of elastic constants [14].The high values of linear compressibility β were obtained along the a and c-directions (β ∼ 0.16 × 10 −11 m 2 /N ).It is a consequence of a weak coupling due to the layered arrangement along the a-axis and the hydrogen bonding of ionic groups along the c-axis.The observed temperature anomalies of elastic constants in the vicinity of the T c [15] were explained as a manifestation of the probably pseudoproper origin of ferroelectricity in GPI due to anharmonic interaction of the order parameter with deformation of the crystal lattice.
Investigations of the IR and Raman spectra [10], the thermal expansion coefficients, spontaneous electrostriction and refraction indices measurements were also performed [16].Relatively large values of the thermal expansion coefficients, especially in the c-axis direction, give some more evidence of the existence of the noticeable lattice anharmonicity.
Although from the present experimental data the conclusion can be made about a significant role of hydrogen bonds and ordering of protons in the appearance of ferroelectric phase in the GPI crystal, an appropriate microscopic model of transition into this state has not been considered so far.Such a model is proposed in this work due to the information obtained last year about the changes in the structure (including the changes in distribution of protons on hydrogen bonds) at the paraelectric-ferroelectric phase transition [2,3].The symmetry analysis of the order parameters is performed.The phase transition into ferroelectric phase with P s b at the presence of the external electric field acting along the c-direction is described based on the model.Attention is paid to the shift of the transition temperature and to the changes of dielectric constants in the c-direction.The results of measurements of temperature and field dependences of dielectric permittivity ε c in the paraelectric phase are presented as well.The predictions of theoretical investigations are compared with the experimental data.Microscopic mechanism of the observed effects is discussed.

Proton ordering model
In our approach the description of the proton subsystem of the GPI crystal will be based on consideration of the long-range and short-range interactions of protons.We start from the structural data for localization of hydrogen bonds and equilibrium positions of protons in paraelectric and ferroelectric phases [2,3].There are four translationally nonequivalent bonds in the unit cell of crystal (the 011−021 and 013−023 bonds as well as other two bonds connected with the mentioned ones by the point symmetry operations).Their projections on the XZ and Y Z planes are shown in figure 1 (Cartesian Y and Z-axes coincide with crystallographic axes b and c, respectively, the cartesian X (or a )-axis makes an angle of 100.43 • with the crystallographic axis c in the (a, c) plane).Directions connected with the ordering of protons and corresponding to the prevailing occupancy of proton positions (n ia > n ib ) are indicated by arrows.
As can be seen, the dipole moments d i arising due to the displacement of protons along the bonds compensate mutually in pairs (d 1 and d 3 ; d 2 and d 4 ) in the low temperature phase in Z and X-directions.At the same time, they sum up along the Y -axis leading to the appearance of the spontaneous polarization.Due to the zigzag structure of the hydrogen bond network we have the ferrielectric type ordering, when the ferroelectricity along b-axis is accompanied by the antiferroelectric-like (antiparallel) arrangement of dipole moments of the neighbouring H-bonded chains along the c-axis (the chains are oriented in this direction).
According to the data given in [2], in the fully ordered state The numerical coefficients are given in Å; q is an effective charge transferred along the bond at the displacement of proton from one equilibrium position into another.
To describe the proton ordering we use the pseudospin variables where n iα is the occupation number of proton in the α position on the bond i, and introduce the dipole moment operators Dα i = d α i S z i .The pseudospin mean values η i = S z i are equal to zero in the paraphase while in the ferroelectric phase The energy of proton subsystem can be written in this representation in the following general form: where n is the unit cell number, k is the number of sublattice, α = x, y, z.The effect of an external electric field E is taken into account.The short-range interactions between protons belonging to the neighbouring bonds in the H-bonded chains are included in the expression (2.3):For a description of the long-range proton-proton interactions we use the mean field approximation.In this approach the corresponding part of the Hamiltonian and the fields H µ 1,2 acting on protons on the bonds s = 1, 2 are introduced where j kk are the constants of the long-range interactions.The self-consistent set of equations for the η i mean values is obtained here within the scheme based on exact consideration of short-range interactions.Such an interaction can be quite large in GPI, since redistribution of protons on the bonds is accompanied by deformation and reorientation of ionic groups connected by the bonds [3].The procedure applied below is well known and accepted for the quasione-dimensional systems (see [17,18]) and corresponds to the two-particle cluster approximation [19,20].

Phase transition in proton subsystem
We use the transfer-matrix method for the calculation of partition function and other thermodynamic functions of H-bonded chains in GPI in the framework of the model considered.In such a case a partition function of the proton subsystem can be written as where λ a,(b) is the maximum value of the transfer matrix for the a(b) chain and β = 1/Θ = 1/k B T .Respectively, the free energy per a unit cell is The self-consistency parameters η i are determined from the minimization procedure which reduces in our case to the following set of equations These equations describe the continuous transition (of the second order) to the phase with nonzero values of the η 1 and η 2 parameters and spontaneous polarization along b-axis.In the presence of the external field E z , the phase transition remains and the order of this transition does not change (the field E z acts along the Hbonded chains perpendicularly to the ferroelectric axis).Such a conclusion follows from the analysis of the free energy behaviour in the vicinity of a transition point.
The temperature of phase transition T c is determined from the equation which is obtained from the condition of the appearance of non-zero solutions of equation (3.5).Here where a ik are the parameters formed by the long-range interaction constants.It can be seen that under the action of the field E z the transition temperature decreases.The shift z at small values of E z .The effect takes place also in the absence of the short-range correlations.In particular, in the case when ϕ where k 1 and k 2 are the constants connected with linear (with respect to field E z ) contributions in the internal fields H µ 1,2 ; T 0 c = T c | Ez=0 .Respective changes take place in the behaviour of dielectric susceptibility χ yy along the ferroelectric axis.In the model, the contribution of the proton subsystem is given by the expression where ε 0 is the electric permittivity of vacuum, v c is the volume of the unit cell.
As it follows from equations (3.5) at E y → 0, E z = 0, in the case of paraphase (T > T c ) 3 )be] .(3.12) In the external field the susceptibility χ yy diverges at T c while at T > T c it decreases proportionally to E 2 z under a small enough field.This effect can be interpreted as a result of the shift of T c to the lower temperature side.
There are also changes in the dielectric susceptibility χ zz ≡ χ c under the effect of the field E z .In the model, it is given by the expression formally obtained from (3.12) by substitution of parameters ϕ i .Calculations show that χ zz also decreases proportionally to E 2 z at T > T c .The reason for such a behaviour may appear to be more complicated than in the case of the χ yy susceptibility (see discussion below).The above mentioned can be illustrated by the results of numerical calculations.In figure 2 the field dependence of T c is shown calculated according to equation (3.6).Respectively, in figure 3 the temperature dependences of χ zz at different values of the field E z are presented.
In our numerical calculations and in figures 2 and 3 all energetic quantities were taken in relation to j + j (the sum of the short-range interaction constants) and dipole moments were given in relation to µ 0 = q √ a 2 + b 2 + c 2 ; respectively, the dimensionless ratios were introduced.For example The long-range interaction has been modelled by two parameters Φ and Φ responsible for the interaction along a separate chain and between chains, respectively: where jik = j ik µ 2 0 /(j+j ).The values Φ = 0.053, Φ = −0.031were taken.This choice, as well as the values j + j = 0.204 eV, T 0 c = 0.094, µ 0 = 7 × 10 −30 C•m, q = 0.17e1 , correspond to the values of the observed critical temperature (T 0 c = 224.7 K [7]) and saturated spontaneous polarization (P s = 5 × 10 −3 C/m 2 [7]).
In this case, for susceptibility we have χzz = 25−35 (figure 3b) at T > T c , which corresponds to the change of the true susceptibility from 7 to 10 (the coefficient µ 2 0 /ε 0 v c (j + j ) is equal to 0.28 at the chosen parameter values).To achieve an agreement with the experimental data (χ zz ≈ ε c ≈ 200 at T = 273 K) we have to suppose that the z-components of the effective dipole moments, arising due to the displacements of protons along the bonds under the field effect, are approximately 5-6 times larger than the initially taken values.The necessary result is obtained at the increase of the c and f parameters.In this case, the dimensionless value Ẽz = 1.5 × 10 −3 corresponds to the field E z = 1 MV/m.The difference ∆ T = 0.002 at Ẽz = 0.006 corresponds to the shift ∆T c = 4.5 K at the field E z = 4.0 MV/m (see the experimental data below).
It should be mentioned that in the proton ordering model, which reduces in the considered case to the four-sublattice Ising model, the phase transition into the ordered phase is of the second order.It concerns both the transition into ferroelectric phase (phase F ) with polarization along b-axis and the transition into the ordered phase (hypothetical phase F ) with the orientation of dipole moments parallel to caxis (and, respectively, antiparallel orientation with respect to b-axis).Such a phase is not realized in GPI, but the example of phase transition with the substitution of ferroelectric and antiferroelectric orderings is known in the literature; in this connection one can mention the DMAGaS crystal [22].
Because of the mentioned property of the model the dielectric susceptibility χ yy decreases in paraphase under the field E z .The orientation of the dipole moments parallel to c-axis (induced by the field E z ), counteracts the appearance of ferroelectric phase with P s b.It leads to the decrease of the critical temperature T c and, as it was mentioned above, to the effective decrease of the χ yy component.On the other hand, the ordering induced along c-axis suppresses the polarization fluctuations in this direction which results in the going down of the χ zz component.Large enough values of the component observed experimentally at T T c in the wide temperature range and the increase of χ zz at the decrease of T can be connected with the proximity to the point of the virtual phase transition to the phase F (the latter could be realized in the GPI crystal in the absence of transition to the phase F at T c ).
In addition, a certain role in the formation of the high polarizability of the crystal in the direction of the H-bonded chains belongs to the short-range proton interactions.This can be seen considering the expression for "nonperturbed" susceptibilities in the paraphase (when the long range interaction is neglected) Here, only the interaction between the neighbouring hydrogen bonds along a chain is taken into account.The ratio χ zz /χ yy changes with temperature approaching to maximum values at low temperatures; formally At high temperatures it goes to the value At T 0.12(j+j ), which corresponds to the condition T 285 K, χ zz /χ yy ∼ 3.0 (see figure 4, χ zz /χ yy = χzz / χyy ).We see that due to the considered mechanism the χ zz component of susceptibility can about three times exceed the χ yy component.The observed ratio of the ε c and ε b permittivities at room temperature is much larger (ε c /ε b ∼ 10 [6,7]).The long-range interaction, which favours the proton arrangement characteristic of the hypothetic F phase, can be the additional reason for such an anomaly.

Experimental
Hereinbelow we present the results of the experimental investigations of the temperature dependences of the dielectric permittivity ε c = 1 + χ zz in the presence of the electric field E z applied in the direction of the c-axis.In the absence of external field the temperature behaviour of ε c was measured earlier [7].
Single crystals of GPI were grown from saturated water solution of deuterated polycrystals by slow evaporation method at temperature of 304 K.The dielectric properties were investigated along the crystallographic c-axis.The capacity of the samples was measured with a precision LCR-meter HP 4284 A at the measuring field frequency equal to 10 kHz and amplitude of 4 V/cm.The relative electric permittivity was measured as a function of temperature with various values of dc electric field ranging from zero to 4 × 10 6 V/m.The rate of temperature changes was equal to 0.5 K/min The temperature dependences of relative electric permittivity ε c for various values of dc field obtained from dielectric studies are presented in figure 5.At the applied field equal to zero, a large value of permittivity is observed at room tem-  perature and its increase is observed on the cooling run.The maximum value is noted at T c while below this temperature the permittivity decreases.Temperature dependences of permittivity under the electric field are similar to that observed at zero field.However, as it can bee seen in figure 5, in the vicinity of phase transition the additional increase of permittivity is observed.The increase becomes higher with the increase of electric field intensity.Besides, its position is shifted to lower temperatures with an increase of the electric field intensity.
The results obtained from figure 5 are presented in the form of temperature dependences of inverse permittivity in figure 6.The fulfillment of the Curie-Weiss law is seen both in paraphase and below the phase transition region in the ferroelectric phase.
The change of permittivity under the field in the vicinity of T c has a form of the jump-like increase.We can introduce, respectively, a jump of 1/ε c (defined as ∆ 1 , see figure 7) in the middle point of the interval at which the transition from the straight line for 1/ε c in paraphase to the straight line representing the strong increase of ε c in ferroelectric phase takes place.Without the field, the jump ∆ 1 is equal to zero.Assuming internal bias field E b existence and introducing the field (±E − E b ) we can notice that the value of the ∆ 1 increases linearly with (±E − E b ) 2 as it follows from the experimental data (E b value is small in comparison with the value of the applied fields).
As a temperature T c of the phase transition we assume the temperature that corresponds to the jump ∆ 1 .The position of the jump changes with the field E c .The shift of phase transition temperature obtained in such a way is presented in figure 8.The shift is satisfactorily described by the parabolic function with the maximum at the temperature T 0 c which is the point of the phase transition into ferroelectric phase at the zero field (the slight asymmetry of the T c versus E z curve is caused by the presence of the small bias field).Independent measurements of the   field dependences of ε c at fixed temperatures (T T 0 c ) reveal the anomalies similar to the described above, at the field values which correspond to the curve in figure 8.This fact confirms the existence of the paraelectric-ferroelectric transition at E z = 0 and points to the possibility of transformation from one phase to another at change of temperature or field E z .
The field dependence of the relative permittivity ε c measured in the paraphase (T = 230.9K) is illustrated by the curve given in figure 9.It was obtained using the procedure of eliminating the hysteresis behaviour of the ε c , caused by the internal field changing its value with the same frequency as the frequency of the applied field but shifted in phase.The effective field introduced as a difference of both fields (E eff = E − E corr ) is used, and in that case the maximum value of permittivity is reached at zero field.As a result, a quadratic dependence ε c (E c ) is obtained; it describes decrease of ε c (and, respectively, the increase of ε c −1 ) with the field at T > T c .Such a behaviour of ε c is characteristic for the whole temperature range in paraelectric phase.

Discussion
The experimentally obtained and theoretically calculated temperature dependences of direct (ε c ) and inverse (ε c −1 ) permittivity as well as their changes in the presence of the field E z are in a good mutual agreement.In particular, the described within the model decrease of ε c (increase of ε c −1 ) with the field in paraphase is in a direct correspondence with the experimental curve in figure 9.The changes of χ zz at the field values from figure 3 (∆χ zz ≈ 1.3 at E z = 1.0 MV/m; ∆χ zz ≈ 5.1 at E z = 2.0 MV/m; T = 239 K) agree with the quadratic dependence of χ zz on E z .
As it was mentioned, the quantitative agreement of the calculated values of χ zz with the data of measurements can be obtained by increase of the effective dipole moments along the c-axis connected with the proton displacements.It is a manifestation of an existence of a certain additional way of the charge transfer (and, respectively, polarization), caused by the redistribution of protons between equilibrium positions on bonds (an effective charge increases: q = 0.17e → q * = (0.85−1.00)e).Most probably, such a mechanism consists in the participation of the intermediate ionic groups (phosphite ions), which are shifted and distorted at the proton ordering.Their contribution can be taken into account implicitly in our model by the renormalization of the corresponding components of dipole moments connected with hydrogen bonds.
The change of T c under the field E z calculated within the proton ordering model is qualitatively the same as the measured one (figure 8).In both cases the shift ∆T c is described by the parabolic function of the field.As it was mentioned above, we can achieve a quantitative agreement between them using the optimal choice of numerical values of the interaction constants in the model.

Figure 1 .
Figure 1.Equilibrium position of O-and H-atoms on hydrogen bonds in the unit cell of GPI crystal in the paraphase (the directions of dipole moments connected with the ordering of protons in the ferroelectric phase are shown by arrows): (a) projection on the XZ (a c) plane; (b) projection on the Y Z (bc) plane.

. 4 )
Here and below the bonds k = 1, . . ., 4 in the unit cell n are characterized by the chain indices µ = a, b and, for the given chain with the index (l, m), by the cell number i and the number of bond s = 1, 2 in the cell (l = 1, . . ., M , i = 1, . . ., L, LM = N ).

)Figure 2 .
Figure 2. The field dependence of the paraelectric-ferroelectric phase transition temperature calculated based on the proton ordering model (2.4) and (2.5).

Figure 3 .
Figure 3.The temperature dependence of (a) longitudinal χ yy and (b) transverse χ zz susceptibility of GPI crystal within the proton ordering model.

Figure 4 .
Figure 4.The ratio χzz / χyy as a function of temperature in the short-range interaction approximation.

3 )Figure 5 .
Figure 5.The temperature dependence of dielectric permittivity ε c measured at different values of dc electric field E z .

Figure 6 .
Figure 6.The temperature dependences of inverse dielectric permittivity ε c −1 at different values of dc electric field E z (experiment).

Figure 7 .
Figure 7.The jump of inverse dielectric permittivity ε c −1 in the point of phase transition.

Figure 8 .
Figure 8.The field dependence of the phase transition point, taken from the measurements at fixed values of electric field E z or temperature T .