Statistical Theory of Two-Component Crystal I-Case of One Dimensional Crystal

Molhem Assef

JSM-Center for Scientific and Engineering Researches, Moscow, Russia.

DOI: 10.4236/oalib.1107433                

Abstract

The Unsymmetrized Self-Consistent Field Method (USCFM) has applied to a linear infinite chain consisting of two different particles. Expanding the self-consistent potential in power series of the lattice displacement, and taking the anharmonicity up to the fifth order, we show that all thermodynamic properties can be found in terms of two universal functions M₁^(α)M₁^(β), which can be implicitly expressed in terms of parabolic cylinder functions. We have applied this approach to study one dimensional system of KCl crystal using Born-Mayer-Huggins pair potential, and have determined some thermoelastic properties of this system.

Keywords

BBGKY, Self-Consistent Theory, Anhormonicity, Alkali Halide Crystals, Thermodynamical Properties, Linear Chain, Expansion Coefficient

1. Introduction

One of the practical methods to study many particles systems is to use many particle distribution functions. To describe the kinetic processes, we should commence from the Liouville’s equation and use the probability density function in the phase space. However, the statistical thermodynamic, which is concerned with the study of equilibrium state of the macroscopic systems, can be studied with the aid of the universal Gibbs ensemble. Indeed, in the two ultimate cases, we obtain an infinite chain of equations for many-particle distribution functions. This chain is called the Bogolyubov-Born-Green-Kirkwood-Yvon (BBGKY) chain of equations. But this system of equations is not closed; the equation for an s-particle distribution function includes s + 1-particle distribution function [1].

To close this system, we usually apply different assumptions, one of the relevant assumptions is based on the principle of spatial attenuation of correlations [2], this is related to the fact that the correlation between particles or groups of particles weakens, if the distance between them increases, and therefore the many-particle distribution function decomposes into a product of distribution functions pertaining to each particle or group of particles, and one can easily obtain a kinetic equation with a self-consistent field, or Vlasov’s equation, which is valid in the case of collisionless low-density warm plasma [3].

In the case of crystalline solid, we can neglect correlations, taking into consideration that atoms (or ions) vibrate in the vicinity of their own lattice points, and the average vibrations’ amplitudes are much smaller than lattice constant, which can expressed by the smallness of the Lindeman’s parameter [4]. Thus, we can represent the many-particle distribution function by a product of individual distribution functions.

The aim of this paper is to develop a method, based on the Unsymmetrized Self-Consistent Field Method (USCFM) [5], to investigate a strongly anharmonic two-component perfect crystal. The first assumption in this approach is the factorizability of the many-particle probability density function, which, as we have argued above, is very reasonable in crystals. The second assumption in the USCFM is the asymmetry of the phase-space probability density with respect to the interchange of the canonical coordinates $\left({\stackrel{\to }{r}}_i,{\stackrel{\to }{p}}_i\right)$ between similar particles.

Taking the anharmonicity up to the fifth order, we have demonstrated that the amplitudes of atomic vibrations about lattice sites can be expressed in term of two universal functions $M_1^{\left(\alpha \right)},M_1^{\left(\beta \right)}$ , which can be implicitly expressed in terms of some famous parabolic cylinder functions; and we show, that with the aid of these two functions, one can calculate Helmholtz free energy and then reproduce all the thermodynamic of crystal.

We have applied this method to linear chain potassium chloride crystal, and have calculated some thermodynamical properties using Born-Mayer-Huggins pair potential. Forthcoming papers are devoted to two- and three-dimensional crystals.

2. Basic Equations

The single self-consistent potential for the particle $\mu$ is given by [5]:

$\begin{array}{l}{u}_{\mu }\left({\stackrel{\to }{q}}_{\mu }\right)={\displaystyle \underset{n\ne 0}{\sum }{\displaystyle \int f}\left({\stackrel{\to }{q}}_{\mu }-{{\stackrel{\to }{q}}^{\prime }}_{\nu }-\hat{A}\stackrel{\to }{n}\right)\omega \left({{\stackrel{\to }{q}}^{\prime }}_{\nu }\right)\text{d}{{\stackrel{\to }{q}}^{\prime }}_{\nu }}\\ -{\displaystyle \underset{n\ne 0}{\sum }{\displaystyle \int {\displaystyle \int f\left({\stackrel{\to }{q}}_{\mu }-{{\stackrel{\to }{q}}^{\prime }}_{\nu }-\hat{A}\stackrel{\to }{n}\right)\omega \left({\stackrel{\to }{q}}_{\mu }\right)\omega \left({{\stackrel{\to }{q}}^{\prime }}_{\nu }\right)\text{d}{\stackrel{\to }{q}}_{\mu }\text{d}{{\stackrel{\to }{q}}^{\prime }}_{\nu }}}}\end{array}$ (1)

The summation over the different neighbors can be done before any other calculations by introducing the functions:

$K^{\mu \nu }\left(x\right)={\displaystyle \underset{n\ne 0}{\sum }f\left({\stackrel{\to }{x}}_{\mu \nu }-\hat{A}\stackrel{\to }{n}\right)}$ (2)

$K^{\mu \nu }\left({\stackrel{\to }{q}}_{\mu }-{{\stackrel{\to }{q}}^{\prime }}_{\nu }\right)={\displaystyle \underset{n\ne 0}{\sum }f\left({\stackrel{\to }{q}}_{\mu }-{{\stackrel{\to }{q}}^{\prime }}_{\nu }-\hat{A}\stackrel{\to }{n}\right)}$ (3)

$\begin{array}{l}{u}_{\mu }\left({\stackrel{\to }{q}}_{\mu }\right)={\displaystyle \underset{\nu }{\sum }{K}^{\mu \nu }\left({\stackrel{\to }{q}}_{\mu }-{{\stackrel{\to }{q}}^{\prime }}_{\nu }\right)\omega \left({{\stackrel{\to }{q}}^{\prime }}_{\nu }\right)\text{d}{{\stackrel{\to }{q}}^{\prime }}_{\nu }}\\ -{\displaystyle \underset{\nu }{\sum }{\displaystyle \int {\displaystyle \int K^{\mu \nu }\left({\stackrel{\to }{q}}_{\mu }-{{\stackrel{\to }{q}}^{\prime }}_{\nu }\right)\omega \left({\stackrel{\to }{q}}_{\mu }\right)\omega \left({{\stackrel{\to }{q}}^{\prime }}_{\nu }\right)\text{d}{\stackrel{\to }{q}}_{\mu }\text{d}{{\stackrel{\to }{q}}^{\prime }}_{\nu }}}}\end{array}$ (4)

where ${\stackrel{\to }{q}}_{\mu },{{\stackrel{\to }{q}}^{\prime }}_{\nu }$ are the lattice displacements, $\hat{A}$ is the lattice matrix, $\stackrel{\to }{n}$ is vector having integer components, $\omega \left({\stackrel{\to }{q}}_{\mu }\right)$ is the probability of presence of μ-th particle.

On the other hand, we accept the following expression for the probability of presence at temperature $\theta =kT$ [5]:

$\omega \left({\stackrel{\to }{q}}_{\mu }\right)=\frac{\exp \left[-u_{\mu }\left({\stackrel{\to }{q}}_{\mu }\right)/\theta \right]}{z^{\mu }\left(\theta ,\stackrel{\to }{a}\right)}$ (5)

where $z^{\mu }\left(\theta ,\stackrel{\to }{a}\right)$ is the normalization constant, which is really the partial partition function for single ion:

$z^{\mu }\left(\theta ,\stackrel{\to }{a}\right)={\displaystyle \int \exp \left[-u_{\mu }\left({\stackrel{\to }{q}}_{\mu }\right)/\theta \right]\text{d}{\stackrel{\to }{q}}_{\mu }}$ (6)

In this approximation, the thermodynamic functions of crystal are based on the individual self-consistent potentials, and we will use a zero order partition function $Z_0\left(\theta ,\stackrel{\to }{a}\right)$ :

$z^{\mu }\left(\theta ,\stackrel{\to }{a}\right)={\displaystyle \int \exp \left[-u_{\mu }\left({\stackrel{\to }{q}}_{\mu }\right)/\theta \right]\text{d}{\stackrel{\to }{q}}_{\mu }}$ (7)

where all the partial partition function $Z_i\left(\theta ,\stackrel{\to }{a}\right)$ are equal:

$Z_i\left(\theta ,\stackrel{\to }{a}\right)=Z_{\mu }\left(\theta ,\stackrel{\to }{a}\right)=h^{-d}{\displaystyle \iint \exp \left\{-\left[{\stackrel{\to }{p}}_{\mu }^2/2m+u_{\mu }\left({\stackrel{\to }{q}}_{\mu }\right)\right]/\theta \right\}\text{d}{\stackrel{\to }{p}}_{\mu }\text{d}{\stackrel{\to }{q}}_{\mu }}$ (8)

h being the Planck constant, and d dimension of the system.

The kinetic part can be integrated out:

$Z_{\mu }\left(\theta ,\stackrel{\to }{a}\right)={\left(\frac{2\pi m_{\mu }\theta }{h^2}\right)}^{d/2}{\displaystyle \int \exp \left\{-u_{\mu }\left({\stackrel{\to }{q}}_{\mu }\right)/\theta \right\}\text{d}{\stackrel{\to }{q}}_{\mu }}={\left(\frac{2\pi m_{\mu }\theta }{h^2}\right)}^{d/2}{z}_{\mu }\left(\theta ,\stackrel{\to }{a}\right)$ (9)

3. One Dimensional System

Let us expand the self-consistent potential of the particle μ up to the fifth order:

$u_{\mu }\left(q_{\mu }\right)={\displaystyle \underset{l=0}{\overset{5}{\sum }}{u}_l\frac{q_{\mu }^l}{l!}}$ (10)

Similarly, we can also expand the functions $K^{\mu \nu }\left(q_{\mu }-{q^{\prime }}_{\nu }\right)$ :

$K^{\mu \nu }\left(q_{\mu }-{q^{\prime }}_{\nu }\right)={\displaystyle \underset{j=0}{\overset{5}{\sum }}{K}_j^{\mu \nu }\frac{{\left(q_{\mu }-{q^{\prime }}_{\nu }\right)}^j}{j!}}$ (11)

where

$K_j={\frac{{\partial }^{j}K\left(q\right)}{\partial q^j}|}_{q=0}$

Hence, the self-consistent potential, in the case of linear chain consisting of two different particles $\alpha ,\beta$ (cation and anion in the case of alkali halides crystals, for example), can be written in the form:

$u_{\mu }\left(q_{\mu }\right)={\displaystyle \underset{\nu =\alpha ,\beta }{\sum }\left\{{\displaystyle \underset{j=0}{\overset{5}{\sum }}\left[K_j^{\mu \nu }{\displaystyle \underset{k=0}{\overset{j}{\sum }}\frac{{\left(-1\right)}^k}{k!\left(j-k\right)!}{q}_{\mu }^{j-k}}\right]{\displaystyle \int {\left({q^{\prime }}_{\nu }\right)}^k\omega \left({q^{\prime }}_{\nu }\right)\text{d}{q^{\prime }}_{\nu }}}\right\}}-u_{00}$ (12)

$u_{00}=\frac{1}{2}{\displaystyle \underset{\nu =\alpha ,\beta }{\sum }{\displaystyle \int {\displaystyle \underset{j=0}{\overset{5}{\sum }}{K}_j^{\mu \nu }{\displaystyle \underset{k=0}{\overset{j}{\sum }}\frac{{\left(-1\right)}^k}{k!\left(j-k\right)!}{q}_{\mu }^{j-k}}{\left({q^{\prime }}_{\nu }\right)}^k\omega \left(q_{\mu }\right)\omega \left({q^{\prime }}_{\nu }\right)\text{d}{q}_{\mu }\text{d}{q^{\prime }}_{\nu }}}}$ (13)

Comparing the relation (12) with (11) and considering the case of alkali halide crystal, one can easily get:

$u_l^{\mu }={\displaystyle \underset{m=0}{\overset{5-l}{\sum }}{\left(-1\right)}^m\left[K_{l+m}^{\mu n}\frac{\langle q_n^m\rangle }{m!}+K_{l+m}^{\mu p}\frac{\langle q_p^m\rangle }{m!}\right]}\text{for}l>0$ (14)

$\begin{array}{l}{u}_0^{\mu }={\displaystyle \underset{m=0}{\overset{5}{\sum }}{\left(-1\right)}^m\left[K_{l+m}^{\mu n}\frac{\langle q_n^m\rangle }{m!}+K_{l+m}^{\mu p}\frac{\langle q_p^m\rangle }{m!}\right]}\\ -\frac{1}{2}{\displaystyle \underset{j=0}{\overset{5}{\sum }}{\displaystyle \underset{k=0}{\overset{j}{\sum }}{\left(-1\right)}^k\left(K_j^{\mu n}\frac{\langle q_n^k\rangle }{k!}+K_j^{\mu p}\frac{\langle q_p^k\rangle }{k!}\right)\frac{\langle q_{\mu }^{j-k}\rangle }{\left(j-k\right)!}}}\end{array}$ (15)

where the term $\langle q_{\mu }^k\rangle$ express the presence probability moment of k-th order for particle $\mu$ :

$\langle q_{\mu }^k\rangle ={\displaystyle \int {\left(q_{\mu }\right)}^k\omega \left(q_{\mu }\right)\text{d}{q}_{\mu }}$ (16)

p denote the positive ions (alkali ions), and n the negative ions (halogen ions).

We can restrict the expansion of $u_{\mu }\left(q_{\mu }\right)$ to the fourth order, in the case of an isotropic interaction, as even terms of K are then zero; in that case $\langle q_{\mu }^k\rangle$ and $u_k^{\mu }\left(q_{\mu }\right)$ are also zero for odd’s k.

$u_{\mu }\left(q_{\mu }\right)=u_0^{\mu }+u_2^{\mu }\frac{q_{\mu }^2}{2}+u_4^{\mu }\frac{q_{\mu }^4}{24}$ (17)

One can easily find:

$\begin{array}{l}{u}_0^{\mu }=\frac{1}{2}\left(K_0^{\mu n}+K_0^{\mu p}\right)+\frac{1}{4}\left(K_2^{\mu n}\langle q_n^2\rangle +K_2^{\mu p}\langle q_p^2\rangle \right)\\ +\frac{1}{48}\left(K_4^{\mu n}\langle q_n^4\rangle +K_4^{\mu p}\langle q_p^4\rangle \right)-\frac{\langle q_{\mu }^2\rangle }{4}\left(K_2^{\mu n}+K_2^{\mu p}\right)\\ -\frac{\langle q_{\mu }^4\rangle }{48}\left(K_4^{\mu n}+K_4^{\mu p}\right)-\frac{\langle q_{\mu }^2\rangle }{8}\left(K_4^{\mu n}\langle q_n^2\rangle +K_4^{\mu p}\langle q_p^2\rangle \right)\end{array}$ (18)

$u_2^{\mu }=K_2^{\mu n}+K_2^{\mu p}+K_4^{\mu n}\frac{\langle q_n^2\rangle }{2}+K_4^{\mu p}\frac{\langle q_p^2\rangle }{2}$ (19)

$u_4^{\mu }=K_4^{\mu n}+K_4^{\mu p}$ (20)

Let us consider the following integral:

$I_{2p}=2{\displaystyle \underset{0}{\overset{\infty }{\int }}{q}^{2p}\exp \left(\frac{-u\left(q\right)}{\theta }\right)\text{d}q}$ (21)

It is easily to verify that:

$I_{2p}={\left(\sqrt{\frac{12\theta }{u_4}}\right)}^{p+\frac{1}{2}}{\text{e}}^{-\frac{u_0}{\theta }}{\displaystyle \underset{0}{\overset{\infty }{\int }}{x}^{p-\frac{1}{2}}\exp \left(-\frac{x^2}{2}-yx\right)\text{d}x}$ (22)

$y=\frac{u_2}{2\theta }\sqrt{\frac{12\theta }{u_4}}=u_2\sqrt{\frac{3}{\theta u_4}}$ (23)

Using the expression of the transcendental parabolic cylinder function [6]:

$D_n\left(y\right)\Gamma \left(-n\right){\text{e}}^{\frac{y^2}{4}}={\displaystyle \underset{0}{\overset{\infty }{\int }}{x}^{-n-\frac{1}{2}}\exp \left(-\frac{x^2}{2}-yx\right)\text{d}x}$ (24)

One can find easily:

$I_{2p}={\left(\sqrt{\frac{12\theta }{u_4}}\right)}^{p+\frac{1}{2}}{\text{e}}^{-\frac{u_0}{\theta }}{D}_{-p-\frac{1}{2}}\left(y\right)\Gamma \left(p+\frac{1}{2}\right){\text{e}}^{\frac{y^2}{4}}$ (25)

Hence the normalized moment of order 2:

$\langle q^2\rangle =\frac{I_2}{I_0}$ (26)

$\sqrt{\frac{u_4}{3\theta }}\langle q^2\rangle =\frac{D_{-1.5}\left(y\right)}{D_{-0.5}\left(y\right)}$ (27)

Then

$M\left(y\right)=\sqrt{\frac{u_4}{3\theta }}\langle q^2\rangle$ (28)

so

$y^{\mu }=\left(K_2^{\mu n}+K_2^{\mu p}+K_4^{\mu n}\frac{\langle q_n^2\rangle }{2}+K_4^{\mu p}\frac{\langle q_p^2\rangle }{2}\right)\sqrt{\frac{3}{\theta u_4}}$ (29)

One can easily find:

$y^{\mu }=X^{\mu }+\frac{3}{2}\frac{K_4^{\mu n}}{\sqrt{K_4^{{}^{\mu }}{K}_4^n}}{M}_1^{\left(n\right)}\left(y^n\right)+\frac{3}{2}\frac{K_4^{\mu p}}{\sqrt{K_4^{{}^{\mu }}{K}_4^p}}{M}_1^{\left(p\right)}\left(y^p\right)$ (30)

$y^n=X^n+\frac{3}{2}\frac{K_4^{nn}}{K_4^n}{M}_1^{\left(n\right)}\left(y^n\right)+\frac{3}{2}\frac{K_4^{np}}{\sqrt{K_4^n{K}_4^p}}{M}_1^{\left(p\right)}\left(y^p\right)$ (31)

$y^p=X^p+\frac{3}{2}\frac{K_4^{pn}}{\sqrt{K_4^p{K}_4^n}}{M}_1^{\left(n\right)}\left(y^n\right)+\frac{3}{2}\frac{K_4^{pp}}{K_4^p}{M}_1^{\left(p\right)}\left(y^p\right)$ (32)

where

$\begin{array}{l}{K}_4^{\mu }=K_4^{\mu n}+K_4^{\mu p}\\ K_2^{\mu }=K_2^{\mu n}+K_2^{\mu p}\\ K_0^{\mu }=K_0^{\mu n}+K_0^{\mu p}\end{array}$

$X^{\mu }=K_2^{\mu }\sqrt{\frac{3}{\theta K_4^{\mu }}}$

Let us denote:

$\begin{array}{l}{x}_1=K_2^n\sqrt{\frac{3}{\theta K_4^n}},x_2=K_2^p\sqrt{\frac{3}{\theta K_4^p}},x_3=\frac{3}{2}\frac{K_4^{nn}}{K_4^n},\\ x_4=\frac{3}{2}\frac{K_4^{pp}}{K_4^p},x_5=\frac{3}{2}\frac{K_4^{pn}}{\sqrt{K_4^p{K}_4^n}}\end{array}$

As we have

$M^{\mu }=\frac{D_{-1.5}\left(y^{\mu }\right)}{D_{-0.5}\left(y^{\mu }\right)}$ (33)

We can deduce the two following implicit equations:

$M_1^{\left(n\right)}=\frac{D_{-1.5}\left(x_1+x_3{M}_1^{\left(n\right)}+x_5{M}_1^{\left(p\right)}\right)}{D_{-0.5}\left(x_1+x_3{M}_1^{\left(n\right)}+x_5{M}_1^{\left(p\right)}\right)}$ (34)

$M_1^{\left(p\right)}=\frac{D_{-1.5}\left(x_2+x_4{M}_1^{\left(p\right)}+x_5{M}_1^{\left(n\right)}\right)}{D_{-0.5}\left(x_2+x_4{M}_1^{\left(p\right)}+x_5{M}_1^{\left(n\right)}\right)}$ (35)

Notice that, the two universal functions $M_1^{\left(p\right)},M_1^{\left(n\right)}$ have five variables, and these variables depend on lattice constant and temperature.

One can easily verify that:

$\langle q_{\mu }^{2p}\rangle =\frac{I_{2p}}{I_0}={\left(\sqrt{\frac{12\theta }{u_4^{\mu }}}\right)}^p\frac{\Gamma \left(p+\frac{1}{2}\right)}{\Gamma \left(\frac{1}{2}\right)}\frac{D_{-p-\frac{1}{2}}\left(y_{\mu }\right)}{D_{-\frac{1}{2}}\left(y_{\mu }\right)}$ (36)

Using the following recursion relation on the parabolic cylinder function D [6]:

$D_{n+1}\left(y\right)-y{D}_n\left(y\right)+n{D}_{n-1}\left(y\right)=0$ (37)

One can easily the following recursion relation for the other moments:

$\begin{array}{c}\langle q_{\mu }^{2p}\rangle =\frac{I_{2p}}{I_0}\\ =\sqrt{\frac{12\theta }{u_4^{\mu }}}\left\{\left(2p-3\right)\sqrt{\frac{3\theta }{u_4^{\mu }}}\langle q_{\mu }^{2\left(p-2\right)}\rangle -y_{\mu }\langle q_{\mu }^{2\left(p-1\right)}\rangle \right\}\end{array}$ (38)

Then, the partial partition function of one unit cell, consisting of two ions, is:

$Z=Z^{\left(n\right)}{Z}^{\left(p\right)}$ (39)

$Z^{\left(n\right)}=\sqrt{\frac{2\pi m^{\left(n\right)}\theta }{h^2}}{\left(\frac{12\theta }{u_4^{\left(n\right)}}\right)}^{\frac{1}{4}}{\text{e}}^{-\frac{u_0^{\left(n\right)}}{\theta }}{D}_{-\frac{1}{2}}\left(y^{\left(n\right)}\right)\Gamma \left(\frac{1}{2}\right){\text{e}}^{\frac{{\left(y^{\left(n\right)}\right)}^2}{4}}$ (40)

$Z^{\left(p\right)}=\sqrt{\frac{2\pi m^{\left(p\right)}\theta }{h^2}}{\left(\frac{12\theta }{u_4^{\left(p\right)}}\right)}^{\frac{1}{4}}{\text{e}}^{-\frac{u_0^{\left(p\right)}}{\theta }}{D}_{-\frac{1}{2}}\left(y^{\left(p\right)}\right)\Gamma \left(\frac{1}{2}\right){\text{e}}^{\frac{{\left(y^{\left(p\right)}\right)}^2}{4}}$ (41)

$\begin{array}{l}Z=\sqrt{\frac{2\pi m^{\left(n\right)}\theta }{h^2}}\sqrt{\frac{2\pi m^{\left(p\right)}\theta }{h^2}}{\left(\frac{12\theta {\pi }^2}{u_4^{\left(n\right)}}\right)}^{\frac{1}{4}}{\left(\frac{12\theta {\pi }^2}{u_4^{\left(p\right)}}\right)}^{\frac{1}{4}}\\ \cdot {\text{e}}^{-\frac{u_0^{\left(n\right)}+u_0^{\left(p\right)}}{\theta }}{D}_{-\frac{1}{2}}\left(y^{\left(n\right)}\right)D_{-\frac{1}{2}}\left(y^{\left(p\right)}\right){\text{e}}^{\frac{{\left(y^{\left(n\right)}\right)}^2+{\left(y^{\left(p\right)}\right)}^2}{4}}\end{array}$ (42)

Hence, the Helmholtz free energy $F_0$ , in this approximation, is [7]:

$F_0=-\theta \ln \left[Z_0\left(\theta ,\stackrel{\to }{a}\right)\right]=-\theta N\ln \left[Z\left(\theta ,\stackrel{\to }{a}\right)\right]$ (43)

Hence, the free energy in zeroth order approximation:

$\begin{array}{l}Z=N\theta \{-\frac{1}{2}\ln \frac{2\pi m^{\left(n\right)}\theta }{h^2}-\frac{1}{2}\ln \frac{2\pi m^{\left(p\right)}\theta }{h^2}-\frac{1}{4}\ln \frac{12\theta {\pi }^2}{u_4^{\left(n\right)}}-\frac{1}{4}\ln \frac{12\theta {\pi }^2}{u_4^{\left(p\right)}}\\ +\frac{u_0^{\left(n\right)}+u_0^{\left(p\right)}}{\theta }-\ln D_{-\frac{1}{2}}\left(y^{\left(n\right)}\right)-\ln D_{-\frac{1}{2}}\left(y^{\left(p\right)}\right)-\frac{{\left(y^{\left(n\right)}\right)}^2+{\left(y^{\left(p\right)}\right)}^2}{4}\}\end{array}$ (44)

Notice that, $y^{\left(n\right)},y^{\left(p\right)}$ are function of $M_1^{\left(p\right)},M_1^{\left(n\right)}$ . Evidently, one can easily deduce all the thermodynamical properties of the system, using free energy.

We can introduce the “linear pressure” $A_0$ as the average conjugate of the length $Na$ of chain:

$A_0=\theta \frac{\partial \ln \left[Z\left(\theta ,\stackrel{\to }{a}\right)\right]}{\partial a}=\theta \frac{\partial \ln \left[Z^{\left(n\right)}{Z}^{\left(p\right)}\right]}{\partial a}$ (45)

Hence,

$\begin{array}{l}{A}_0=\theta \{-\frac{1}{4}\frac{1}{u_4^{\left(n\right)}}\frac{\partial u_4^{\left(n\right)}}{\partial a}-\frac{1}{4}\frac{1}{u_4^{\left(p\right)}}\frac{\partial u_4^{\left(p\right)}}{\partial a}-\frac{1}{\theta }\left(\frac{\partial u_0^{\left(n\right)}}{\partial a}+\frac{\partial u_0^{\left(p\right)}}{\partial a}\right)\\ -\frac{1}{2}\left(M_1^{\left(n\right)}\frac{\partial y^{\left(n\right)}}{\partial a}+M_1^{\left(p\right)}\frac{\partial y^{\left(p\right)}}{\partial a}\right)\}\end{array}$ (46)

$\begin{array}{l}{A}_0=-\theta \{\frac{1}{4}\frac{1}{u_4^{\left(n\right)}}\frac{\partial u_4^{\left(n\right)}}{\partial a}+\frac{1}{4}\frac{1}{u_4^{\left(p\right)}}\frac{\partial u_4^{\left(p\right)}}{\partial a}+\frac{1}{2\theta }\frac{\partial }{\partial a}\left(K_0^n+K_0^p\right)+\frac{1}{2}{M}_1^{\left(n\right)}\frac{\partial x_1}{\partial a}\\ +M_1^{\left(p\right)}\frac{\partial x_2}{\partial a}+\frac{1}{2}{M}_1^{\left(n\right)}{M}_1^{\left(p\right)}\frac{\partial x_5}{\partial a}+\frac{1}{4}{M}_1^{\left(n\right)2}\frac{\partial x_3}{\partial a}+\frac{1}{4}{M}_1^{\left(p\right)2}\frac{\partial x_4}{\partial a}\}\end{array}$ (47)

The main problem is to find the values of $M_1^{\left(p\right)},M_1^{\left(n\right)}$ for every set of value of the five variables $x_1,x_1,x_3,x_4,x_5$ . Fortunately, the two Equations (34) and (35) are written in appropriate form, such that we can easily apply the method of iterations; but the principal difficulty in this approach is to express the function $D_{\nu }$ through some elementary functions to simplify the calculations process; one can circumvent this problem using the hyper-geometrical function $\Phi \left(a,c;x\right)$ , thus we can write [6]:

$D_{\nu }\left(z\right)=2^{\frac{\nu }{2}}{\text{e}}^{-\frac{z^2}{4}}[\frac{\Gamma \left(\frac{1}{2}\right)}{\Gamma \left(\frac{1-\nu }{2}\right)}\Phi \left(-\frac{\nu }{2},\frac{1}{2};\frac{z^2}{2}\right)+\frac{z}{\sqrt{2}}\frac{\Gamma \left(-\frac{1}{2}\right)}{\Gamma \left(-\frac{\nu }{2}\right)}\Phi \left(\frac{1-\nu }{2},\frac{3}{2};\frac{z^2}{2}\right)]$ (48)

where

$\Phi \left(a,c;x\right)=\frac{\Gamma \left(c\right)}{\Gamma \left(a\right)\Gamma \left(c-a\right)}{\displaystyle \underset{0}{\overset{1}{\int }}{\text{e}}^{xu}{u}^{a-1}{\left(1-u\right)}^{c-a-1}\text{d}u}$ (49)

One can easily find for:

$\frac{D_{-1.5}\left(z\right)}{D_{-0.5}\left(z\right)}=\frac{1}{\sqrt{2}}\frac{{\beta }_1\Phi \left(\frac{3}{4},\frac{1}{2};\frac{z^2}{2}\right)+{\beta }_2\Phi \left(\frac{5}{4},\frac{3}{2};\frac{z^2}{2}\right)}{{\beta }_3\Phi \left(\frac{3}{4},\frac{1}{2};\frac{z^2}{2}\right)+{\beta }_4\Phi \left(\frac{5}{4},\frac{3}{2};\frac{z^2}{2}\right)}$ (50)

where

${\beta }_1=\frac{\Gamma \left(\frac{1}{2}\right)}{\Gamma \left(\frac{5}{4}\right)},{\beta }_2=\frac{\Gamma \left(-\frac{1}{2}\right)}{\sqrt{2}\Gamma \left(\frac{3}{4}\right)},{\beta }_3=\frac{\Gamma \left(\frac{1}{2}\right)}{\Gamma \left(\frac{3}{4}\right)},{\beta }_4=\frac{\Gamma \left(-\frac{1}{2}\right)}{\sqrt{2}\Gamma (\; 1\; 4\; )}$

The functions $\Phi \left(a,c;x\right)$ are very well studied and tabulated in many monographs (see for example [6] ). Hence, to apply the method of iterations, we can accept some initial values of $M_1^{\left(p\right)},M_1^{\left(n\right)}$ and after several iterations, the difference between the left side and the right side of the Equations ((34), (35)) must be zero.

It should be mentioned here that the variables $x_1,x_1,x_3,x_4,x_5$ are functions of lattice constant and temperature, so someone can say it is enough to express the universal functions $M_1^{\left(p\right)},M_1^{\left(n\right)}$ in term of $\left(a,T\right)$ , but to apply this approach we should determine the form of crystal and interionic potential. However the universality of our approach appears in the use of $M_1^{\left(p\right)},M_1^{\left(n\right)}$ , the main challenge is to express them as functions of $x_1,x_1,x_3,x_4,x_5$ .

In the case of weak anharmonicity, i.e. when:

$K_2^n\sqrt{\frac{3}{\theta K_4^n}},K_2^p\sqrt{\frac{3}{\theta K_4^p}}\gg 1$

we can expand the function $D_{\nu }\left(z\right)$ as following

$D_{\nu }\left(z\right)=z^{\nu }{\text{e}}^{-\frac{z^2}{4}}\left[1-\frac{\nu \left(\nu -1\right)}{2z^2}+\frac{\nu \left(\nu -1\right)\left(\nu -2\right)\left(\nu -3\right)}{8z^4}+\cdots \right]$ (51)

restricting at the first term of the expansion, one can obtain:

$\frac{D_{-1.5}\left(z\right)}{D_{-0.5}\left(z\right)}=\frac{1}{z}$ (52)

Hence equations (34) and (35) can easily be solved algebraically; and knowing that $M_1^{\left(p\right)},M_1^{\left(n\right)}$ are always positive, we can obtain:

$M_1^{\left(p\right)}=0.5\left(-x_2+\sqrt{x_2^2+4x_4+4x_5}\right)/\left(x_4+x_5\right)$ (53)

$\begin{array}{l}{M}_1^{\left(n\right)}=\frac{1}{4x_3\left(x_4+x_5\right)}[-2x_1{x}_4-2x_1{x}_5+x_2{x}_5-x_5\sqrt{x_2^2+4x_4+4x_5}\\ +\sqrt{16x_3{\left(x_4+x_5\right)}^2+{\left(2x_1{x}_4+2x_1{x}_5-x_2{x}_5+x_5\sqrt{x_2^2+4x_4+4x_5}\right)}^2}]\end{array}$ (54)

These results enable us to investigate some aspects of the problem and prove the good usability of our theory to obtain thermodynamical properties of two-component crystals.

4. Results and Discussion

In the present calculations the Born-Mayer-Huggins potential was used to describe the forces between two ions [8]:

${\Phi }^{ij}\left(r\right)=\frac{e_i{e}_j}{4\pi {\epsilon }_{0}r}+b{\beta }_{ij}\exp \left(\left(R_i+R_j-r\right)/\rho \right)-\frac{c_{ij}}{r^6}-\frac{d_{ij}}{r^8}$ (55)

where the factors ${\beta }_{ij}$ originate from certain quantum-mechanical considerations and they depend only on the charges of interacting ions, namely: ${\beta }_{pp}=1.25,{\beta }_{pn}=1,{\beta }_{nn}=0.75$ . The other parameters are listed in Table 1 and Table 2.

Concerning Table 1 and Table 2, it should be mentioned here that $c_{ij},d_{ij}$ are Van der Waals coefficients characteristics of the dipole-dipole and dipole-quadrupole interactions respectively, and its value are taken from [9], the values of ionic radii $R_n,R_p$ represents Pauling set of radii and taken from [10], while the parameters $b,\rho$ are taken from [8]. In Table 2, parameters (Set. 1) has been computed using static crystal approach and experimental data at 0K; (Set. 2) has been computed using Hildebrand equation of state and so called first set of experimental data at room temperature; (Set. 3) has been computed using Hildebrand equation of state and so called second set of experimental data at room temperature; (Set. 4) has been computed using Mie-Gruniesen equation of state and so called first set of experimental data at room temperature; (Set. 5) has been computed using Mie-Gruniesen equation of state and so called second set of experimental data at room temperature (For details see our work [8] ).

Table 1. The parameters of pair potential (Equation (55)).

Table 2. The parameters of pair potential (Equation (55)).

The interionic distance can be obtained by equating the expression of the linear pressure to a known external pressure, a good approximation in many cases:

$A_0=0$ (56)

Solving numerically the Equation (56), one can obtain the nearest-neighbor distance $a\left(T\right)$ at various temperatures; the results of our calculations are re-presented in Figure 1 and Figure 2.

Figure 1. Solutions of Equation (56) for nearest-neighbor distance of KCl linear chain, as function of temperature.

Figure 2. Nearest-neighbor distance as function of temperature of 1D KCl crystal.

One can easily notice that various set of potential parameters give close results for the nearest neighbor distance at temperature ranges.

In all cases, Equation (56) has two roots $a_1\left(T\right)\le a_2\left(T\right)$ , these two roots coincide at temperature Ts (spinodal point), 4 times higher than melting temperature, the lower branch of the isobar represents the stable thermodynamic states, while the upper branch represents the unstable ones, at the spinodal point of the isobar ${\left(\partial a/\partial T\right)}_T\to \infty$ .

So, along the lower branch of the isobar, one easily can calculate all the thermodynamic properties of linear chain of potassium chloride. For example, we have calculated the linear thermal expansion coefficient along the isobar using the following formula:

$\alpha =\frac{1}{a}\frac{\partial a}{\partial T}$ (57)

and the results are represented in Figure 3.

Unfortunately, there are no experimental data about one-dimensional KCl crystal to compare with our theoretical results; and in Figure 3 we have presented experimental data of real three-dimensional KCl crystal to demonstrate the order of thermal expansion coefficient; so it is not reasonable to treat the discrepancy between theoretical results and experimental data.

Additionally, in the present work, all calculations have been achieved at zeroth order approximation, i.e. when the probability density function is factorizable. So, to obtain a real picture of the crystal, all thermodynamic functions should be corrected:

Figure 3. Linear expansion coefficient KCl crystal. [11] [12] [13].

1) At low temperature, where quantum effects can arise and we can include it using quantum perturbation theory [5] [7];

2) At high temperature, where correlative effects can take place and we can include it using Bogoliubov’s statistical principle [14].

Furthermore, that the parameters of the Born-Mayer-Huggins potential have been determined using experimental data of real crystal, i.e. three dimensional crystal, so the one-dimensional alkali-halide crystals can’t be exactly described using the same parameters as in the three-dimensional crystal, because the environment of the ions are different in the two cases.

Moreover, strong anharmonicity should be taken into consideration to correctly describe some properties of crystals.

All these problems will be studied in future works.

5. Conclusions

Basing on the Unsymmetrized Self-Consistent Field method, we have developed a statistical theory for two-component crystal. The USCFM approximation does not take account of any correlation or quantum effects, which can be included as corrections.

We have found that, all thermodynamical functions can be expressed through two universal functions $M_1^{\left(p\right)},M_1^{\left(n\right)}$ , which can be implicitly expressed in terms of parabolic cylinder function $D_m$ (Webber function). We show that the functions $M_1^{\left(p\right)},M_1^{\left(n\right)}$ can be calculated algebraically using the hypergeometrical function $\Phi \left(a,c;x\right)$ .

In the case of weak anharmonicity, we have calculated the interionic distance of KCl linear chain using the equation of state $A_0=0$ , as function of temperature, and we have deduced the linear thermal expansion coefficient along the isobar.

In summary, the use of functions $M_1^{\left(p\right)},M_1^{\left(n\right)}$ enables us to calculate exactly all the thermodynamic properties of crystal, which can be considered as promised and perspective approach to study some relevant problems as: surfaces, adsorption, impurities and defects, solid mixtures.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Diu, B., Guthman, C., Lederer, D. and Roulet B. (1989) Elements De Physique Statistique. Hermann, Paris. (in French)
[2]	Akhiezer, A.I. and Peletminskii, S.V. (1981) Methods of Statistical Physics. Translated from Russian by M. Schukin, Pergamon Press, Oxford.
[3]	Akhiezer, A.I. (1974) Electrodynamics of Plasma. Nauka, Moscow. (in Russian)
[4]	Bazarov, I.P. (1971) Statistical Theory of Crystalline State. Moscow State University, Moscow. (in Russian).
[5]	Zubov, V.I. (2003) Unsymmetrized Distribution Functions and the Self-Consistent Theory of Strongly Anharmonic Crystals. Vestnik RUDN, Seria Physica, No. 11, 119-141. (in Russian)
[6]	Bateman, H. and Erdelyi, A. (1953) Higher Transcendental Functions: Vol. II. McGraw-Hill Book Company, Inc., New York.
[7]	Landau, L.D. and Lifshitz, E.M. (1950) Statistical Physics. Pergamon Press, London.
[8]	Molhem, A. (2021) Interionics in Alkali Halide Crystals. (in Print)
[9]	Mayer, J.E. (1933) Dispersion and Polarizability and the van der Waals Potential in the Alkali Halides. The Journal of Chemical Physics, 1, 270. https://doi.org/10.1063/1.1749283
[10]	Born, M. and Huang, K. (1954) Dynamical Theory of Crystal Lattices. Clarendon Press, Oxford.
[11]	Anderson, O.L. (1995) Equation of State of Solids for Geophysics and Ceramic Science. Oxford University Press, Oxford.
[12]	Sirdeshmukh, D.B., Sirdeshmukh, L. and Subhadra, K.G. (2001) Alkali Halides: A Handbook of Physical Properties. Springer-Verlag, Berlin. https://doi.org/10.1007/978-3-662-04341-7
[13]	Salimäki, K.E. (1960) The Thermal Expansion Coefficients of KCl and KBr, and Their Solid Solutions in the Temperature Range -180° to 400°C. Annales Academiae Scientiarum Fennicae. Series A. VI, Physica, 56, 1-40. (in German).
[14]	Bogoliubov, N.N. (1946) The Problems of Dynamical Theory in Statistical Physics. GosTekhIzdat, Moscow-Leningrad. (in Russian).