Investigation Of Locally Anharmonic Models Of Structural Phase Transitions. Seminumerical Approach

In the context of description of the structural phase transition in solids, the anharmonic potential eeecting on the ion in the crystal lattice is selected in the form that in addition to quadric includes cubic anharmonicity, thus possessing more general non-symmetrical form, and in the symmetrical form with the Gaussian barrier. The approach based on the numerical diagonalization of the one-site Hamiltonian matrix is applied. The peculiarity of the energy spectrum of the system with the Gaussian barrier potential, particularly the existence of the parts of the spectrum with the quasidoublets or quasi-equally distanced levels structure which crossover takes place in the narrow energy interval, is established. The dependencies of the order parameter , free energy and dielectric susceptibility on the externally applied eld and the order parameter and free energy on the temperature are derived for various parameters of both models. The phase diagram (temperature-eld) is derived and the innuence of the anharmonicity on the form of the phase coexistence curve and location of the critical point is studied. It is shown that the asymmetry of the local potential results in the deviation of the phase coexistence curve from the vertical line in the plane (temperature-eld) and in the appearing of the rst order phase transition occurring as the temperature varies.


Introduction
While studying the mechanism of the structural phase transitions in crystals, the problem of choosing the type of the local potential V (q) in the model Hamiltonian H = X i ( p 2 i 2M + V (q i )) ? 1 2 X ij ' ij q i q j ; (1) describing a certain mode of the lattice vibrations which is active at the phase transition, is raised.We use the concept of the local normal coordinate denoted here by q i (see, for example, 1]).
In the harmonic approximation (V (q) = 2 q 2 ; = M! 2 0 ), the crystal lattice may become unstable in the case of the negative value of .As it is mentioned in a number of papers, including 1], the simplest stabilizing interaction may be chosen in the form c I.V. Stasyuk V (q i ) = ?j j q 2 i 2 + q 4 i 4 : (2) In this case the system undergoes the structural phase transition which can be of the displacive or order-disorder type, depending on the values of parameters in the initial Hamiltonian.Such a model is often applied to the description of the phase transition in ferroelectrics.
A large number of papers was dedicated to studies of the thermodynamics and dynamics of the model described by the Hamiltonian (1), (2).Along with the application of the analytical methods of various degrees of approximation, including the method of self-consistent phonons, the method of renormalized group, etc., the attempts to take into account the e ects caused by the local anharmonicity of the model by using numerical methods were made.Particularly, in ref . 2] the energy spectrum of the single-ion system described by (2) was calculated, while the interaction between the ions in di erent cells was taken into account in the mean eld approximation; the comparison with the results obtained in the self-consistent phonons approximation was made.This approach is justi ed in the case of long-range interaction ' ij ; the obtained results can be improved by including of the higher order corrections when the expansion in powers of the inverse radius of interaction is used.
In addition to the description of ion motion in ferroelectrics, the model local anharmonic concept is applied to the description of the lattice anharmonicity in the high temperature superconductors.As it is pointed out in 3], such an anharmonicity is inherent to the motion of the apex oxygen ions in Y BaCuO and other superconductive compounds.To take into consideration such an anharmonicity in the case of local double-well potential, the approach which introduces pseudospin variables describing the vibrational degrees of freedom was used.The pseudospin-electron model derived in this way, describing also the interaction of the conducting electrons with local anharmonic mode, was the subject of intensive study in the last years [4][5][6].The asymmetry of the anharmonic potential which is characteristic for the systems of this kind was taken into account by including into the Hamiltonian the term which described the interaction with some internal eld hi , depending on the occupancy of electronic states ( hi = h + gn i ).Structural phase transitions in HTSC systems, including Hg-based superconductors, have been recently of wide consideration, in the context of the reported connection between the lattice softening in these compounds and transition to the superconducting state.Particularly, in ref . 7] it was mentioned that near the transition point an anomalous abrupt mode softening was observed.The peculiarities related to the presence of the lattice anharmonicity were observed for the Y BaCuO series on the temperature dependencies of the lattice constant c, coe cients of thermal expansion in the direction of anharmonic vibrations of oxygen ions, as well as speci c heat, thermoconductivity and velocity of ultrasonic waves in the form of jump-like and hysteresis behaviour (see 8]).
For the description of such situations, the model potential possessing more general non-symmetrical shape due to the presence of the cubic term, was considered in 8-9]: (so called "' 3 + ' 4 " model 1 ).Basing on (3) and applying the method of self-consistent phonons it was shown that the behaviour of the apex oxy-gen in Y Ba 2 Cu 3 O 7? compound is bistable.The dependence of the order parameter on the temperature was shown to be of a hysteresis character.However, the question of the applicability of the method of self-consistent phonons to the description of the phase transition in the systems described by (3), including critical areas arises, (see 1]).
Another possibility to describe the anharmonic potential acting on the ion in the crystal lattice is to write it as a superposition of a Gaussian and a harmonic oscillator potential: V G (q i ) = q 2 i 2 + Ce ?Bq 2 i ; (4) which has two symmetrically located minima if 2BC < 1.
This potential was applied, for example, in 11] and 12] to consider the crystal dynamics and phase transition in ferroelectrics within the selfconsistent phonon eld approximation.Since the role of the anharmonic interaction in (3) is not essential if the temperature is high, or when hq 2 i >> (1= ), the author of 6] identi es potential (3) as the "low-temperature anharmonicity potential".
The potential of V G (q i ) type has been used also in the theory of molecules.In 13] the numerical treatment of (4), similar to 2], was carried out; the low-lying energy levels were derived by solving the secular equations for various sets of potential parameters.The latter one were then derived for certain con gurations of some molecules.
The main aim of this paper is the further development of the numerical approaches mentioned above.We model the potential acting on the ion in the form (3) and ( 4) and consider the anharmonic character of the ion motion by the numerical diagonalization of the Hamiltonian matrix (1).The interaction between ions in di erent cells we allow for in the mean-eld approximation, assuming the long-range character of the intercell interaction.Basing on this approach, we derive the dependencies of the order parameter, free energy and dielectric susceptibility on the externally applied eld and of the order parameter and free energy on the temperature for various model parameters.Analysing these dependencies we investigate the phase transition as the external eld varies as well as the temperature; in the case of the potential in form (3) it is of the rst order; we also consider the behaviour of the dielectric susceptibility near the phase transition point.Phase transition in the system described by ( 4) is of the rst order, as the external eld varies, and second order, at the temperature change.In addition, we construct the phase diagrams (external eld, temperature) and study the in uence of the anharmonicity of each of potentials (3), (4) on the form of the phase diagrams and the location of the critical points.

Model Hamiltonian
In this section we consider the case of local anharmonic potential in nonsymmetric form (3), where = M! 2 0 > 0. Transformation of the last term in (1) in accordance with the mean-eld approximation leads to ? 1 2 X ij ' ij q i q j !'hqiq ? 1 2 'hqi 2 ; (5) This results in the following form of H: + dq i + 'hqiq i ? 1 2 'hqi 2 : (6) Here the external eld described by parameter d is introduced.On the basis of eigenfunctions of harmonic oscillator the operators p and q have the form of the following in nite-dimensional matrices: After substituting q and p in (6) according to (7), the Hamiltonian matrix has the following form (all terms are divided by ~!0 ): (8) where For the purpose of the numerical treatment of the Hamiltonian matrix (8) we limit the dimension of this matrix to the required size.As the calculation of the dependencies of thermodynamical functions on the model parameters shows, the limitation of the size of Hamiltonian matrix (this corresponds to allowing for the nite number of harmonic oscillator levels) to N = 25 is su cient when kT ~!0 6 5.All calculations in Section 2 are made within this approximation.

Order parameter and free energy
For calculation of the mean value of coordinate q i which has a meaning of the order parameter we use the expression hqi = Sp (qe ?H ) Sp e ?H : (11) After the unitary transformation is made H d = V ?1 HV; (12) which diagonalizes the Hamiltonian matrix (8), we get hqi = Sp (qe ?Hd ) Sp e ?Hd ; q = V ?1 qV: (13) Note that hqi is contained in the Hamiltonian (8).Denoting f = d + 'hqi, we derive the self-consistent system of equations: Numerical solution of this system allows to obtain the dependence hqi = hqi(d) ( gure 1a) for various values of kT h!0 (in all calculations here, C 1 =0.157,C 2 =0.025, and C 3 '=-20; these values of parameters correspond to the ones used in 5]).Substituting this function in the Hamiltonian (8), we can calculate the dependence of the free energy on the external eld, according to the expression below: F(d; T) = ?T ln Sp e ?Hd ? 1  2 'hqi 2 : ( Figure 1b represents the obtained result for F = F(d).C 1 = 0, the phase transition in the accepted approximation remains having a pure displacive nature.Thus, the peculiar feature of this model consists in the signi cant enhancement of the critical temperature due to the presence of the anharmonicity of third order.Figure 3 illustrates the dependence of the critical temperature on the value of cubic anharmonicity.? @qi @d . Te dependencies ?1 = ?1(d) are shown on gure 4a -4b.While the value of (or ?1 ) is nite at kT ~!0 = 1:4 ( gure 4a), which corresponds to the moving across the phase coexistence curve in the middle part of it (see diagram on gure 2a), it tends to in nity at kT ~!0 ' 3:9 (correspondingly, ?1 reaches zero), i.e. as it approaches the critical point at the phase diagram ( gure 4b).

Dependence of the order parameter and free energy on the temperature
Solving numerically system (14) for certain values of d we derive the dependence hqi = hqi( kT ~!0 ).Using this dependence in the Hamiltonian (8) and diagonalizing it we obtain the free energy as a function of the temperature, according to (15).The obtained dependencies of the order parameter and free energy are shown in Figs.5-7  as the temperature varies (on moving across the phase coexistence curve), the minimal free energy corresponds to the jump of the order parameter from branch 3 to branch 1 ( gure 6).At C 1 = 0:157 this takes place at the value kT ~!0 = 2:5, which corresponds to the rst order phase transition (this transition may be interpreted as a point of lattice bistability 9]).When d < 1:293 (the left part of the phase diagram), the order parameter varies smoothly along branch 1 in gure 7; no phase transition occurs in this case.
There is an another way which allows to determine the range of values We plot the dependencies of the free energy on the order parameter for d = 1:4 and d = 1:2 and kT ~!0 = 0:2 in Figs.9a -9b.Comparing the dependence in gure 9a and the dependence of the order parameter on the temperature shown in gure 5a (case d = 1:4) one concludes that the lowest minimal value of the free energy (at hqi ' ?0:07) corresponds to branch 3 of the gure 5a.The same is true for all values of kT ~!0 , therefore, the order parameter varies smoothly along this branch.The same smooth behaviour of the order parameter in the case of d = 1:2 can be concluded by comparing gure 9b and gure 7a; in this case the lowest minimal value of the free energy in the gure 9b (at hqi ' 0:2) corresponds to the order parameter varying along branch 1 of the gure 7a.
The case of d = 1:3 for two values of kT ~!0 = 2 and kT ~!0 = 3 is represented on Figs.9c -9d.At kT ~!0 = 2 the minimal value of the free energy in gure 9c corresponds to the value of the order parameter hqi ' 0:18, or branch 1 of the gure 6a.However, at kT ~!0 = 3 the minimal value of the free energy in gure 9d corresponds to the value of hqi ' ?0:02, or branch 3 of the gure 6a.Between these two cases, at kT ~!0 = 2:5, a jump-like transition of the order parameter from branch 1 to branch 3 occurs ( gure 6a).This con rms the conclusion resulting from the analysis of the dependencies of the order parameter and free energy on the temperature.
We note that in the analysis presented above we did not raise the issue of how the possibility of the occurrence of the phase transition depends on the constant of the long-range interaction '.One can see that as well as in the case of pure quadric model (C 1 = 0), the phase transition disappears if ' is less than the certain critical value (j ' j<j ' j).The calculation of the wider class of the phase diagrams, taking into account the varying of ', should be the subject of a separate investigation.

Model with the Gaussian barrier-type potential 3.1. Model Hamiltonian and spectra
While in the previous sections the potential considered possessed non-symmetrical form due to the presence of the cubic anharmonicity, we consider now the potential of the symmetrical form (4): V G (q i ) = M! 2 0 q 2 i 2 + Ce ?Bq 2 i : (16) We use the same approach as we did in the previous sections, that is we construct the Hamiltonian matrix H mn on the basis of the normalized eigenfunctions of the harmonic oscillator (with the potential energy given by the rst term in ( 16)), using relations ( 1) and ( 5), with V G (q i ) in the form above. Then e ?(1+BM! 0 )x 2 H n (x)H m (x)dx: Here H n (x) are the Hermite polynomials, and q mn are de ned in (7).
Below are the relations used for calculation of the integral appearing in (19): e ?px 2 H 2m (x)H 2n+1 (x)dx = 0: Using relations ( 16) and ( 18), we write V G (q) C = q 2 + e ?Bq 2 : (20) The condition of VG(q)   C to have two minima is < B. Assuming in the calculations below B = 1, this condition rewrites as < 1.Then the meaning of is that the closer is the value of to 1 the smaller is the height of the Gaussian barrier (if = 1, the potential (4) has only one minimumthe barrier disappears).
Transforming the matrix H to the diagonal form in the absence of the external eld d, we derive the spectra for the various values of .Approximately 35 lowest levels of the spectra are plotted on gure 10 at = 0:4 and = 0:7; 32 energy values ~!0 at = 0:4 and = 0:7 are also given in tables 1 and 2, respectively.(We note that within the range of the model parameters considered, particularly when = 0:4 ?0:7, the number of harmonic oscillator levels which were needed to be taken into account to provide the required accuracy ( " " = 7 10 ?5 ) of calculations of this part of spectrum or the required size of the Hamiltonian matrix was N ' 100).Approximately ten lowest levels of the spectrum on gure 10a ( = 0:4) are doublets (which is a well-known consequence of the presence of doublewell potential), although the rst three doublets can not be visually resolved (for example, the width of the rst doublet " 1 is much smaller than the distance between the rst and the second doublet " 12 : "1 "12 ' 4 10 ?4 , see table 1).As one moves to the middle and upper parts of the spectrum, one can see that the spectrum gradually transforms into the system of equally distanced energy levels which is typical for the harmonic oscillator, thus proving that the potential in form ( 4) is the low-temperature anharmonicity potential indeed.If we decrease the height of the Gaussian barrier ( gure 10b, = 0:7), the number of doublets in the spectrum decreases.The crossover between doublet levels and nearly equally distanced parts of the spectrum takes place in this case at the closer to ground state values of energy.Finally, if > 1, or if the potential (20) has only one minimum, all spectrum becomes more and more closer to the system of equally distanced energy levels.

Order parameter and free energy. Phase diagram
Applying the procedure described in Section 2.2 (relations (11)(12)(13)(14)(15)) to the Hamiltonian (17), we derive the dependencies of the order parameter and free energy on the external eld which are shown on gure 11.Unlikely to the case of non-symmetrical potential (3), where the value d = d at which the jump-like change of the order parameter occurs was di erent from zero (see gure 1 and gure 2), the phase transition takes place now at d = 0.This is the result of the symmetrical form of potential (4).The temperature change does not in uence the value of d .Therefore, the phase coexistence curve shown on gure 12 is the vertical line.Its end's position (the critical point) depends on the value of .Figure 13 illustrates the connection between the critical temperature kT ~!0 and the values of parameter .As it follows from gure 13, for the smaller value of , i.e., the larger height of the Gaussian barrier, the larger value of the critical temperature kT ~!0 is needed to destroy the hysteresis-type dependence of the order parameter on d, or to destroy the phase transition.We note that even when > 1, i.e. when the potential (4) transforms from two-minima to the oneminimum potential, the system still undergoes a phase transition (which in this case is displacive in character, while it is of order-disorder type in case of < 1, when the potential has two minima): as it follows from gure 13, the phase transition in the system occurs until ' 4.

Dependence of the order parameter and free energy on the temperature
The dependencies of the order parameter and the free energy on the temperature are derived by the calculations similar to those in Section 2.3.We We also plot the dependencies F = F(hqi) for d = 0:1 and d = 0 in gure 16.For d = 0:1 the minimal value of F corresponds to branch 3 of hqi on gure 14a (hqi ' ?0:54), thus showing that hqi varies along this branch; for d = 0 F = F(hqi) has two equivalent minima at hqi ' 0:54, which correspond to branches 1 and 3 on gure 15a.This is in agreement with the above given conclusion and con rms that in this approximation the phase transition at the temperature change is of the second order.The self-consistent phonons approximation results, however, in the phase transition of the rst order 11].The reason for such a discrepancy is the di erent character of both approaches used: the rst approach does not fully takes into account the long-range interaction, and the second one treats approximately the local anharmonicity.More consistent description of the phase transition in the framework of the proposed in this paper scheme certainly requires the consideration of quantum uctuations relatively to the mean eld.

Conclusions
The approach based on the numerical solution of the single ion problem in the crystal lattice is developed for the investigation of the local anhar--0.10-0.05 0.00 0.05 0.10 0  3) which e ects on the ions, as well as the symmetrical potential of a special form (4).The long-range interactions are taken into account in the mean-eld approximation.The dependencies of the dynamic and thermodynamic functions on the external eld and the temperature are derived.In the case when the system is described by the potential (3), these dependencies testify that the system undergoes the rst order phase transition, with the jump-like change of the order parameter hqi, as the external eld varies as well as the temperature.It is established that the increase of the cubic anharmonicity constant leads to the increase of the critical temperature.In the case when the system is described by the potential (4), it undergoes the rst order phase transition as the external eld varies and the second order phase transition as the temperature varies at zero value of the eld.The increase of the Gaussian barrier height results in the increase of the critical temperature.The phase diagrams (d; T) of the models described by (3) and (4) are built.
We would like to point out that the results obtained from the consideration of the problem with non-symmetrical potential (3) regarding the character of the phase transitions generally comply with the conclusions of 6] and 14], in which the similar investigations were related to the pseudospinelectron model of HTSC systems.Both in the those cases, and in the model (3), the potential asymmetry leads to the deviation of the phase coexistence curve from the vertical line on the plane (temperature-eld) and to the appearance of the rst order phase transitions (which may manifest themselves in the instability of the dielectric susceptibility 6]) as the temperature varies in the certain range of the parameter which de nes mentioned asymmetry.Existence of this phase transition corresponds to the bistability of the system with the local potential (3) mentioned in 9].However, the character of the dependence hqi = hqi( kT ~!0 ) (shown on gure 6) di ers from the one derived in 9] which is the result of using di erent approximations.
The substantial conclusion can be made from the investigation of the single-ion spectra made in Section 3: the existence of two parts of the spec-  trum of di erent character, with quasidoublet (at the lower values of energy) and quasi equally distanced (at the higher values of energy) structures.The crossover from one part to another is relatively sharp and takes place in the range of energies which are of the height of potential barrier by the order of magnitude (the possibility of such a transformation in the spectrum was foreseen in 15] and was formulated in terms of simple model).As it follows from the results obtained, the model (4) incorporates phase transitions of displacive (1 < < 4), as well as order-disorder types ( < 1).
The approach used in this paper may be applied to the investigation of the dynamics of the considered models and thermodynamics of the phase transitions using approximations higher than the mean eld approximation.In doing it, the capabilities of this seminumerical approach can be demonstrated and its comparison with approximate analytical schemes can be made.

Figure 1 .Figure 2 .
Figure 1.Dependencies of the order parameter and free energy on the external eld, kT ~!0 =0.17.The dependencies shown in gure 1 are typical for the rst order phase transition.The abscissa of the self-crossing point of the curve F corresponds to the value d at which the phase transition occurs, causing a jump-like change of the order parameter on gure 1a.Raising the temperature leads

Figure 3 .
Figure 3. Dependence of the critical temperature on the value of cubic anharmonicity On the basis of the obtained dependence hqi = hqi(d), one can calculate

Figure 4 .
Figure 4. Inverse dielectric susceptibility as function of the external eld, for the values of kT ~!0 =1.4 and kT ~!0 =3.9 which correspond to the middle part and to the left end of the phase coexistence curve shown in gure 2a, respectively.

Figure 9 .
Figure 9. Dependencies of the free energy on the order parameter: d = 1:4 (a), d = 1:2 (b), ( kT ~!0 = 0:2), and d = 1:3, kT ~!0 = 2 (c), kT ~!0 = 3 (d), C 1 = 0:157.of d where the phase transition occurs.We can calculate the dependence of the free energy on the order parameter, according to (15), for the xed values of kT ~!0 , and di erent values of d, which fall, respectively, to the right of, to the left of, and in the range of the phase coexistence curve.We plot the dependencies of the free energy on the order parameter

Figure 11 .
Figure 11.Dependence of the order parameter and free energy on the external eld, kT ~!0 = 0:2, = 0:4.plot dependencies of hqi = hqi( kT h!0 ) and F = F( kT h!0 ) for d = 0:1 and d = 0 (see Figs. 14 -15).In the range d > 0 ( gure 14) branch (3) of hqi = hqi( kT ~!0 ) corresponds to the minimal value of the free energy which therefore means that hqi varies smoothly along this branch.The same is true with interchanging of roles of 1 and 3 branches in the case of d < 0. These two cases (d 6 = 0) are similar to those represented in Figs. 5 and 7.The Case of d = 0, which corresponds to the phase diagram of gure 12, is presented on gure 15.Branches 1 and 3 of hqi = hqi( kT ~!0 ) ( gure 15a) correspond to one "degenerated" branch of F = F( kT ~!0 ) ( gure 15b) and therefore are equivalent in terms of the possible values of hqi (see gure 15a).At the branching point ( kT ~!0 ' 10:4) the appearance of the non-zero value of hqi corresponds to the second order phase transition.We also plot the dependencies F = F(hqi) for d = 0:1 and d = 0 in gure 16.For d = 0:1 the minimal value of F corresponds to branch 3 of hqi on gure 14a (hqi ' ?0:54), thus showing that hqi varies along this branch; for d = 0 F = F(hqi) has two equivalent minima at hqi ' 0:54,

Figure 12 .
Figure 12.Phase diagram for two di erent values of : = 0:5 and = 2:5: monicity model of structural phase transitions in solids for the more general non-symmetrical form of the potential (3) which e ects on the ions, as well as the symmetrical potential of a special form (4).The long-range interactions are taken into account in the mean-eld approximation.The dependencies of the dynamic and thermodynamic functions on the external eld and the temperature are derived.In the case when the system is described by the potential (3), these dependencies testify that the system undergoes the rst order phase transition, with the jump-like change of the order parameter

Figure 13 .
Figure 13.Dependence of the critical temperature on :