A new approach to the calculation of the thermodynamic potential of inhomogeneous electron gas

The main task of the equilibrium statistical theory of inhomogeneous systems is to calculate the thermodynamic functions and statistical distribution functions [1]. Application of the functional integration method for these calculations [2,3] enables us to obtain cluster expansion for these functions [1], based on the screened potential of interaction g(r1, r2). The equation for g(r1, r2) is an integral equation of convolution type whose analytical solution is a non-trivial problem.


Introduction
The main task of the equilibrium statistical theory of inhomogeneous systems is to calculate the thermodynamic functions and statistical distribution functions [1].Application of the functional integration method for these calculations [2,3] enables us to obtain cluster expansion for these functions [1], based on the screened potential of interaction g(r 1 , r 2 ).The equation for g(r 1 , r 2 ) is an integral equation of convolution type whose analytical solution is a non-trivial problem.This paper deals with the problem of calculating the thermodynamic potential Ω of an spacebounded electron gas using the method of functional integration.A new approach is proposed to calculate the thermodynamic potential, which consists in reducing the non-Gaussian functional integral to its Gaussian form with a renormalized "density-density" correlator.An integral equation has been obtained for the "density-density" correlator in an approximation which follows the approximation of random phases (RPA); in a low temperature case, its analytical solution has been obtained in "constant density" approximation.(

The Hamiltonian of the system
The Hamiltonian of such a system is as follows: where is operator of kinetic energy of electrons (m is mass of electron), is potential energy of electron-electron interaction, is potential energy of positive charge, is energy of interaction of electrons with positive charge, ϕ(r j ) is the potential created by positive charge.Let us assume that this charge is distributed as follows: where 0 = 2N SL e, θ(−Z) is Heavyside function, ∆ (R) describes positive charge distribution in the transitional layer from 0 inside the metal to 0 outside the metal.It is clear that this deviation is concentrated near the plane Z = 0 and, in fact, it forms a single-particle potential for electrons, which we call a surface potential and model it as a potential barrier.Taking this into account, we obtain: where . Therefore, we suggest that the surface potential V (r) (which is related to ∆ (R)) should be modelled as some model potential barrier, thus making it energetically unfavourable for the electron to leave the metal.
In order to calculate a grand statistical sum of this system, it is convenient to present the Hamiltonian (2) as a secondary quantization representation.For this purpose, we introduce single-particle wave functions Ψ a (r) and corresponding energies E a of the electron in the field of the surface V (r): which we use to construct the representation of the secondary quantization.Due to the symmetry of our problem it can be assumed that V (r) = V (z), where z is the electron coordinate normal to the plane XOY , r = (r || , z), r || is two-dimensional radius-vector in a plane parallel to XOY .In this case, the energy of the electron may be represented as follows: where p is the electron impulse in a plane parallel to XOY , α is a quantum number that depends on a specific form of the surface.The wave functions are as follows: We model the surface potential as a potential step of the height W since this model is relatively simple and adequately reflects the true surface potential.That is, we assume that The eigenfunctions and eigenvalues of this potential are given by the following expressions: where the domain of the values of the coordinate normal to the plane XOY is halfclosed interval [−L/2, +∞); α is quantum number, it assumes values that can be determined from the following equation: which follows from the condition that ϕ α (−L/2) = 0; quantity s is related to barrier height W : and it is as follows: In all the real problems, the potential barrier height is greater than the Fermi level, that is W > µ [4].

Functional representation of a grand statistical sum
The grand statistical sum that determines the thermodynamic potential of the system: Ω = − ln Ξ/β and other thermodynamic functions, in the interaction picture becomes: where where T is the symbol of chronological ordering of "times" β = 1/θ, θ is thermodynamic temperature.
For further calculations it is convenient to go over to spectral representation: where ν = 2π β n (n = 0, ±1, ±2, . ..) are Bose frequencies.Then (23) becomes: In order to simplify approximation S(β) according to (22) we go over to functional representation for S(β) [2,3], using Stratonovich-Hubbard identity [5] exp − 1 2 where i is imaginary unit, A = ||A n,m || is symmetric matrix, then for (26) we obtain: where (dω) is phase space element Note that due to the fact that operator variables ρ k (q|ν) are in (28) under sign T -ordering, it is impossible to integrate by β in (25).
For convenience we introduce the following notation: x = (q, ν), then averaging S(β) according to (22), we obtain: where ..+xn,0 are the socalled irreducible mean values (cumulants); the stroke near the sum also means that q = 0, and consequently M k 1 (x 1 ) = 0. Calculation of M k 1 ,...,kn (x 1 , . . ., x n ) is made in [6].For convenience we introduce the following notations: , then, the integral from expression (32) becomes: Calculation of integral (32) is a complicated problem due to the exponential index having terms with n 3. Their elimination results in a Gaussian approximation (or the so-called approximation of random phases).There are various approaches that aim to go beyond this approximation, that is to take into account the terms from n 3 [1,7,8].As a rule, this is done by means of series expansion of the non-Gaussian part of the integral (32) with subsequent averaging in Gaussian distribution and partial summing up of the terms which give the most important contribution.In order to avoid these complications we approximate the non-Gaussian integral (32) with the Gaussian one, introducing the unknown function D 1,2 : which we shall seek from the condition that the mean value from ω 1 ω 2 , calculated with non-Gaussian distribution J 1 (ω) is equal to the mean value calculated with Gaussian distribution J 2 (ω): where we introduce the following notations: . . .= (dω)J 1 (ω) . . .
The mean value of ω 1 ω 2 G is quite simple to calculate since the averaging is done with Gaussian distribution: where G 1,2 is effective potential (with an accuracy of up to factor β/S), Calculation of the mean value of ω 1 ω 2 is much more complicated, and for this we use the approach developed in [9,10].The mean value of ω 1 ω 2 is directly related to the two-particle correlator M 1,2 = i 2 T ρ 1 ρ 2 , where the averaging is done with the Hamiltonian of the whole system as follows: Function ω −1 ω −2 is unknown, and in order to find it in [9,10] a system of equations was proposed which are treated in the next chapter.
In matrix form, (38) can be written as: and from (37) we find that From the condition of equality of mean values of ωω G and ωω (expression (35)) we obtain: or where I is the unit matrix.Thus, an important result has been obtained: in order to go beyond RPA one may be confined to the Gaussian form of the functional integral, substituting the two-particle "density-density" correlator (in which the averaging was made with the Hamiltonian of an ideal system) with the renormalized "densitydensity" correlator in which averaging was made with the Hamiltonian of the whole system.
Then, we obtain the following expression for the thermodynamic potential: where

Calculation of effective potential
As shown in [11,7,12] for a correct description of exchange-correlation effects it is sufficient to be confined to four-particle correlations, i.e. in expression (33) to put M 0 1,...,n = 0, n 5. Therefore, in our work we shall confine ourselves to the same approximations in order to calculate the mean value of ω 1 ω 2 .
In this case the mean value of ω 1 ω 2 is in a complicated way expressed through the mean values of ω 1 ω 2 ω 3 and ω 1 ω 2 ω 3 ω 4 , for which in [9] a chain of equations was proposed, there being a certain arbitrarity in the choice of the concrete forms of those equations, which is to a certain extent similar to the ambiguity of splitting of Green's function.In [10] a regular procedure is proposed for break and closure of this chain of equations.Similar to [10] we choose the following equations: where (see in (33)).
We solve the system of equations (44) seeking solutions for ω 1 ω 2 , ω 1 ω 2 ω 3 and ω 1 ω 2 ω 3 ω 4 as: where notation . . . 0means solution in the approximation of zero sums in x = (q, ν), all the sums in k being taken into account (we cannot confine ourselves to a finite number of sums in k since index k is responsible for the inhomogeneity of the system, and the series will not coincide in the domain of rapid change of electron density).Unlike this, vector q is responsible for the homogeneity of the system (q is a vector parallel to the interface) and expansions will be made in plasma parameter λ = τ /τ TF , which is equal to the ratio of the mean distance between particles to the Thomas-Fermi screening radius (mean distance τ = (3/4πn) 1/3 , Thomas-Fermi screening radius τ TF = 1/κ TF , n is concentration of electrons).
For further discussion it is convenient to introduce the following notations: In the approximation of zero sums in x equation ( 44) for l = 2 in those notations becomes: where P 1,2 (1, 2) is the operator of cyclic permutation.It acts on arbitrary function f 1,2 (1, 2) as follows: ).We solve equation ( 50) by iterations, assuming that where notation . . . 0,i means the following: approximation of zero sums in x and i sum in k.Substituting (51) into (50) and collecting the whole infinite series we obtain: where function R 1,2 (1, −1) satisfies the following integral equation: The expression for ω −1 (−1)ω −2 (−2) 1 (59) includes the mean value from three field variables ω in the approximation of zero sums in x.In order to find this mean value we consider equation ( 44) for l = 3 in this approximation: where the following notation is introduced: Solving equation ( 60) by means of iterations, we find an expression for in the approximation of zero sums in x: Similarly, one can find ω −1 (−1)ω −2 (−2) 2 ; for this purpose one must consider equation (44) for l = 2 in the approximation of two sums in x, find 1 from equation (44) for l = 3 and 0 from equation (44) for l = 4 and so on.
Substituting expression , into (63) we obtain: then, the first two terms from expression (64) can be written as follows: or as: is effective potential of electron-electron interaction in a common Gaussian approximation, that is, dropping the terms with n 3 in the exponential of expression (33) we obtain: Finally the effective potential of electron-electron interaction becomes: where the first term corresponds to the approximation of zero sums in x, the second term corresponds to one sum in x.
In particular, in a case of a homogeneous system we have the following effective potential: ) is effective potential of electron-electron interaction in RPA.

Calculation of two-particle correlator in approximation of zero sums in x
As shown above (see ( 55)), two-particle correlator M 0 k 1 ,k 2 (x, −x) (in the approximation of zero sums in x) satisfies the integral equation Proceeding to (q, z)-representation, we obtain: where the following notations are introduced: In a case of modelling the surface potential as potential barrier of finite height and in the approximation of "constant density" (which consists in substituting function |ϕ α (z)| 2 with (1 − p α )θ(−z) + p α θ(z), where p α = sin(2γ α )) for β → ∞ secondorder integral equation (69) transforms into first-order integral equation (boundary transition L → ∞ has been made within the integral) [6]: where Calculation of the function M 0 (q, ν = 0|z 1 , z 2 ) was done in [6], its analytical expression in "constant density" approximation being: Equation ( 72) can be solved analytically, its solution being: where g 0 (q, ν = 0|z 1 , z 2 ) is screened potential in the approximation of "constant density" (g 0 = β S G 0 ), which is the solution of integral equation (65) In [6] it is solved as where Q 1 = q 2 + κ 2 TF − ∆κ, Q 2 = q 2 + ∆κ.
Function M(q, ν = 0|z 1 , z 2 ) describes electron-electron correlations, one of which is in the point z 1 , and the other in z 2 .There being no Coulomb repulsion between electrons M(q, ν = 0|z 1 , z 2 ) (75) has an abrupt minimum -this corresponds to the impossibility to find two electrons in one point.Coulomb interaction taken into account (expression (77)) there occurs additional repulsion, which results in blurring of this abrupt minimum (the electrons feel each other at a distance).
Thus, the problem has been basically solved concerning the calculation of the effective potential in a low temperature case, and concerning the thermodynamic potential based on it, since substituting the found function M 0 1,2 (1, −1) into expansion (67) we have got the effective potential of electron-electron interaction.

Conclusions
A new approach is proposed to calculate the thermodynamic potential, which consists in reducing the relevant non-Gaussian functional integral to its Gaussian form with a renormalized "density-density" correlator.This study is the first to have obtained a general expression for this correlator (42) for the case of an inhomogeneous system.For the case of a homogeneous system, our results coincide with those obtained in [12].It is shown that the knowledge of the effective potential of electron-electron interaction is sufficient to calculate the thermodynamic potential in this approach.
An expression for the effective potential has been obtained in the form of an infinite series and it is sufficient to find a two-particle correlator in approximation of zero sums in x to calculate it.For this correlator an integral equation has been obtained which was solved analytically for a case of low temperatures in "constant density" approximation.Thus, the problem has been basically solved regarding the calculation of the thermodynamic potential in an approximation in which four-particle correlations are taken into account.