On collective variables method in the microscopic theory of alloys

Ab initio approach is developed for thermodynamical investigations of disordered binary alloys. It is based on utilizing the collective variables method. The explicit expression for the free energy and the equation determining the alloy components’ chemical potentials are obtained within the rigid lattice approximation and for the case when atomic static displacements (ASD) are taken into consideration. An ASD drastic effect on the behaviour of the binary correlation function Fourier components in the first Brillouin zone is observed. The ASD is shown to favour the ordering tendency in alloys of Ca-Ba and K-Cs systems. An approach in which configurational and vibrational degrees of freedom are considered at the microscopic level in the grand partition sum calculation is advanced. The role of the atomic thermal vibrations in alloy forming is discussed. The condition when configurational and vibrational effects can be treated separately is formulated.


Introduction
The method of the collective variables (CV) was introduced in [1] and developed by I.Yukhnovskii [2].The efficiency of this method was demonstrated during the last three decades in solving quite different problems of condensed matter physics: theory of liquids [3][4][5], statistical physics [6] and theory of multicomponent systems [7][8][9].The modified CV method (named as displacements and collective variables method) was successfully used for describing the interacting electron gas [10][11][12] and bose-systems [13,14].
The theory of the second order phase transition [6,15] and the study of liquid-gas critical point [16][17][18] are the most striking examples of achievements within the CV method.

Z.Gurskii
The partition sum of the three-dimensional Ising model was first rigorously obtained by I.Yukhnovskii, see [19,20].It allows one to get an equation for the temperature of the second order phase transition and perform numerical calculations of the Ising model thermodynamical properties in the critical region [15,[21][22][23][24].To put it in another way, the nonuniversal characteristics of the three-dimensional Ising model have been calculated in the critical region for the first time.This significant result stimulates investigations of phase transitions in different physical objects that can be described by the Ising-like models [25][26][27][28][29].
The CV method is also a powerful tool in constructing the microscopic theory for substitutional binary alloys [7,26,27].To show the main stages of this theory development within the CV method is the aim of the present paper.
Calculation of the alloy free energy as a function of temperature and alloy components' concentration is the key point in constructing alloy phase diagrams as well as in investigating thermodynamical properties.Three trends in the problem of alloy free energy determination can be distinguished.
1.In the "constructive method" [30] the free energy components (E, the internal energy, and S, the entropy) are calculated separately within quite different approaches.For example, one can get E by the pseudopotential method [7,31] and S using combinatorial arguments [30].
2. The analytical theory proceeds from the partition function Z: Here k B is the Boltzmann's constant and T is the temperature.The analytical theory is usually based on model assumptions [32].
3. Computer simulations have been exploited intensively for investigating different properties of alloys [33].
Each trend has its advantages and restrictions.For instance, the Ising model of the nearest neighbours (or of the nearest neighbours and next-to-nearest neighbours) is commonly used in the analytical theory [32].However such models cannot describe real alloys adequately because the pair interatomic effective interactions in alloys belong neither to the short-range potentials nor to the long-range ones [7,30].They need to take into account the interatomic interactions within 5 or 6 coordination spheres at least to reach reasonable results for atomic properties of metals and alloys [30,31].
In computer simulations one deals with systems of finite size and the problem always exists of extrapolating the results obtained to the thermodynamical limit when N → ∞, Ω → ∞, Ω/N = const, where N is the number of particles and Ω is the system volume.
What should one know to calculate the alloy thermodynamical properties at the microscopic level?In the first place we need information on interatomic potentials in the alloy.Besides, one should take into consideration local static displacements of the lattice sites and the atomic thermal vibrations [31].The CV method is efficient in taking into account all these factors [7,34,35].
The paper is organized as follows.The problem of deriving the interatomic potentials in metals and alloys is considered in section 2. Calculation of the binary alloy grand partition sum using the CV method within the rigid lattice approximation is a subject of section 3. Modification of the developed theory if the atomic static displacements (ASD) are taken into account is shown in section 4. The role of the atomic thermal vibrations (ATV) in binary alloy thermodynamics is discussed in section 5. Some possible ways regarding the further development of the given theory are considered in Conclusions.Acknowledgements complete the paper.

Many-body interatomic potentials in metals and alloys
Calculation of interatomic potentials is one of the central problems in modern condensed matter physics.It is caused by utilizing these potentials in numerous problems of physics and materials science in construction of alloy statistical theory on the microscopic level [7,31,32], computer simulations that use the molecular dynamics method [31], etc.In some cases an adequate theoretical description of metal and alloy properties could not be achieved within the pair interatomic interaction model [7].It means the many-body interactions (three-, four-particles) play an important role in these objects [36].
One can distinguish several trends in deriving interatomic potentials in metals.2.1.A rigorous determination of interatomic potentials in quantum mechanics is based on the pseudopotential method [37,38,7].Working within the framework of perturbation theory in pseudopotential, we can represent the metal total energy in the form of expansion in the effective multi-ion interactions [7,31,[36][37][38]: The first term E 0 depends on the atomic volume Ω 0 only and includes the energy of the interacting electron gas as well as all one particle interatomic contributions to the E tot .The interatomic potentials V 2 , V 3 etc.are implicitly volume dependent but explicitly structure independent and thus rigorously transferable at a given volume Ω 0 to the all bulk structures either ordered or disordered.To put it in another way the potentials themselves are independent of the absolute ion positions R i and depend only on relative separations 36,38].For example, the angular force triplet potential V 3 is a three-dimensional function: [36].The prime on each summation in (2.1) denotes exclusion of all the self-interaction terms where two indices are equal.The structure dependence of the total energy (2.1) appears through the summation in (2.1) over all N ion positions.One should stress the important feature of the equation (2.1) obtained within the pseudopotential theory.Expansion (2.1) is rapidly convergent in the sense that the four-ion quadruplet potential V 4 is smaller than V 2 , V 3 and higher potentials (V 5 , V 6 , . ..) appear to be negligible [36].

2.2.
A somewhat another approach to deriving the effective interatomic potentials in metals has been advanced in [39,40].It is based on utilizing the Hohenberg-Kohn theorem [41] and on treating the Kohn-Sham total energy functional [42] in terms of the many-body interatomic interactions.
Let us consider an electron-ion system within the pseudopotential concept.It means that pseudoions are placed in positions R and their interactions with the valence electrons are described by a pseudopotential w(r − R) which is an operator in a general case.
The Kohn-Sham total energy functional of such a system has the following form [39,40,42]: (all equations are written in the atomic units in the present paper).Here T [ρ] is a functional of the kinetic energy, the second and the third terms in (2.2) describe the energy of the electron subsystem in the external field of the pseudoions and the Hartree energy respectively.The last terms in (2.2) are E xc , the exchangecorrelation energy functional and E i−i , energy of the ion-ion direct interactions.The total energy E[ρ] (2.2) is the universal functional of an electron density ρ.
The present theory provides knowledge regarding the electron density ρ(r), that minimizes E[ρ(r)] (2.2), with high accuracy, see [43,44].In [39] the known quantity of ρ(r) has been presented in the form where ρ 0 = z/Ω 0 (Ω 0 = Ω/N) is the density of the uniform electron distribution, z is the ion valency, Ω 0 is the atomic volume and ρ i (r − R) is the electron density related to the pseudoion at the site R. Equation (2.4) indicates that the total valence electron density is represented as the linear superposition of contributions from the individual pseudoions, embedded in the uniform background of the density ρ 0 .The focal point of the linear superposition assumption (2.4) is an accurate real-space treatment of the metal total energy in terms of well-defined interatomic potentials [39,40].The local density approximation [42] and equation (2.4) make it possible to represent the T [ρ] and E xc [ρ] functionals in terms of appropriate contributions to the many-body interatomic interactions [45] and get the total energy functional (2.2) in the form of equation (2.1), see for details [39,40,45].The explicit expressions for the pair and triplet interatomic potentials are given in [39,45].One should stress that these results are obtained without utilizing the perturbation theory in pseudopotential.Unlike the perturbation theory the advanced method permits one to accurately calculate the so-called reducible contributions to the pair and triplet potentials arising from the n-particle interactions when (n − 2) or (n − 3) indices of ion co-ordinates coincide, respectively.The analytical expressions for these contributions are presented in [45].Thus, the present stage of theory allows one to determine the interatomic potentials in metal systems with high accuracy.

The alloy partition function within the collective variables method
Consider a substitutional binary alloy with one atom per elementary cell.Atoms of two kinds A and B are placed arbitrarily on N crystal lattice sites.Their configuration is given by the set of numbers σ R that equals +1 if the site R is occupied by an A-type of atom, and equals −1 otherwise.
One can separate the electronic and ionic degrees of freedom in non-transition metal alloys [7].The alloy Hamiltonian within the rigid lattice and pair interatomic interaction approximations has the following form after summing over electron states [7] H{σ Here V i,j (q) (i, j = A, B) is the Fourier transform of the effective interaction between ions of i-and j-types, V AB (q) = V BA (q).The explicit expressions for V ij (q) are given in [7], see also [46].
We proceed from the partition function to find the free energy (see equation (1.2)) The symbol Tr {σ R } denotes the sum over all possible values of the occupation numbers {σ R }, where β = (k B T ) −1 is the inverse temperature.This procedure is equivalent to summing over all possible configurations of two types of atoms on N lattice sites at the given alloy concentration c i = N i /N, where N i is the number of i-type atoms, i = A, B.
The calculation of Z within the canonical ensemble is a hard problem even in the case of finite small systems because the occupation numbers σ R are not independent.To avoid this difficulty, and to perform the Tr {σ R } operation on each site R, independent of the specific configuration, one has to work within the grandcanonical ensemble, introducing µ i , the chemical potentials of the alloy components, as Lagrangian multipliers.Then in view of (3.1) and the following conditions the grand partition function takes the form [7,46] The chemical potentials {µ i } are found from the condition The explicit expressions for V 0 , V 1 , V 2 (R, R ′ ) are given in [7].They have the following physical meaning: V 0 is the part of alloy energy which does not depend on the atomic configuration, V 1 (R) indicates the difference between alloy component atomic characteristics, and is the ordering potential; see [7] for details.In crystals with one atom per elementary cell V 1 (R) = const [7].
The next relation exists between the grand potential F = k B T ln Z and the free energy F in the thermodynamic limit [32] (3.7) The grand partition function Z, equation (3.4) is calculated by the collective variables (CV) method [6,7].The CV space {ρ k } is introduced in the following way [6] ρk = ρ k J(ρ k , σ r ) k∈BZ dρ k . (3.8) is the Fourier transform of the occupation numbers σ R .The values ρk depend on the atomic configuration in the alloy.
is the transition operator from the set {ρ k }, that is the Fourier components of the occupation numbers {σ R }, to the collective variables {ρ k }.The symbol k ∈ BZ denotes that the wave vector k takes N values in the first Brillouin zone (BZ), δ is the Dirac delta function.
Equation (3.4) for the grand partition function is rewritten as follows within the CV method where is the transition Jacobian to the collective variables {ρ k } space, V 2 (k) is the Fourier transform of the ordering potentials V 2 (R).One should not mix the variables ρ k with ρk (3.9) that are the configuration-dependent quantities.
The general ideas for calculating J(ρ) are given in [6] and [7].
The transition Jacobian to the CV space can be presented in the following way with where a n are the coefficients of J(ρ) and As follows from equations (3.14) and (3.15), the cumulants M n and the coefficients of the transition Jacobian a n appear to be complicated functions of the potentials Finally, in view of (3.11) and (3.13), the problem of calculating the alloys' grand partition function is reduced to the following integral form [46] A natural question arises.Which distribution function should one take in J(ρ) (3.13) to achieve an adequate description of alloy thermodynamic properties using the CV method?To put it another way, how many coefficients a n (n = 1, 2, 3, 4, . ..) does one have to take into account in the series (3.13)?
It was shown in [6,7] that the Gaussian distribution (a 1 , a 2 = 0, a 3 = a 4 = . . .= 0) is appropriate in a wide temperature range both above the temperature of the phase transition T c , and below T c (T > T c and T < T c ).They had to use the measure ρ 4 (a 3 , a 4 = 0) in equation (3.13) just in a narrow temperature interval τ = |T − T c |/T c ≈ 0.01 to get the correct values of the critical exponents [6,7].
Let us analyse the alloy thermodynamic properties at T > T c .All calculations can be performed analytically within the Gaussian approximation of the CV method T > T c .One has to restrict oneself to the cumulant M 2 in (3.14), that is put M 3 = M 4 = . . .= 0.Then, Further on, the index G will be omitted in expressions (3.17) and (3.18) to simplify notations, because the next analysis concerns the Gaussian approximation only.
Calculation of the alloy's grand partition function (3.16), in view of (3.17) and (3.18), does not encounter any difficulties.Then, the grand potential F per atom is The equation, see (3.6) determines the difference of alloy component chemical potentials at the given alloy concentration.In view of (3.19), equation (3.20) is written as follows where x ≡ βV 1 (µ), and see [46].
One has to solve a system of equations (3.19)-(3.21) to find the alloy free energy as a function of temperature and alloy concentration.The binary correlation function σ R σ R ′ is an important characteristic of disordered alloys.It is related to the short-range order parameter [7].The calculation of σ R σ R ′ Fourier components is fairly trivial within the Gaussian approximation of the CV method [7,27,46] One can obtain the equations for T c , the temperature of the order-disorder phase transition proceeding from the fact that ρk * ρ−k * diverges at T = T c .Here {k * } are the points of the Brillouin zone where the ordering potential V 2 (k) has the absolute minimum [7,25].The equation determining T c follows from (3.23) At first glance, equation (3.24) is very much similar to the one obtained by the method of static concentration waves [47] 1 The cumulant M 2 in (3.24) is a complex function of temperature, alloy concentration, and the potential V 1 .Therefore (3.24) is a complicated non-linear equation with respect to the temperature.It is reduced to the form of equation (3.25) in the hightemperature limit only [7,46] However, inequalities (3.26) become true just at T > T L for many alloys, where T L is the liquidus temperature.Thus, conditions (3.26) ascertain the limits of equation (3.25) applicability and allow one to understand why the static concentration waves method does not work in describing the order-disorder phase diagrams of polyvalent metal alloys [46].The theory given above was exemplified in [46] by investigations of Ca-Ba and K-Cs systems.Both systems are characterized by a great mutual solubility of components.The system of equations (3.19)-(3.21)are solved to find the free energy as a function of temperature and concentration c.The equilibrium atomic volume is found from the condition The heat of alloy formation calculated at different values of temperature and concentration for K-Cs and Ca-Ba systems in [46] indicates the tendencies which perfectly agree with the experimental data [48].In (3.28) F i (i = A, B) denotes free energy of the i-type of pure metal.
One can construct the order-disorder phase diagram plotting T c against alloy concentration.The T c temperature is determined from equation (3.24).As the cumulant M 2 is a complicated function of the chemical potentials and temperature, equations (3.19)-(3.21)and (3.24) and (3.27) are solved simultaneously.The reasonable agreement between theoretical and experimental data for T c = f (c i ) is achieved for K-Cs and Ca-Ba systems in [46] (always within 30 %).
All numerical results presented in [46] demonstrate the efficiency of the approach which is based on the CV method.

Role of atomic static displacements in alloy thermodynamics
It is well known that the formation of metal solid solutions is accompanied by arising of local lattice distortions.The latter are characterised by atomic static displacements (ASD) with respect to the ideal mean lattice sites.The ASD have a drastic effect on the X-ray (neutron) diffuse scattering [49].They determine a lattice parameter dependence on alloy concentration.That is why including the ASD into consideration is of great importance in constructing a consistent microscopic theory of binary alloys, see for more details [35] and [50].
Let us take into account the fact that the local ASD are present in an alloy.Then, the coordinates of the lattice sites are the following ones: where δR are the ASD with respect to the sites R 0 of the ideal mean lattice.Assume that δR does not depend on the type of an atom and perform the Fourier transformation of δR: The ASD (δR; as well as δR k ) are random quantities in a disordered alloy.Let us separate A k , a configurationally independent part of δR k , by means of the following relationship [35] Here, compare with (3.9),  [35,49].
Let us expand the factor exp(iqR) in (3.1) in power series of the static displacements δR, restricting ourselves to the square of δR.The alloy Hamiltonian H(σ) (3.1) with allowance for (4.1) to (4.4) in the harmonic approximation takes the form [35] H where is the Hamiltonian of an ideal mean lattice without displacements, see section 3.
The addends H 1 (k, A k , ρk ) and H 2 (k, A k , ρk ) are linear and quadratic in A k amplitudes, respectively.The explicit equations for them are given in [35].
We proceed from the grand partition sum to find the free energy, see (3.4)One can rewrite equation (3.4) with a view of (4.5) and using the rigid ideal lattice of an alloy as a reference system, as follows: (4.7) Details are given in [35].Here are the addends of the alloy Hamiltonian (4.6) renormalized by the ASD and The following notations are accepted in (4.8)-(4.11):G are the reciprocal lattice vectors, and Φ (0) -the force constant matrix of the reference system.The correlated average crystal (CAC) in the rigid lattice approximation is used as a reference system.One can get familiarized with the CAC term value in [7,34].The expression for Φ (0) is given in [34], also see appendix 2 in [35].
The grand partition sum (4.7) is calculated by the CV method.Equation (4.7) is rewritten in the following way within the CV method [50], 12) compare (4.12) with (3.11).
Calculation of the grand partition sum (4.12) can be performed analytically in the Gaussian approximation.Details of the consideration are omitted because they are similar to those given in section 3. Then the grand potential per one atom equals Here M n are cumulants, see (3.15).Equation (3.20) determines the difference of alloy components chemical potentials at the given alloy concentration.
One has to perform the Legandre transformation and solve equation (3.20) to find the alloy free energy F (T, C) as a function of temperature and component concentration.Then, see [50] for details where and x = β Ṽ1 (µ), see (4.9), is the solution of the system of equations and (3.20).Solution of equation (4.16) is given in [35].We present the final result omitting details Here ε kλ and ω 2 kλ are eigenvectors and eigenvalues of the force constant matrix Φ (0) , respectively, λ = 1, 2, 3 -the polarization index and  is the average ion mass, see [7] for details.Analyse result (4.17).One can conclude from (4.11) and (4.17) that the ASD amplitudes A k are small if the pair interatomic potentials V AA and V BB Fourier components are similar: V AA (q) ≈ V BB (q).Really, P k ≡ 0 at V AA (q) = V BB (q) and then A k = 0.This conclusion allows one to clear up the nature of the wellknown phenomenological Hume-Rothery rules [51] on the microscopic level.Using equations (4.9), (4.15) and condition (4.16) one can prove that It means that the potential Ṽ1 does not depend explicitly on the ASD amplitudes A k .This result greatly simplifies the calculation of the alloy free energy (4.The Fourier components of the binary correlation function are important alloy characteristics.They are needed to calculate the X-ray (neutron) diffuse scattering intensity [49,52].Besides, they are related to the SRO parameter (4.22)Here α R is the value of the SRO parameter on the R-coordination sphere, ρ k ρ −kthe Fourier components of the binary correlation function.Calculation of ρ k ρ −k does not face any difficulties within the Gaussian approximation of the CV method [7], see also section 3 (4.24)It is seen from (4.24) that Ṽ2 < V 2 in the whole first Brillouin zone except for the points of high symmetry where vector P k = 0 [35].One can notice from (4.23), (4.24) and (4.17) that the binary correlation function Fourier components directly depend on the ordering potential, temperature and the ASD.Besides, they are complicated functions of potential V 1 and the alloy concentration via cumulant M 2 and the equilibrium atomic volume.In the present section the theoretical results are illustrated by numerical calculations performed for the alloys of K-Cs and Ca-Ba systems.Solid solutions of the body centred cubic (bcc) structure exist in wide ranges of temperature and alloy concentration in both systems [48].The renormalized potential Ṽ2 (k, A k ) (4.24) for K-Cs and Ca-Ba alloys are shown by dashed lines in figures 1 and 2, respectively.The bare ordering potentials Ṽ2 (k) are depicted by full lines.Details of the calculations are omitted because they are described in [7,35,46] The given results can be summarized in the following statements.
1.The ASD have a drastic effect on the binary correlation function Fourier components ρ k ρ −k behaviour in the first Brillouin zone.They smooth the dispersion of the ρ k ρ −k = f (k)-function.
2. Tendency to ordering becomes more pronounced owing to the ASD in the alloys of K-Cs and Ca-Ba systems.
3. Dependence of the SRO parameter on temperature obtained numerically in [50] agrees with the conclusions drawn from the treatment of X-ray diffuse scattering experiments [52].

Lattice-vibrations' contribution to the free energy of disordered binary alloys
The most calculations of the alloy free energy have been performed within the rigid lattice approximation, see, for example [7,26,27,[30][31][32].The first attempts to incorporate the vibrational contributions to the alloy free energy were based on empirical methods or phenomenological models, see review of literature in [54].Ab initio approach to take the atomic thermal vibrations (ATV) into account was advanced in [54].The main ideas of [54] are presented shortly in this section.
There are two kinds of freedom degrees in disordered alloys: configurational and vibrational ones.They should be incorporated into computational scheme on the same microscopic level to construct a consistent theory [31,54].The main aim of the present section is to calculate the substitutional binary alloy free energy if the configurational and vibrational degrees of freedom and interrelations between them are taken into account explicitly.Some features of it are worthy of attention.
(i) We start from the problem for the grand partition sum of substitutional binary alloy.One can perform summing over vibrational and configurational states separately because of the two quite different time scales characterizing them [55].
(ii) The correlated average crystal (CAC) model [7,34] is used here to sum over the vibrational degrees of freedom.This microscopic model is more adequate than the Einstein one or phenomenological Debye-Gruneisen approach.The shortrange order effect on the phonon density of states, the Debye temperature and other quantities can be naturally investigated within this model [7,34].Analysis of the relationships between CAC approach and the known virtual crystal and coherent potential approximations is given in [34].
(iii) Conditions when interrelations between the configurational and vibrational degrees of freedom can be ignored have been formulated for the first time in the present section, see also [54].
Substitutional binary alloy Hamiltonian with the ATV included into consideration has the following form in the pair interatomic interaction approximation [7,54]. and Here T and U are the kinetic and potential energies, respectively which depend on configurations of two types (A, B) atoms on N lattice sites R via sets {σ R }.Each configuration is characterized by a certain set of {σ R }.Such notations are used in ( is the mass of i-type atom, δR atomic displacement from the ideal lattice site R 0 caused by the thermal vibrations, symbol δ Ṙ means the derivative of δR with respect to time, and V ij (q) (i, j = A, B) is the Fourier transform of the effective pair interaction between ions of i-and j-types.We assume in (5.2) that δR do not depend on the atom type.
Introduce the Fourier components of σ R and δR ) It should be emphasized that δR k (5.5) implicitly depend on an alloy configuration.
To put it another way, they are the configurationally dependent quantities.Equation (5.1) with allowance for (5.2) to (5.5) takes the form in the harmonic approximation [54] where U 0 (ρ) is the potential energy of the rigid lattice, T (k) + U(k) are Fourier components of the atomic thermal vibrations Hamiltonian H th .
Expression for U 0 (ρ) will be given below, equations for f (i) (k) functions (i = 0, 1, 2) are presented in [7].The terms T (1) (k), U (1) (k) and U (2) (k) proportional to ρk and ρk ρk ′ are concerned with the long-range order and the short-range order effects accordingly, while T (0) (k) and U (0) (k) depend on an alloy configuration implicitly via δR k δR k ′ .Interrelations between δR k with different indices k are present in the T (1) (k), U (1) (k) and U (2) (k) terms because of the absence of the translational invariance in a disordered alloy, see (5.7), (5.8).The same takes place for ρk .It means that the wave vector k is the adequate ("good") quantum number just in the case of an ideal crystal.Nevertheless, the Fourier transformation of the Hamiltonian is also advantageous in the case of a disordered alloy, so long as it allows to treat interatomic interactions via V ij (q) formally exactly, see (5.3) [7,47].As is well known, the pair interatomic potentials in alloys neither belong to the short-range potentials nor to the long-range ones [7,30,31].They need to count interatomic interactions within 5 or 6 coordination spheres at least to reach reasonable results for metal and alloy atomic properties [7,30,31].This problem is naturally solved by utilizing the Fourier transformed alloy Hamiltonian and the collective variables (CV) or static concentration waves methods [7,47].Let us search the unknown amplitudes of the atomic thermal vibrations δR k in the form of expansion in the complete set of some functions.The polarization vectors e kλ of the CAC vibrations generate such a complete set of the basis functions [21].They are the eigenvectors of the following eigenvalues problem [34] where Φ (0) (k) is the CAC dynamic matrix, M = i=A,B C i M i the average ion mass, ωkλ the frequencies of vibrations and λ = 1, 2, 3 is polarization.Equation for Φ (0) (k) is presented in [7]. Then with a kλ the expansion coefficients.
A correct transition to the quantum mechanical formulation of the problem under consideration can be performed based on the complete set {e kλ }.The operator is put in correspondence with the amplitude δR k (5.10).Here B kλ are unknown dimensionless expansion coefficients, b + −kλ and b kλ operators of creation and anni-hilation of the vibrational excitation ωkλ , accordingly.Obviously B kλ ≡ 1 for the CAC case.
Momentum P = M δ Ṙ is the conjugate quantity to δR.The corresponding operator P is written as (5.12) for the coefficients B kλ follows from the commutators [δ R, P ] = ih and The problem of the partition sum calculation within the grand canonical ensemble reads as where symbols Tr {σ R } and Tr {P h} denote summing diagonal matrix elements of the statistical operator exp(−βH) over all possible configurational and vibrational alloy states, respectively.The configurational and vibrational degrees of freedom in (5.14) are characterized by two quite different time scales.The typical times for the lattice vibrations are of the order of 10 −13 sec., while substitutional interchanges occur in times that are several orders of magnitude longer [55].That is why one can assume that vibrational states are ergodic over the time scale of substitutional configurations.Therefore the Tr {P h} operation should be performed at first in equation (5.14) at a definite alloy configuration.As a result, the effective Hamiltonian renormalized by the atomic thermal vibrations will be received.
The complete system of |n kλ -functions which are the eigenfunctions of b + kλ , b kλ operators is used for performing the Tr {P h} procedure: where (5.17) Here n kλ is the number of vibrational excitations characterized by the quantum numbers k and λ.One should remark that |n kλ are the eigenfunctions of the CAC (reference system) Hamiltonian and they are not the eigenfunctions of the H th (5.7), (5.8) written in the secondary quantization presentation.Nevertheless calculation of the H th diagonal matrix elements ( n kλ | Ĥth |n kλ ) does not face any difficulties with utilizing (5.16), (5.17) and orthonormality of the basis functions [54].
The diagonal matrix element of Ĥth has the following form [54] Appearance of the second term in the equation (5.18) is caused by the fact that |n kλ are not the eigenfunctions of H th .Consider such temperatures when βH th ≪ 1 and they may restrict themselves to the first order of the thermodynamical perturbation theory (TPT).Then, see equation (5.15) Tr (5.20) Equation for the grand partition sum takes the form after performing the Tr {n kλ } procedure [54].
and n kλ | Ĥth |n kλ given by equations (5.18) and (5.19).The terms entering equation (5.22) have the same physical meaning as for the rigid lattice case, see section 3. Equation (5.21) can be rewritten as follows [54] where Replace for a moment n kλ in (5.19) by its mean value nkλ to carry out such an analysis.First of all the absolute value of ϕ kλ function (5.19) is much smaller than |V 2 (k)| [34].Therefore, renormalization of the ordering potential is a rather small effect and one should regard the rigid lattice as a reasonable zero approximation in the alloy thermodynamics [7,[30][31][32].Analysis of the f (2) (k) function entering equation (5.19) for ϕ(k, λ) shows that renormalization of the ordering potential V 2 (k) by the atomic thermal vibrations is more pronounced in alloys where there is a significant difference in the interatomic pair potentials [7].Thermal vibrations should play a more important role in the thermodynamics of such alloys, especially in the phase diagram calculations.Numerical results of [56,57] perfectly confirm this conclusion obtained within the analytical theory.
We shall calculate the grand partition sum (5.23) using the CV method [6,7].This procedure does not differ from the cases described in section 3 and 4. That is why we present the final result for the grand potential per one atom obtained within the Gaussian approximation of the CV method.See for details [54]. where and The terms V 0 , V 2 (0) in (5.27) contribute to the alloy internal energy while the terms proportional β −1 contribute to the alloy entropy.The last summand in (5.27) is caused by the atomic thermal vibrations.At the first glance it has the similar structure as for the case of an ideal metal [58].Let us consider the L(k, λ) function (5.29) to stress the features of the advancing approach.
The second term in the parentheses in the right-hand side of equation (5.29) is caused by interrelations between configurational and vibrational degrees of freedom in an alloy.If one neglects this term equation (5.27) becomes F (T, µ) = F (0) (T, µ) + F ph (5.30) with F (0) the grand potential of the rigid lattice and the vibrational free energy.Therefore inequality formulates the condition when configurational and vibrational contributions to the alloy free energy may be treated separately.The right-hand side of (3.7) becomes as follows, with a view of (5.27) and (5.24) The last term in (5.34) describes the contribution of the atomic thermal vibrations to the chemical potentials difference (µ A − µ B ). Solving equation (5.34) at specified (C A − C B ) value we find x = βV 1 (µ), that is (µ A − µ B ), see (5.24).Let us analyse equations (5.32) and (5.34) more in detail to elucidate relationships between the well known alloy theories [30][31][32]47] and the approach developed.
Present some exact formulae following from (5.24) and (5.33) which will be used further on (5.36) dependence of the disordered alloy formation energy on the two effects: the partial order and the size mismatch.Their numerical results are confirmed completely by the analysis of equation (5.40).By the way, calculations in [59] were carried out at temperatures which meet the inequality (5.37).Consider now the vibrational free energy in equations (5.32) and (5.40).The argument of the ln sinh[βL(kλ)/2] function must satisfy the following inequality βL(kλ) ≪ 1 (5.42) as far as the expression for F (T, C) has been obtained within the first order of the thermodynamical perturbation theory, see (5.20).It means that including the atomic thermal vibrations into consideration stabilizes a disordered alloy in comparison with the rigid lattice model: ln sin h[βL(kλ)/2] < 0. This deduction is in complete agreement with the results of [56] and also [30].The Fourier components of the binary correlation function are defined by the equation with allowance for (4.23) and (5.32) [54] ρk ∂L(k, λ) ∂βV 2 (k) , k ∈ BZ, k = 0.
(5.43)The second term in (5.43) is due to interrelations between configurational and vibrational degrees of freedom.One should expect that it is a rather small effect, see the explicit expression for the L(k, λ) function (5.29).
Equation for T c the temperature of the order-disorder transition proceeds from the fact that ρk * ρ −k * diverges at T = T c if transformation takes place as the phase transition of the second order [6,7].
The equation determining T c 1 + β c V 2 (k * )D 2 = 0 (5.44) follows from (5.43).Here {k * } are the points of the Brillouin zone where the ordering potential V 2 (k) has the absolute minimum [7,25].The cumulant D 2 in (5.44) is a complex function of temperature, alloy concentration, potential V 1 and hω kλ .Therefore (5.44) is a complicated non-linear equation with respect to the temperature.It is reduced to the form (3.25) in the high temperature limit only: D 2 → 4C A C B if βV 2 (k) ≪ 1.By the way, the construction of the order-disorder phase diagram for Ca-Ba system in [46] indicated a satisfactory agreement with the experimental data.But theoretical values of T c , calculated according to equation (3.24) within the rigid lattice model, exceed the experimental ones [46].
It is of interest to answer a question: how will the inclusion of the atomic thermal vibrations into computational scheme vary results of calculations?We have to investigate the effect of the atomic thermal vibrations on a D 2 cumulant value.Consider equation (5.34) determining D 1 and D 2 , see (5.33) and analyse the L(kλ) function (5.29) for this purpose.Taking into account the explicit expression for f (2) (k) function and the fact that V 2 (k) < 0 one can derive the following inequalities f (2) (0) − f (2) (k) = < 0, k * = 0, > 0, k * = 0 . (5.45) The order-disorder phase transition occurs in alloys if the absolute minimum of the V 2 (k) potential is in the star {k * } with |k * | = 0 [7,30,47].The case when the ordering potential V 2 (k) has the absolute minimum in the centre of the Brillouin zone is related to decomposition of the homogeneous disordered alloy on two separated phases with the decrease of temperature [7,47].This case will not be examined in the present paper.Analysing the effect of the thermal vibrations on solutions of equation (5.34) with respect to D 1 , one can state with a glance to (5.45) and (5.33) that D 2 (r.l.) > D 2 (vib.). (5.46) Here notations D 2 (vib.) and D 2 (r.l.) mean the cumulant D 2 values obtained with and without thermal vibrations, respectively.An important conclusion follows from (5.44) and (5.46): the inclusion of the lattice vibrations reduces the temperature of the second order phase transition in binary alloys.
One can present several examples that show a significant role of the lattice vibrations in determining T c and confirm the trend stated above, see [31,56,60].A noticeable improvement of the results of the Bragg-Williams theory was achieved in [60] by the inclusion of the thermal vibrations even within a very simple Einstein model.The lowering of calculated T c values was observed for MgCd alloys [60].
Finally, the unknown correlation functions B * kλ B kλ present in the expression for the L(kλ) function (5.29), see also (5.40), are to be found to complete the discussion.A system of equations ∂F (T, C) ∂ B * kλ B kλ = 0, k ∈ BZ (5.47) determines the B * kλ B kλ quantities.Actually, B * kλ B kλ are dimensionless weight factors of the Fourier components of the Green function "displacement-displacement" averaged over configurational states [61].Equations (5.47) are solved numerically.

Conclusions
The results presented above can be summarized as follows.The CV method is successfully applied for investigating the thermodynamics of substitutional binary alloys.A new approach to calculating binary alloy thermodynamical properties, when the atomic static displacements (ASD) and the atomic thermal vibrations (ATV) are taken into consideration, is advanced.It gives good prospects for an adequate description of real metal systems on a microscopic level.
The following problems need to be developed.
1. Role of many-body interatomic interactions (three-and four-particles) in the alloy thermodynamics.
2. Inclusion of the ASD and the ATV into consideration simultaneously.
3. Analysis of the ASD and the ATV effect on alloy phase diagram within the "ρ 4 " approximation of the CV method.

Figure 1 .
Figure 1.Behaviour of the ordering potential Fourier transform V 2 (k) in the [111] direction in alloys of K-Cs system at T = 250 K. Dashed and full curves show results obtained, respectively, with and without the ASD taken into account.Curves 1 refer to alloy K 0.7 Cs 0.3 while the curves 2 correspond to alloy K 0.1 Cs 0.9 .

Figure 2 .
Figure 2. Behaviour of the ordering potential Fourier transform V 2 (k) in the [111] direction in alloys of Ca-Ba system at T = 750K.Notations are the same as in figure 1. Curves 1 refer to alloy Ca 0.5 Ba 0.5 and the curves 2 correspond to alloy Ca 0.2 Ba 0.8 .
14).Let us more carefully analyse equation (4.14) for the alloy free energy.The third term in (4.14) proportional to β −1 is entropy (S), while the rest of the terms define the alloy internal energy (E).One can get the next formulae for E and S considering equations (4.14) and (3.20) in the high temperature limit: βV 2 (k) ≪ 1.

. 21 )
Equation (4.20) determines the energy of an average crystal: all the lattice sites are occupied by mean ions which interact via the mean potential.v= v A C A + v B C Bwith v i -the potential of an i-type ion.Equation (4.21) defines the configurational entropy of an ideal binary solution.Thus, the high temperature limit of the CV method Gaussian approximation is equivalent to the well-known W.Bragg -E.Williams theory which ignores the pair atomic correlations.By the way, the difference ∆F = F (T, C) − E id + T S id with F (T, C) (4.14) indicates the contribution of the short-ranger order (SRO) effects to the alloy free energy.

Figure 3 .
Figure 3. Temperature effects on the binary correlation function Fourier components in the K 0.7 Cs 0.3 alloy.Dashed and full curves show results obtained, respectively, with and without the ASD taken into account.Curves 1 refer to T = 300 K and curves 2 correspond to T = 200 K.

Figure 4 .
Figure 4.The same for the K 0.1 Cs 0.9 alloy.

Figure 5 .
Figure 5.The same for the Ca 0.5 Ba 0.5 alloy.Curves 1 refer to 850 K while curves 2 correspond to T = 650 K.

Figure 6 .
Figure 6.Dependence of the binary correlation function on atomic concentration in the Ca-Ba system alloys.Direction [111] of the Brillouin zone.Notations are the same as in figure 3. Curves 1 refer to the Ca 0.2 Ba 0.8 alloy and curves 2 correspond to the Ca 0.8 Ba 0.2 alloy.
. The potential Ṽ2 (k, A k ) has the absolute minimum in the [111] direction in the alloys of the systems investigated.The ASD smooth the dispersion of the ordering potential V 2 (k) in the first Brillouin zone, especially in the [100] direction.It is seen from figures 1 and 2 that an additional minimum appears owing to the ASD in the [111] direction in the alloys of K-Cs and Ca-Ba systems.Dependence of Ṽ2 (k, A k ) on the atomic concentration is more pronounced in K-Cs alloys than in Ca-Ba ones.Behaviour of the binary correlation function Fourier components ρ k ρ −k in some principal symmetry directions has been investigated according to equation (4.23) for K-Cs and Ca-Ba alloys.Calculations have been performed with the ASD taken into consideration (dashed curves), and without them: A k ≡ 0 for k ∈ BZ (full curves), see figures 3-6.Drastic effect of the ASD on the ρ k ρ −k = f (k) behaviour is observed, especially for Ca-Ba alloys, see figures 3-6.The ASD encourage the gaining of the ρ k ρ −k values in the whole first Brillouin zone.They smooth dispersion of ρ k ρ −k in the alloys studied.Thus, the alloys become more similar to the ideal solutions owing to the ASD, especially at high temperatures.The ρ k ρ −k = f (k) functions strongly depend upon the atomic concentration and temperature in alloys of the both systems investigated, see figures 3-6.Dispersion of ρ k ρ −k becomes more visible at the decrease of temperature, see figures 3-5.The ρ k ρ −k = f (k) functions display the most interesting behaviour in the [111] direction.The largest effect of the ASD on ρ k ρ −k is observed in the [110] direction of the first Brillouin zone, see figures 3-6.
Component B is regarded as an "impurity".One should emphasize that the approximation (4.3) works very well in the whole region of C B values (0 < C B < 1) in such alloys where the dependence of the mean lattice parameter on C B is close to the linear one .4) is the k-th Fourier component of the occupation numbers, δ 0,k , the Kronecker symbol.Equations (4.3) and (4.4) indicate that the Fourier components of the local lattice distortions are caused by fluctuations of the impurity concentration waves with respect to the average value C B = N B /N.
.26) In general, B * kλ B kλ is a functional of the alloy short-range order.We replace B * kλ B kλ by their average values B * kλ B kλ in (5.18), (5.19) to simplify the problem (5.21).Equation that determines B * kλ B kλ value will be received below.Equation (5.26) defines the ordering potential renormalized by the atomic thermal vibrations.Let us analyse the obtained result (5.26) more in detail.