Simplified Quantum Statistical Model for Electron Equation of State of Materials

1. Introduction

Orbital-free density functional theories of electronic structure provide the simplest way for determining the electron-components of thermodynamic properties of materials [1]. The Thomas-Fermi (TF) model with its finite-temperature generalization [2], together with the spherical atomic cell concept, is the workhorse for computing electron-properties of hot dense matter [3,4]. The excellent compilation of approximate solutions of TF equation [5] shows the need for efficient and accurate numerical methods in this field. Addition of the quantum mechanical exchange correction to the TF model leads to Thomas-Fermi-Dirac (TFD) model [6], while the Thomas-Fermi-Dirac-Weizsäcker (TFDW) model, also known as quantum statistical model (QSM), is obtained when electron-density gradient corrections to kinetic energy are also accounted [7]. Numerical computations of electron-component of equation of state (EOS) using the QSM is somewhat involved [8,9]. Moreover, there is no scaling of thermodynamic properties with atomic (Z) and mass (A) numbers [10,11]. Therefore, it would be valuable to have improvements to the TF model, with the addition of electron exchange-correlation and gradient effects, however, at less computational efforts.

In a seminal work, Schwinger showed [12] that bound electron contribution to the energy of a cold free atom within the TF model is divergent, and the addition of the correct component yields the leading correction, proportional to $0.5Z^2$, to TF-energy, which scales as $(7/3)Z^{7/3}$. All energies quoted here are in atomic unit $(e^2/a_0)$ (where e is the electron charge and $a_0={\hslash }^2/\left(m{e}^2\right)$ is the Bohr radius ). He also showed [13] that a second correction, which arises from the gradient of electron density in the atom (as in QSM) is very accurately expressed as two by ninth (2/9) of the exchange correction (which scales as $Z^{5/3})$. Of course, the energy of strongly bound electrons has to be corrected in this procedure. An important step in this work is to use TF-density to compute the magnitude of gradient contribution. This is a very plausible assumption as the gradient contribution is only the first quantum correction ($\sim {\hslash }^2,\mathrm{where}\hslash$ is the reduced Planck’s constant) to TF-energy. So this contribution may be accounted as an accurate correction term. More elaborately, Englert and Schwinger [14] demonstrated that a new energy functional (hereafter referred as Englert-Schwinger (ES) functional) could be introduced, in terms of pseudo electron density and electrostatic potential. Then, the contribution of the gradient term is found expressible as (2/9) of exchange contribution plus the $0.5Z^2$ correction. Their computations of electron properties of cold isolated atom showed significant improvements over the results of TF model. In fact, their results of atomic biding energies differ from the corresponding results of Hartree-Fock method by about 1.1 percent (for $2\le Z\le 32$) and 0.1 percent for $Z>32$. The new model was also applied in EOS studies to compute (the zero-temperature) pressure of Al and Cu, in compressed volume states, and the results showed good agreement with those of detailed QSM [15].

The present paper is an extension of ES model to finite temperatures, and is motivated by their comment [14] on QSM proper: ’The resulting forth order differential equation for (electrostatic potential) V in QSM is too complicated, and there is the danger of taking too seriously the initial terms of an infinite series of quantum corrections’ (below Eq.3 in page-23 ). So, a relation between the electron density gradient and exchange contributions to free energy functional (at finite temperature) is established first in this paper. This is done without resorting to the properties of TF-density, as was done earlier by Englert and Schwinger. The proportionality factors between the two components are (2/9) and (3/9) in the zero-temperature and high-temperature limits, respectively. A temperature-dependent correction to the energy of strongly bound electons also arises. Thus, the approximation procedure yields a modified Thomas-Fermi-Dirac (mTFD) model, which incorporates the density gradient component accurately. Additionally, the exchange-correlation term is explicitly accounted in the present paper. Then, a new method to solve the mTFD model is presented, which is an extension of a recently developed algorithm for solving the TF model [16]. The algorithms make use of very accurate rational function approximations of Fermi-Dirac integrals [17].

The paper is organized as follows. The relevant formulas related to the QSM are collected in Section 2. The relation between the density-gradient and exchange components, leading to the finite-temperature mTFD model, is derived in Section 3. Thermodynamic functions and virial theorem are discussed in Section 4. The Newton-Ralphson algorithm to solve the mTFD equation is outlined in Section 5, and numerical results and their comparisons are provided in Section 6. Finally, a summary of the work is provided in Section 7. To make the paper totally self-contained, derivation of the various terms in the QSM functional are provided in sections Section 8.1 (kinetic energy term), Section 8.2 (density-gradient term), Section 8.3 (exchange term) , Section 8.4 (exchange-correlation-free-energy) of the Appendix. In addition, the stationary property of the functional, around its minimum, is proved in Section 8.5. The finite-temperature corrections to energy of strongly bound electrons is obtained in Section 8.6, and details of the algorithm is provided in Section 8.7 .

2. Quantum Statistical Model

Mermin [18] extended the Hohenberg-Kohn theorem by showing that, for given temperature $\left(T\right)$ and chemical potential $\left(\mu \right)$, there exists a unique functional $G[n,T]$ of electron density $n\left(\overrightarrow{r}\right)$ such that $\Omega [n,T]=G[n,T]+\int U\left(\overrightarrow{r}\right)n\left(\overrightarrow{r}\right)d\overrightarrow{r}-\mu \int n\left(\overrightarrow{r}\right)d\overrightarrow{r}$ is a minimum, where U is an external potential energy per electron. At the minimum, $\Omega$ equals the thermodynamic grand potential and $n\left(\overrightarrow{r}\right)$ is the equilibrium density in presence of U. The explicit form of $G[n,T]$ is not provided by the theorem. The QSM starts with an assumed form for $F[n,T]=G[n,T]+\int U\left(\overrightarrow{r}\right)n\left(\overrightarrow{r}\right)d\overrightarrow{r}$ in the atomic cell. This functional of n, for the specific Coulomb case $U=-Z{e}^2/r$, is expressed as [9,10]

$$F[n,T]=ʃd\overrightarrow{r}\{f_{kin}[n,T]+f_{grad}[n,T]+f_{xc}[n,T]\}-ʃd\overrightarrow{r}{\displaystyle \frac{Z{e}^2}{r}}n\left(\overrightarrow{r}\right)+{\displaystyle \frac{e^2}{2}}ʃd\overrightarrow{r}ʃd{\overrightarrow{r}}^{\prime }{\displaystyle \frac{n\left(\overrightarrow{r}\right)n\left({\overrightarrow{r}}^{\prime }\right)}{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}}.$$

(1)

The corresponding $\Omega$ has a minimum is proved in the Appendix (see Section 8.5 ). In the above expression, the volume integrals are over the atomic cell of radius R. The last two terms represent electron potential energy from the nuclear charge $Ze$ at the center of the cell, and self repulsive energy of electrons, respectively. Note that e is the magnitude of electron charge and so that $-Z{e}^2/r$ is the attractive component of electrostatic energy of one electron at r. The first three terms represent contributions from kinetic energy, density-gradient and exchange-correlation effects of electron states. Although well known, derivation of explicit expressions of these terms, from basics, are provided in different sections of the Appendix. However, for the sake of completeness, these are also briefly discussed below. Thus, the kinetic energy density is

$$f_{kin}=C_{kin}{\beta }^{-5/2}\{-(2/3)I_{3/2}\left({\eta }^{\ast }\right)+{\eta }^{\ast }{I}_{1/2}\left({\eta }^{\ast }\right)\},C_{kin}=(\sqrt{2}/{\pi }^2){[m/{\hslash }^2]}^{3/2}.$$

Here, and in what follows, $\beta ={\left(k_{B}T\right)}^{-1}$ and $I_k\left(x\right)$ is the Fermi-Dirac function of order k (see below). Further, ${\eta }^{\ast }\left(\overrightarrow{r}\right)$ is the energy parameter of the model, m is electron mass, $k_B$ is Boltzman’s constant and ℏ is the reduced Planck’s constant . The energy parameter ${\eta }^{\ast }$ is introduced in electron density as

$$n\left(\overrightarrow{r}\right)=C_n{\beta }^{-3/2}{I}_{1/2}\left[{\eta }^{\ast }\left(\overrightarrow{r}\right)\right],C_n=C_{kin}=(\sqrt{2}/{\pi }^2){[m/{\hslash }^2]}^{3/2},$$

(2)

This follows from the phase-space density $n(\overrightarrow{r},\overrightarrow{p})=[2/{(2\pi \hslash )}^3]{\{exp[\beta p^2/2m-{\eta }^{\ast }\left(\overrightarrow{r}\right)]+1\}}^{-1}$. Next, the gradient-correction to kinetic energy is

$$f_{grad}=\sigma \left(n\right){\displaystyle \frac{{(\nabla n)}^2}{n}},\sigma \left(n\right)={\displaystyle \frac{1}{24}}\left[I_{1/2}\left({\eta }^{\ast }\right)I_{1/2}^{\prime \prime }\left({\eta }^{\ast }\right){\left\{I_{1/2}^{\prime }\left({\eta }^{\ast }\right)\right\}}^{-2}\right][{\hslash }^2/m].$$

(3)

The symbol ′ in the expression for $\sigma \left(n\right)$ denotes derivative with respect to the argument. This complicated formula for the coefficient $\sigma \left(n\right)$ is derived in the Appendix. Using the asymptotic forms $I_{1/2}\sim (2/3){\eta }^{3/2}$ and $I_{1/2}\sim (\sqrt{\pi }/2)exp\left(\eta \right)$, it is easily verified that $\sigma (T\to 0)=(1/72)[{\hslash }^2/m]$ and $\sigma (T\to \infty )=(1/24)[{\hslash }^2/m]$ in the appropriate limits.

The exchange-correlation free energy density $f_{xc}$ arises from the requirement of anti-symmetry of the wave function, which includes the effects of Coulomb repulsion of electrons. Computation of this effect requires solution to an interacting many-particle quantum system. Just as in the case of $f_{kin}$, the expressions obtained within the uniform electron gas (UEG) model are employed here also. Results of extensive quantum Monte Carlo simulations are fitted in terms of electron density and temperature so as to obtain analytical formulas for $f_{xc}$ [19]. Alternately, a mapping of the quantum electron fluid model to the classical fluid model, based on Ornstein-Zernike equations, also provides good results and fits for $f_{xc}$ [20]. All these important points are discussed in the Appendix. However, note that Monte Carlo representations of $f_{xc}$ is of the form

$$f_{xc}=-(n/r_s)[G_1+G_2{\sqrt{r}}_s+G_3{r}_s]{[1+G_4{\sqrt{r}}_s+G_5{r}_s]}^{-1}.$$

(4)

Here, $r_s$ is the average inter-electron distance, defined by $(4\pi /3)r_s^3=n^{-1}$. The parameters $G_1,G_2,G_3,G_{4}and{G}_5$ are complicated (fitted) functions of T and n (see Appendix for explicit expressions of the parameters ).

Finally, note that the pure exchange term is given by

$$f_{ex}=-C_{ex}{\beta }^{-2}\underset{-\infty }{\overset{{\eta }^{\ast }}{ʃ}}{\left[I_{1/2}^{\prime }\left(u\right)\right]}^{2}du,C_{ex}=(2/{\pi }^3)e^2{[m/{\hslash }^2]}^2,$$

(5)

which is also derived in the Appendix. This is a part of $f_{xc}$, but it is based only on free-electron wave functions. However, for any finite temperature T, electrons become fully degenerate as $r_s\to 0$, and hence $f_{xc}\to f_{ex}$ in that limit.

The Fermi-Dirac functions, occurring in different expressions, are defined as

$$I_k\left(x\right)=\underset{0}{\overset{\infty }{ʃ}}dy{y}^k{[exp(y-x)+1]}^{-1}=\underset{-x}{\overset{\infty }{ʃ}}dy{(y+x)}^k{[exp\left(y\right)+1]}^{-1}.$$

These satisfy the recursion relation $I_k^{\prime }\left(x\right)=k{I}_{k-1}\left(x\right)$ (where the symbol ′ denotes derivative). This formula readily follows by differentiating the second representation given above.

Note that the energy parameter in the free gas model is ${\eta }_{\mathrm{f}}=\beta \mu$ where as it is ${\eta }_{\mathrm{tf}}=\beta [eV\left(\overrightarrow{r}\right)+\mu ]$ in the TF model. The chemical potential is $\mu$ ( which is also known as Fermi-energy at zero-temperature) and $-eV$ the electrostatic potential energy per electron. The chemical potential $\mu$ in each of the models is determined via the normalization condition $Z=d\overrightarrow{r}n\left(\overrightarrow{r}\right)$ with appropriate energy parameter $\eta$ in the density $n\left(\overrightarrow{r}\right)$. In QSM or TFD models, $\eta \left(\overrightarrow{r}\right)$ is to be determined via minimizing the functional $F[n,T]$ over the local density $n\left(\overrightarrow{r}\right)$ (Mermin’s theorem quoted above). The minimization is to be done under the constraint $Z=d\overrightarrow{r}n\left(\overrightarrow{r}\right)$, and $\mu$ enters the models as Lagrange multiplier. In all cases, $\eta \left(\overrightarrow{r}\right)$ and $n\left(\overrightarrow{r}\right)$ are related via the Fermi-Dirac function $I_{1/2}\left(\eta \right)$ as given by $n\left(\overrightarrow{r}\right)=C_n{\beta }^{-3/2}{I}_{1/2}\left[\eta \left(\overrightarrow{r}\right)\right]$, and the phase space density is $n(\overrightarrow{r},\overrightarrow{p})=[2/{(2\pi \hslash )}^3]{\{exp[\beta p^2/2m-\eta \left(\overrightarrow{r}\right)]+1\}}^{-1}$.

The condition that the functional derivative of $F[n,T]-\mu n\left(\overrightarrow{r}\right)d\overrightarrow{r}$ with respect to $n\left(\overrightarrow{r}\right)$ is zero ($\mu$ is the Langrange multiplier) yields the Euler-Lagrange equation, to determin $n\left(\overrightarrow{r}\right)$ in QSM, as

$${\displaystyle \frac{\delta }{\delta n}}{f}_{kin}-{(\nabla n)}^2{\displaystyle \frac{\partial }{\partial n}}\left[{\displaystyle \frac{\sigma }{n}}\right]-2\left[{\displaystyle \frac{\sigma }{n}}\right]{\nabla }^{2}n+{\displaystyle \frac{\delta }{\delta n}}{f}_{xc}=eV+\mu .$$

(6)

Here, the electrostatic potential V is defined as

$$V={\displaystyle \frac{Ze}{r}}-eʃd{\overrightarrow{r}}^{\prime }{\displaystyle \frac{n\left({\overrightarrow{r}}^{\prime }\right)}{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}}.$$

(7)

so that $-eV$ is the potential energy per electron. Of course, V satisfies the Poisson’s equation and boundary conditions

$${\nabla }^{2}V=-4\pi (-en),{\left[rV\right]}_{r\to 0}=Ze,and,{[\nabla V]}_R=0.$$

(8)

While the first boundary condition imposes the nuclear charge, the last one mimics the effect of the surrounding medium. The atomic cell is electrically neutral. Since $f_{kin}$ is expressed as functions of ${\eta }^{\ast }$, use of the expression $(\delta /\delta n)[\cdots ]=(\delta {\eta }^{\ast }/\delta n)(\delta /\delta {\eta }^{\ast })[\cdots ]$ reduces Eq.(6) to

$${\displaystyle \frac{{\eta }^{\ast }}{\beta }}-{(\nabla n)}^2{\displaystyle \frac{\partial }{\partial n}}\left[{\displaystyle \frac{\sigma }{n}}\right]-2\left[{\displaystyle \frac{\sigma }{n}}\right]{\nabla }^{2}n+{\displaystyle \frac{\delta }{\delta n}}{f}_{xc}=eV+\mu .$$

(9)

The derivative term $(\delta f_{xc}/\delta n)$ is an effective potential energy per electron, and it needs to be computed numerically using Eq.(4). Extensive comparison of this effective potential energy, using different approximations to $f_{xc}$, has been done [21] because of its use in average atom models.

The energy parameter ${\eta }_{\mathrm{tf}}=\beta (eV+\mu )$ for the TF model follows on neglecting the exchange-correlation and gradient terms. However, the non-linear equation for ${\eta }_{\mathrm{tfd}}$, viz.

$${\eta }_{\mathrm{tfd}}-[C_{ex}/C_n]{\beta }^{1/2}{I}_{1/2}^{\prime }\left({\eta }_{\mathrm{tfd}}\right)=\beta (eV+\mu ).$$

(10)

is to be solved in the TFD model [22], as it accounts only for exchange effects. This follows from $\delta f_{ex}/\delta n=-[C_{ex}/C_n]{\beta }^{-1/2}{I}_{1/2}^{\prime }\left({\eta }^{\ast }\right)$ where the expression $(\delta n/\delta {\eta }^{\ast })=C_n{\beta }^{-3/2}{I}_{1/2}^{\prime }\left({\eta }^{\ast }\right)$ is used . This way of solving the TFD model [22], together with the Poisson’s equation for V, is much simpler than the original method due to Cowan and Askin [6]. In lieu of Eq.(2), these authors used a generalized form of TF phase-space distribution function, viz. $n(\overrightarrow{r},\overrightarrow{p})=[2/{(2\pi \hslash )}^3]{\{exp\left[\beta (p^2/2m+E_{ex}(\overrightarrow{p},\overrightarrow{r})-eV-\mu )\right]+1\}}^{-1}$, where $E_{ex}(\overrightarrow{p},\overrightarrow{r})$ is the momentum-dependent exchange energy. While the TF and TFD models leads to divergence of $n\left(\overrightarrow{r}\right)$ as $r\to 0$, solution of full equation (9) in QSM is known [9,10] to provide a finite value of $n\left(0\right)$ due to the occurrence of derivative terms in it. Use of Poisson’s equation reduces Eq.(9) to a highly non-linear fourth order differential equation.

By deriving the relation $ʃd\overrightarrow{r}{f}_{grad}\left(\overrightarrow{r}\right)=(2/9)ʃd\overrightarrow{r}{f}_{ex}\left(\overrightarrow{r}\right)+F_c(r\to 0)$, via using the TF approximation for $n\sim n_{tf}$, and hence ${\eta }^{\ast }\sim {\eta }_{tf}=\beta (eV+\mu )$, Schwinger simplified the QSM in the T $\to 0$ limit. But he found that $F_c(r\to 0)$ is divergent because of the use of TF-density $n_{tf}$ in the procedure. However, Schwinger argued that $F_c$ can be omitted together with the substitution of the correct energy of strongly bound electrons in place of the wrong TF-energy of these electrons. He had already obtained the correction term that is needed to remedy the latter defect as $0.5Z^2(e^2/a_0)$. Thus, according to Schwinger, $(2/9)d\overrightarrow{r}{f}_{ex}\left(\overrightarrow{r}\right)$ is the non-divergent part of the gradient term that should be used in the QSM functional. This development is quite important because the contribution of the gradient term to thermodynamic properties is relatively small, particularly for temperatures where electron contribution to thermal energy and pressure are significant. Furthermore, the the gradient term is only the leading quantum correction to the semi-classical (TF model) free energy functional (see also Section 8.2). Lastly, it reduces the order of the non-linear differential equation of the model from four to two. Schwinger’s approximation scheme for finite temperatures is explored next.

3. Englert-Schwinger Functional

A relation between gradient and exchange terms is obtained below following the method of Schwinger for zero-temperature [13]. The starting point in the derivation is that the Poisson’s equation for V is also expressed as ${\nabla }^2[{\eta }^{\ast }+{\eta }_{gxc}]=4\pi \beta e^{2}n$, where ${\eta }_{gxc}$ represents the parts in Eq.(9) due to gradient and exchange-correlation effects (see also third para in Section 8.5). Now, on multiplying this equation with $-(\sigma /n)g\left({\eta }^{\ast }\right)$, where $g\left({\eta }^{\ast }\right){=}_{-\infty }^{{\eta }^{\ast }}{\left[I_{1/2}^{\prime }\left(y\right)\right]}^{2}dy$, and integrating over the atomic cell leads to

$$-ʃd\overrightarrow{r}(\sigma /n)g\left({\eta }^{\ast }\right){\nabla }^2[{\eta }^{\ast }+{\eta }_{gxc}]=-4\pi e^2\beta ʃd\overrightarrow{r}\sigma \left({\eta }^{\ast }\right)g\left({\eta }^{\ast }\right)=C_{\eta }{\beta }^3ʃd\overrightarrow{r}\sigma \left({\eta }^{\ast }\right)f_{ex}\left({\eta }^{\ast }\right)$$

where $C_{\eta }=\left(4\pi e^2\right)/C_{ex}$. Next, the term in the integral on the left is written as

$$\begin{array}{ccc}\hfill (\sigma /n)g\left({\eta }^{\ast }\right){\nabla }^2[{\eta }^{\ast }+{\eta }_{gxc}]& =& \nabla ·\{(\sigma /n)g\left({\eta }^{\ast }\right)\nabla [{\eta }^{\ast }+{\eta }_{gxc}]\}\hfill \\ & -& [{(\sigma /n)}^{\prime }g\left({\eta }^{\ast }\right)+(\sigma /n)g^{\prime }\left({\eta }^{\ast }\right)][{(\nabla {\eta }^{\ast })}^2+(\nabla {\eta }^{\ast })·(\nabla {\eta }_{gxc})]\hfill \end{array}$$

where the symbol ′ denotes derivative with ${\eta }^{\ast }$. Conversion of the divergence term to surface integrals yields

$$\begin{array}{ccc}\hfill \left[S\right(0)-S(R\left)\right]& +& ʃd\overrightarrow{r}[{(\sigma /n)}^{\prime }g\left({\eta }^{\ast }\right)+(\sigma /n)g^{\prime }\left({\eta }^{\ast }\right)][{(\nabla {\eta }^{\ast })}^2+(\nabla {\eta }^{\ast })·(\nabla {\eta }_{gxc})]\hfill \\ & =& C_{\eta }{\beta }^3ʃd\overrightarrow{r}\sigma \left({\eta }^{\ast }\right)f_{ex}\left({\eta }^{\ast }\right).\hfill \end{array}$$

(11)

The surface term $S\left(\overrightarrow{r}\right)=\{4\pi r^2(\sigma /n)g\left({\eta }^{\ast }\right)\nabla [{\eta }^{\ast }+{\eta }_{gxc}]$ is zero at R as $\nabla [{\eta }^{\ast }+{\eta }_{gxc}]\}\sim \nabla V$ vanishes there. But it tends to a finite limit at the origin as $\nabla [{\eta }^{\ast }+{\eta }_{gxc}]=\beta e\nabla V\sim -\beta Z{e}^2/r^2$. Note that n, and hence ${\eta }^{\ast }$, are finite at the origin in QSM. But, this term at $r=0$ is to be ignored together with the replacement of free energy of strongly bound electrons from a quantum mechanical model (see Eq.(13) given below). Now, noting that $\nabla n=C_n{\beta }^{-3/2}{I}_{1/2}^{\prime }\nabla {\eta }^{\ast }$, the term proportional to $(\sigma /n)$ is expressed as $(\sigma /n)g^{\prime }\left({\eta }^{\ast }\right){(\nabla {\eta }^{\ast })}^2=\left[C_n^{-2}{\beta }^3\right]f_{grad}\left({\eta }^{\ast }\right)$. Thus, also using the result $C_{\eta }{C}_n^2=4[m/{\hslash }^2]$, the relation in Eq.(11) reduces to

$$ʃd\overrightarrow{r}(1+\psi )(1+\zeta )f_{grad}\left({\eta }^{\ast }\right)=4[m/{\hslash }^2]ʃd\overrightarrow{r}\sigma \left({\eta }^{\ast }\right)f_{ex}\left({\eta }^{\ast }\right).$$

(12)

Here, $\psi =\left[{(\sigma /n)}^{\prime }g\right]/\left[(\sigma /n)g^{\prime }\right]$ and $\zeta =(\nabla {\eta }^{\ast })·(\nabla {\eta }_{gxc})/{(\nabla {\eta }^{\ast })}^2$. The expression in Eq.(12) is a generalization of Schwinger’s relation, between gradient and exchange terms at zero-temperature. The correction term $\zeta$ involving ${\eta }_{gxc}$ vanishes if TF approximation is used for n or ${\eta }^{\ast }$, and is neglected hereafter. The parameter $\psi$ tends to (-3/4) and (-1/2) in the limits $T\to 0$ and $T\to \infty$, respectively. But $\sigma \left(\eta \right)$ reduces to $(1/72)[{\hslash }^2/m]$ and $(1/24)[{\hslash }^2/m]$ in these limits. So, it is clear that the proportionality factors between $ʃd\overrightarrow{r}{f}_{grad}\left(\overrightarrow{r}\right)$ and $ʃd\overrightarrow{r}{f}_{ex}\left(\overrightarrow{r}\right)$ terms are (2/9) and (3/9) in the two limits. Thus, the Englert-Schwinger functional is obtained (in the limits of zero and high temperature) by using $f_{ex}^{\ast }$, in lieu of $f_{grad}$ with the definition $C_{ex}^{\ast }=C_s{C}_{ex}$. The constant $C_s$ takes values (2/9) and (3/9) in the limits $T\to 0$ and $T\to \infty$, respectively. Since the range of these limits is very small, all computations reported in this paper employed the geometric average value $C_s=0.2722$. The resulting free energy functional is a modified form of TFD functional (denoted as mTFD), the modification being the addition of $f_{ex}^{\ast }$ and the exchange-correlation component. The functional in mTFD is given by

$$F_{es}[n,T]=\Delta E+ʃd\overrightarrow{r}\{f_{kin}[n,T]+f_{ex}^{\ast }[n,T]+f_{xc}[n,T]\}-ʃd\overrightarrow{r}{\displaystyle \frac{Z{e}^2}{r}}n\left(\overrightarrow{r}\right)+{\displaystyle \frac{e^2}{2}}ʃd\overrightarrow{r}ʃd{\overrightarrow{r}}^{\prime }{\displaystyle \frac{n\left(\overrightarrow{r}\right)n\left({\overrightarrow{r}}^{\prime }\right)}{|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|}}.$$

(13)

Here, $\Delta E={\int }_{r_0}d\overrightarrow{r}[e_{qm}\left(\overrightarrow{r}\right)-e_{tf}\left(\overrightarrow{r}\right)]$ is the correction to energy of strongly bound electrons within a radius $r_0$. It refers only to kinetic and electrostatic energies (i.e. TF model) because the wrong contribution from gradient term is already omitted. Further, quantum Monte Carlo simulations accurately model the exchange-correlation term in the free energy functional. The quantity $\Delta E$ is taken as correction in energy as it is computed by conserving the number of states within $r_0$. As mentioned earlier, $\Delta E(T\to 0)$ equals $0.5Z^2(e^2/a_0)$ at zero-temperature when all bound states are considered. It is computed in the Appendix for finite-temperature (see Section 8.6). The scaled factor $\Delta E/(Z^2{e}^2/a_0)$ verses scaled temperature parameter $\beta Z^2(e^2/a_0)$, for the lowest 10 bound states, is shown in Figure 1. The magnitude of the hump in the figure decreases, and the asymptotic value approaches 0.5, if more levels are considered. The symbols in this figure depict the fit function $0.46\times {\{exp(-2.5\ast \beta Z^2{E}_0)+1\}}^{-1}$ ( with error $\le 4\%$) for the correction factor. The difference in the occupation probability of these levels versus $\beta Z^2(e^2/a_0)$ (see the inset figure) is less than 1 percent in the entire range of temperature.

The energy parameter ${\eta }^{\ast }$ obeys the equation ${\eta }^{\ast }\left(x\right)+\beta [\delta f_{gxc}/\delta n]=\beta [eV+\mu ]$. For the sake of compact notation, a combined fenergy density $f_{gxc}=f_{ex}^{\ast }+f_{xc}$ is introduced here. While $f_{ex}^{\ast }$ represents a part of the gradient term, the second is for exchange-correlation contribution. The first term is explicitly evaluated as $[\delta f_{ex}^{\ast }/\delta n]=-[C_{ex}^{\ast }/C_n]{\beta }^{-1/2}{I}_{1/2}^{\prime }\left({\eta }^{\ast }\right)$. Since electron density $n\left(x\right)=C_n{\beta }^{-3/2}{I}_{1/2}\left[{\eta }^{\ast }\left(x\right)\right]$, Poisson’s equation for V, written in terms of the TF function $\varphi \left(x\right)$, is

$${\displaystyle \frac{d^2}{d{x}^2}}\varphi =\mathcal{A}x{I}_{1/2}\left[{\eta }^{\ast }\left(x\right)\right].$$

(14)

where $\varphi \left(x\right)/x=\beta \left[eV\right(r)+\mu ]$. The independent variable is $x=r/R$ and $\mathcal{A}=\left[4\pi R^2{e}^2\right]C_n{\beta }^{-1/2}$. Boundary conditions on $\varphi \left(x\right)$ are reduced to $\varphi \left(0\right)=\beta Z{e}^2/R$ and ${\varphi }^{\prime }\left(1\right)=\varphi \left(1\right)$. After solving the Eq.(14), the chemical potential is obtained as $\mu =\varphi \left(1\right)$ because of the normalization $V\left(R\right)=0$.

4. Thermodynamic Functions

Derivations of the thermodynamic functions within the TF model [3,16] are extended due to the presence of modified-exchange term (representing the density-gradient effects) and exchange-correlation term. Electron pressure is obtained using the definition $P=-(\partial F/\partial \mathcal{V})$ and Eq.(13). In this step T and ${\left[n\right]}_b$ are kept constant, where the subscript $’b’$ denotes cell-boundary. Using the relation $\partial [\cdots ]/\partial \mathcal{V}={\left(4\pi R^2\right)}^{-1}\partial [\cdots ]/\partial R$ for derivative with spherical cell-volume, it readily follows that pressure is given by $P=-{[f_{kin}+f_{gxc}]}_b+{\left[neV\right]}_b-\mu \int d\overrightarrow{r}(\partial n/\partial \mathcal{V})$. Last term arises from $-\int d\overrightarrow{r}(\partial [\cdots ]/\partial \mathcal{V})=-\int d\overrightarrow{r}(\partial [\cdots ]/\partial n)(\partial n/\partial \mathcal{V})=-\int d\overrightarrow{r}\mu \partial n/\partial \mathcal{V})$. Also, note the result $\partial Z/\partial \mathcal{V}=0={\left[n\right]}_b+\int d\overrightarrow{r}(\partial n/\partial \mathcal{V})$. Substitution of free energy densities shows that the terms $\mu {\left[n\right]}_b$ and ${\left[neV\right]}_b$ cancel out, so that P is expressed as

$$P={\displaystyle \frac{2}{3}}{C}_n{\beta }^{-5/2}{I}_{3/2}\left[{\eta }^{\ast }\left(1\right)\right]+{[n\delta f_{gxc}/\delta n-f_{gxc}]}_b.$$

(15)

Note that ${[n\delta f_{ex}^{\ast }/\delta n-f_{ex}^{\ast }]}_b$ is evaluated as $C_{ex}^{\ast }{\beta }^{-2}[{\int }_{-\infty }^{{\eta }_1^{\ast }}{\left[I_{1/2}^{\prime }\left(y\right)\right]}^{2}dy-I_{1/2}^{\prime }\left[{\eta }_1^{\ast }\right]I_{1/2}\left[{\eta }_1^{\ast }\right]]$, where ${\eta }_1^{\ast }={\eta }^{\ast }\left(1\right)$. Accurate approximations exists [23] for the integral term involving ${\left[I_{1/2}^{\prime }\left(y\right)\right]}^2$.

The potential energy (per cell) of electrons is $-e\int d\overrightarrow{r}[Ze/r+(1/2)V_e]n\left(\overrightarrow{r}\right)$ where $V_e$ is the electrostatic potential due to electrons. This is the same as $-(e/2)\int d\overrightarrow{r}[Ze/r+V]n\left(\overrightarrow{r}\right)$, so potential energy is given by

$$E_{pot}=3\mathcal{V}{C}_n{\beta }^{-5/2}(-1/2)\underset{0}{\overset{1}{ʃ}}dx{x}^2[\varphi \left(0\right)/x-\beta \mu +\varphi \left(x\right)/x]I_{1/2}\left({\eta }^{\ast }\right)$$

On integrating the TF equation [see Eq.(14)], it is found that ${\int }_0^{1}dxx{I}_{1/2}$ equals $[\varphi \left(0\right)/\mathcal{A}][{\varphi }^{\prime }\left(1\right)-{\varphi }^{\prime }\left(0\right)]$. Similarly, multiplying the TF equation with x and integrating shows that the second integral equals $[-\beta \mu /\mathcal{A}]\varphi \left(0\right)$. Thus, $E_{pot}$ is obtained as

$$E_{pot}=3\mathcal{V}{C}_n{\beta }^{-5/2}(1/2)\{[\varphi \left(0\right){\varphi }^{\prime }\left(0\right)/\mathcal{A}]-\underset{0}{\overset{1}{ʃ}}dxx\varphi \left(x\right)I_{1/2}\left({\eta }^{\ast }\right)\}$$

(16)

The last integral needs to be computed after obtaining the solution $\varphi \left(x\right)$ and ${\eta }^{\ast }\left(x\right)$

Kinetic energy of electrons in the cell $\int d\overrightarrow{r}\int d\overrightarrow{p}(p^2/2m)n(\overrightarrow{r},\overrightarrow{p})$ is expressed in the integral

$$E_{kin}=3\mathcal{V}{C}_n{\beta }^{-5/2}\underset{0}{\overset{1}{ʃ}}dx{x}^2{I}_{3/2}\left[{\eta }^{\ast }\left(x\right)\right].$$

The phase space density $n(\overrightarrow{r},\overrightarrow{p})=[2/{(2\pi \hslash )}^3]{\{exp[\beta p^2/2m-{\eta }^{\ast }\left(\overrightarrow{r}\right)]+1\}}^{-1}$ provides this expression. One integration by parts, and the equation for ${\eta }^{\ast }$ yields

$$\begin{array}{ccc}\hfill E_{kin}& =& 3\mathcal{V}{C}_n{\beta }^{-5/2}\{(1/3)I_{3/2}\left[{\eta }^{\ast }\left(1\right)\right]-(1/2){\int }_0^{1}dx{x}^3{I}_{1/2}\left({\eta }^{\ast }\right)[(x{\varphi }^{\prime }-\varphi )/x^2]\}+\hfill \\ & +& (3/2)\mathcal{V}{\int }_0^{1}dx{x}^{3}n(d/dx)[\delta f_{gxc}/\delta n].\hfill \end{array}$$

The TF equation [see Eq.(14)] shows that ${\int }_0^{1}dxx{I}_{1/2}\left({\eta }^{\ast }\right)(x{\varphi }^{\prime }-\varphi )=(1/2)[\varphi \left(0\right){\varphi }^{\prime }\left(0\right)/\mathcal{A}-{\int }_0^{1}dxx\varphi \left(x\right)I_{1/2}\left({\eta }^{\ast }\right)dx]$. Then, $E_{kin}$ is obtained as

$$\begin{array}{ccc}\hfill E_{kin}& =& 3\mathcal{V}{C}_n{\beta }^{-5/2}\{(1/3)I_{3/2}\left[{\eta }^{\ast }\left(1\right)\right]-(1/4)[\varphi \left(0\right){\varphi }^{\prime }\left(0\right)/\mathcal{A}-{\int }_0^{1}dxx\varphi \left(x\right)I_{1/2}\left({\eta }^{\ast }\right)dx]\}+\hfill \\ & +& (3/2)\mathcal{V}{\int }_0^{1}dx{x}^{3}n(d/dx)[\delta f_{gxc}/\delta n].\hfill \end{array}$$

(17)

To use this expression, the last integral also needs to be done numerically. More importantly, Eq.(17) together with Eq.(16) show that

$$2E_{kin}+E_{pot}=2\mathcal{V}{C}_n{\beta }^{-5/2}{I}_{3/2}\left[{\eta }^{\ast }\left(1\right)\right]+3\mathcal{V}{\int }_0^{1}dx{x}^{3}n(d/dx)[\delta f_{gxc}/\delta n].$$

(18)

With out the last integral term (see also Eq.(15)), this equation provides the virial theorem within TF model. Now, one integration by parts provides the result

$$3\int d\overrightarrow{r}[n\delta f_{gxc}/\delta n-f_{gxc}]=3\mathcal{V}\{[n\delta f_{gxc}/\delta n-f_{gxc}){]}_b-{\int }_0^{1}dx{x}^{3}n(d/dx)[\delta f_{gxc}/\delta n]\}$$

(19)

The result $(d/dx)[n\delta f_{gxc}/\delta n-f_{gxc}]=n(d/dx)[\delta f_{gxc}/\delta n]$ is used here. The three expressions in Eq.(15), Eq.(18) and Eq.(19) yield the virial theorem in mTFD model, viz.

$$2E_{kin}+E_{pot}+3\int d\overrightarrow{r}[n\partial f_{gxc}/\partial n-f_{gxc}]=3\mathcal{V}P$$

(20)

The same relation is found [10] in the original QSM because Eq.(18) is also valid there. In the zero-temperature limit, $(n\partial f_{gxc}/\partial n-f_{gxc})\sim n^{4/3}$. Then the virial theorem is reduced to $2E_{kin}+E_{pot}+E_{gxc}=3\mathcal{V}P$. Feynman et al has argued that the last relation is valid in general [2]. It is noted that Cowan et al have used the same for finite temperatures as well [6].

It is necessary to obtain total energy of electrons using the relation $E={[\partial \left(\beta F\right)/\partial \beta ]}_{\mathcal{V},n}$. The density profile $n\left(\overrightarrow{r}\right)$ is kept constant so as to avoid any contribution due to its variation. Multiplication of Eq.(13) by $\beta$ and differentiation with $\beta$ yields

$$E=\Delta E+E_{kin}-\int d\overrightarrow{r}n{\eta }_{\beta }^{\ast }+\int d\overrightarrow{r}{\left[n{\eta }^{\ast }\right]}_{\beta }+[\partial \left(\beta F_{gxc}\right)/\partial \beta ]+E_{pot}+\beta {\left[E_{pot}\right]}_{\beta }$$

Here, the subscript $\beta$ stands for partial derivative with $\beta$ at constant $\mathcal{V}$ and n. The second and third terms cancel and the last term vanishes. Thus, E reduces to

$$E=\Delta E+E_{kin}+E_{pot}+[\partial \left(\beta F_{xc}\right)/\partial \beta ]+[\partial \left(\beta F_{ex}^{\ast }\right)/\partial \beta ]$$

Substitution of the explicit form for $F_{ex}^{\ast }$ yields

$$E=\Delta E+E_{kin}+E_{pot}+E_{xc}+C_{ex}^{\ast }{\beta }^{-2}\int d\overrightarrow{r}[(-3/2)I_{1/2}^{\prime }{I}_{1/2}+{\int }_{-\infty }^{{\eta }^{\ast }}du{\left(I_{1/2}^{\prime }\right)}^2]$$

(21)

The result $n_{\beta }=0=C_n[-(3/2){\beta }^{-5/2}{I}_{1/2}+{\beta }^{-3/2}{I}_{1/2}^{\prime }{\eta }_{\beta }^{\ast }]$ is used in the last step. The integral terms reduces to $E_{ex}^{\ast }$ at zero-temperature. A fitted formula for $E_{xc}$, obtained from liquid state theory [20], is given in the Appendix (see Section 8.4). A fitted function for the correction term for energy of lowest 10 strongly bound states, $\Delta E=0.46Z^2(e^2/a_0)\times {\{exp(-2.5\ast \beta Z^2{E}_0)+1\}}^{-1}$, is found adequate even for all temperatures (see Figure 1).

Using the phase-space density function $n(\overrightarrow{r},\overrightarrow{p})$, the number of free electrons (with positive energy) is expressed as

$$Z^{\ast }=C_n{\beta }^{-3/2}3\mathcal{V}{\int }_0^{1}dx{x}^2{I}_{1/2}^{\ast }[{\eta }^{\ast }\left(x\right),Y\left(x\right)].$$

(22)

Here, $I_{1/2}^{\ast }(x,y)={\int }_y^{\infty }dz\sqrt{z}{[exp\{z-x\}+1]}^{-1},(y\ge 0)$ is the incomplete Fermi-Dirac integral [24] of order one-half. The positive quantity $Y\left(x\right)=({\eta }^{\ast }\left(x\right)-\beta \mu )$ arises from the condition [22] of positive total energy of an electron, viz. $\beta p^2/\left(2m\right)>({\eta }^{\ast }-\beta \mu )$. Note that this energy includes the contributions from gradient and exchange-correlation terms.

5. Numerical Scheme

The numerical scheme needed to solve Eq.(14) is easily obtained via a slight modification of the recently proposed algorithm [16] for the TF model. First of all, it is rewritten as an integral equation

$$\varphi \left(x\right)+\mathcal{A}\underset{0}{\overset{1}{ʃ}}\mathcal{K}(x,y)y{I}_{1/2}\left[{\eta }^{\ast }\left(y\right)\right]dy-\varphi \left(1\right)x=\varphi \left(0\right).$$

(23)

The boundary conditions are incorporated above and the kernel is $\mathcal{K}(x,y)=x$ for $y\ge x$ and $\mathcal{K}(x,y)=y$ for $y\le x$. Further, ${\eta }^{\ast }\left(x\right)$ obeys the non-linear equation

$${\eta }^{\ast }\left(x\right)-\mathcal{B}{I}_{-1/2}\left({\eta }^{\ast }\right)+\beta U_{xc}=\varphi \left(x\right)/x,\mathcal{B}=(1/2)[C_{ex}^{\ast }/C_n]{\beta }^{1/2}$$

(24)

Here, $U_{xc}=[\delta f_{xc}/\delta n]$ denotes the effective potential energy (per electron) due to exchange-correlation term. As already pointed out, $U_{xc}$ is to be computed via a finite difference formula. The second term above is $U_{ex}^{\ast }=[\delta f_{ex}^{\ast }/\delta n]$. In comparison to the algorithm for TF model [16], the modification required in the present case is just finding the solution of Eq.(24). But this is easily done via Newton’s method.

To use Newton’s method for Eq.(23), guess-solutions $\tilde{\varphi }$ and ${\tilde{\eta }}^{\ast }$ are assumed first. Then, from Eq.(24), it is noted that $[\partial {\eta }^{\ast }/\partial \varphi ]=x^{-1}{[1-\mathcal{B}{I}_{-1/2}^{\prime }\left({\tilde{\eta }}^{\ast }\right)+\beta U_{xc}^{\prime }]}^{-1}$. Just like in $I_{-1/2}^{\prime }$, the symbol ′ on $V_{xc}^{\prime }$ indicates derivative with ${\eta }^{\ast }$. The non-linear term in Eq.(23) is approximated via the expansion $I_{1/2}\left[{\eta }^{\ast }\right]\sim I_{1/2}\left[{\tilde{\eta }}^{\ast }\right]+(1/2)I_{-1/2}\left[{\tilde{\eta }}^{\ast }\right][\partial {\eta }^{\ast }/\partial \varphi ](\varphi -\tilde{\varphi })$. This follows from Taylor’s expansion of $I_{1/2}\left[{\eta }^{\ast }\right]$ (using the relation $I_n^{\prime }=n{I}_{n-1}$) around the function ${\tilde{\eta }}^{\ast }$. Substitution into Eq.(23) yields the linear integral equation

$$\begin{array}{ccc}\hfill \varphi \left(x\right)& +& \mathcal{A}\underset{0}{\overset{1}{ʃ}}\mathcal{K}(x,y)Q\left({\tilde{\eta }}^{\ast }\right)\varphi \left(y\right)dy-\varphi \left(1\right)x=S\left(x\right),\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill Q({\tilde{\eta }}^{\ast }& )& ={[1-\mathcal{B}{I}_{-1/2}^{\prime }\left({\tilde{\eta }}^{\ast }\right)+\beta U_{xc}^{\prime }\left({\tilde{\eta }}^{\ast }\right)]}^{-1}(1/2)I_{-1/2}\left[{\tilde{\eta }}^{\ast }\right].\hfill \end{array}$$

(26)

The known function $S\left(x\right)$ is given by

$$S\left(x\right)=\varphi \left(0\right)+{\mathcal{A}}_0^1\mathcal{K}(x,y)\{Q\left({\tilde{\eta }}^{\ast }\right)\tilde{\varphi }\left(y\right)-y{I}_{1/2}\left({\tilde{\eta }}^{\ast }\right)\}dy,$$

(27)

If newton’s method converges, the terms involving Q cancel out and the original Eq.(23) for $\varphi \left(x\right)$ is recovered. Therefore, to simplify algorithm, $U_{xc}^{\prime }$ is approximated as $U_{ex}^{\prime }$ in the in the Q-function. This amounts to modifying the definition of $\mathcal{B}$ in Eq.(24) with the constant $C_{ex}^{†}=(C_s+1)C_{ex}$ in place of $C_{ex}^{\ast }$. Once the linear equation is solved to obtain the improved solution $\varphi$, it is necessary to update ${\eta }^{\ast }$ also. In the numerical computations, it is found that this updating is best achieved via Picad’s iterations $(\sim 5)$ of Eq.(24). The linear approximation connecting ${\eta }^{\ast }$ and ${\tilde{\eta }}^{\ast }$ is found to slowdown convergence for expanded volume states and very low temperatures.

It is evident from Eq.(24) that ${\tilde{\eta }}^{\ast }\sim {\eta }^{\ast }\to [\varphi \left(0\right)/x]$ as $x\to 0$. However, the integral and the source term (see 25 and 27) in the linear Fredholm integral equation (FIE) exist for $x\ge 0$ because $I_{-1/2}[\varphi /y]\sim {[\varphi \left(0\right)/y]}^{1/2}$ and $I_{1/2}[\varphi /y]\sim (2/3){[\varphi \left(0\right)/y]}^{3/2}$ as $y\to 0$. So, starting with guess solution $\tilde{\varphi }$, and corresponding ${\tilde{\eta }}^{\ast }$, the FIE is solvable employing well known approaches, say, the Nyström’s method [25]. Thereafter, ${\eta }^{\ast }\left(x\right)$ is updated via Picad’s iterations. A few repetitions of this whole procedure, together with the replacements ${\tilde{\eta }}^{\ast }\leftarrow {\eta }^{\ast }$ and $\tilde{\varphi }\leftarrow \varphi$ provide the required solution, if Newton’s method converges. But for the occurrence of Q-function [see Eq.(26)] in lieu of $(1/2)I_{-1/2}\left({\tilde{\eta }}^{\ast }\right)$ and Picad’s iterations, the algorithm is identical to that for the TF model [16]. Therefore, additional details of the algorithm are provided in Section 8.7 of the Appendix.

After obtaining the converged profiles of $\varphi \left(x\right)$ and ${\eta }^{\ast }\left(x\right)$, the integral $I_1={\int }_{-\infty }^{{\eta }^{\ast }\left(1\right)}{\left[I_{1/2}^{\prime }\left(y\right)\right]}^{2}dy$ is needed (see Eq.(15)) for pressure. An accurate numerical fit [23] is available for $I_1$. For various components of energy, the integrals $I_2{=}_0^{1}dxx\varphi \left(x\right)I_{1/2}\left({\eta }^{\ast }\right)$, $I_3={\int }_0^{1}dx{x}^3{I}_{1/2}\left({\eta }^{\ast }\right)(d/dx)[\delta f_{xc}/\delta n]$ and $I_4={\int }_0^{1}dx{x}^2{I}_{1/2}^{\prime }\left({\eta }^{\ast }\right)I_{1/2}\left({\eta }^{\ast }\right)]$ are required (see equations (16) ), (17) and (21). These integrals are evaluated by dividing the interval [0,1] to sub-intervals and using piece-wise quadratic interpolation functions, as done for $\varphi \left(x\right)$ (see Section 8.7). The same technique is used for $[\delta f_{xc}/\delta n]$ and $I_{1/2}^{\ast }[{\eta }^{\ast }\left(x\right),Y\left(x\right)]$ to facilitate the computation of $I_3$ and ionization $Z^{\ast }$.

6. Numerical Results

The method described above is now applied to compute thermodynamic properties of Cu and Al. Results of electron pressure for Cu, obtained in this paper, are given in Figure 2. The symbols in this figure denote the results of QSM reported by Perrot [10]. The curves marked 1, 2 and 3 correspond to $\rho /{\rho }_0$ values 10, 1 and 0.1, respectively, where ${\rho }_0=8.93{\mathrm{g}/\mathrm{cm}}^3$. The agreement is good at much less computational effort (in comparison to QSM proper), thereby demonstrating that mTFD is an excellent approximation to the detailed QSM. Similar results for electron energy are given in Figure 3 for three values (10, 1 and 0.1) of density ratio $\rho /{\rho }_0$. The binding energy of Cu at $\rho /{\rho }_0=1$ is subtracted from the computed values of total energy, as done in the results of Perrot [10]. These comparisons show that the main parts of corrections to TF model, due to exchange-correlation and gradient effects are captured in the mTDF model. It extends the arguments of Schwinger [13], that the gradient effects can be accurately treated as an exchange term, to finite temperature cases. Next, the different components of energy (kinetic, gradient-exchange part, exchange-correlation and potential) for Cu for the case $\rho /{\rho }_0=1$ are shown in Figure 4.

Similar set of results are obtained for Al and are shown in Figure 5 for pressure and Figure 6 for energy. The result for energy corresponding to QSM for the case of 50 eV temperature is from Fromy et al [11], which agree better with those of mTFD. Agreement with results of QSM is again good. Finally, Figure 7 shows the variation of ionization (${\mathrm{Z}}^{\ast }$) versus temperature for density rations 10, 2, 1 and 0.1. Note that ${\mathrm{Z}}^{\ast }$ at 300 K are 5.06 and 0.966 for $\rho /{\rho }_0=10$ and 2, respectively. An alternate definition ${\mathrm{Z}}^{\ast }$$=\mathcal{V}n\left(R\right)$ is also used quite often [21]. The values based on this relation and those using Eq.(22 are compared (for $\rho /{\rho }_0=1$ ) in the inset figure. All computations leading to the results given above used 20 sub-intervals in [0,1]. Newton’s method converged in less than 10 iterations (criteria $|\varphi -\tilde{\varphi }|<{10}^{-10}$) for all temperatures in the compressed volume states. For expanded volume states and low temperatures ($\sim 0.01$ eV) around 20 iterations were needed.

7. Summary

A simplified version of the QSM for computing the equation of state of electrons in materials is presented. To this end, Englert-Schwinger approximation scheme (at zero-temperature) of QSM is extended to finite temperatures. This procedure leads to a modified Thomas-Fermi-Dirac model (mTFD); the modification being the addition of a term proportional to the exchange free energy in lieu of the gradient term. Of course, the usual exchange-correlation effects are accounted accurately. An efficient numerical algorithm based on Newton’s method to solve the mTFD equation is developed in this paper. The basic idea is to rewrite the mTFD equation as a non-linear Fredholm integral equation and use Newton’s method. Explicit formulas for thermodynamic functions are derived which facilitates computation of pressure and various components of energy of electrons. Numerical results obtained for Cu and Al, using mTFD model, are compared with those based on detailed QSM. Good agreement is found for compressed and expanded volume states. Derivation of explicit expressions of different terms in the free energy functional of QSM are provided in the Appendix. Additionally, the stationary property of QSM is established and temperature-dependent corrections to energy of strongly bound electrons is obtained.

8.1. Kinetic Free Energy

The starting point is definition of electron density n in the uniform free-gas model (see Eq.(2)) given by $n=C_n{\beta }^{-3/2}{I}_{1/2}\left({\eta }_f\right)$ where $C_n=C_{kin}=(\sqrt{2}/{\pi }^2){[m/{\hslash }^2]}^{3/2}$. The thermodynamic definition of chemical potential is $\partial f_{kin}/\partial n=\mu ={\beta }^{-1}{\eta }_f$. It should be noted that temperature (i.e. $\beta$) and volume V are kept unaltered here. Use of the expression $(\partial /\partial n)[\cdots ]={\left[C_n{\beta }^{-3/2}{I}_{1/2}^{\prime }\left({\eta }_f\right)\right]}^{-1}(\partial /\partial {\eta }_f)[\cdots ]$ gives $f_{kin}=C_n{\beta }^{-5/2}\int {\eta }_f{I}_{1/2}^{\prime }\left({\eta }_f\right)d{\eta }_f$ . This is re-written as $f_{kin}=C_n{\beta }^{-5/2}\int {\eta }_{f}d\left[I_{1/2}\left({\eta }_f\right)\right]$. Integration by parts yields $f_{kin}=C_n{\beta }^{-5/2}\{{\eta }_f{I}_{1/2}\left({\eta }_f\right)-\int I_{1/2}\left({\eta }_f\right)d{\eta }_f\}$. On substituting $I_{3/2}^{\prime }\left({\eta }_f\right)=(3/2)I_{1/2}\left({\eta }_f\right)$, integration readily yields $f_{kin}=C_n{\beta }^{-5/2}\{{\eta }_f{I}_{1/2}\left({\eta }_f\right)-(2/3)I_{3/2}\left({\eta }_f\right)\}$. The same expression with general ${\eta }^{\ast }$ is used in Eq.(1).

An alternate method is to use the Gibbs-Helmholtz equation $-T^2\partial (F/T)/\partial T=\partial \left(\beta F\right)/\partial \beta =E$. The expression for total energy is $E=\int d\overrightarrow{r}[C_n{\beta }^{-5/2}{I}_{3/2}\left({\eta }_{tf}\right)-eUn\left(\overrightarrow{r}\right)]$, where the spatially varying energy parameter is ${\eta }_{tf}=\beta (eU+\mu )$ with U denoting an external potential. Integration over $\beta$ yields $\beta F=\int d\overrightarrow{r}\int [C_n{\beta }^{-5/2}{I}_{3/2}\left({\eta }_{tf}\right)-eUn\left(\overrightarrow{r}\right)]d\beta$. Integration by parts (over $\beta$) and use of the relation $I_{3/2}^{\prime }=(3/2)I_{1/2}$ yields

$$\beta F=\int d\overrightarrow{r}\{C_n[(-2/3){\beta }^{-3/2}{I}_{3/2}\left({\eta }_{tf}\right)+\int {\beta }^{-3/2}{I}_{1/2}\left({\eta }_{tf}\right)d{\eta }_{tf}]-\beta eUn\left(\overrightarrow{r}\right)+\int \beta eUd\left[n\right]\}.$$

Another integration by parts (over ${\eta }_{tf}$) shows that

$$\beta F=\int d\overrightarrow{r}\{C_n[(-2/3){\beta }^{-3/2}{I}_{3/2}\left({\eta }_{tf}\right)+{\beta }^{-3/2}{\eta }_{tf}{I}_{1/2}\left({\eta }_{tf}\right)]-\int {\eta }_{tf}d\left[n\right]-\beta eUn\left(\overrightarrow{r}\right)+\int \beta eUd\left[n\right]\}$$

where the definition of n is used. Now, the terms consisting $\int d\left[n\right]$ cancel because of the normalization $Z=\int d\overrightarrow{r}n$. These steps also show that terms corresponding to TF-model in Eq.(1) follow with the addition of electron-electron potential energy.

8.2. Gradient Energy Density

The approach outlined in this section is along the work of Brack et al [28], however, readers may find that it is more complete and self-contained. The computation of the gradient contribution, to the finite-temperature free energy functional, starts with the definition of density distribution of a single-particle quantum system

$$\rho \left(\overrightarrow{r}\right)=2\sum _{n-all}{\psi }_n^{\ast }\left(\overrightarrow{r}\right){\psi }_n\left(\overrightarrow{r}\right){\rho }_{\mu }(\beta ,{ϵ}_n).$$

(28)

where ${\rho }_{\mu }(\beta ,{ϵ}_n)={\{1+exp\left[\beta ({ϵ}_n-\mu )\right]\}}^{-1}$ is the Fermi-Dirac occupation probability with ${ϵ}_n$ denoting the energy of $n^{th}$ state. The summation in Eq.(28) is over all states (bound and free) and the factor 2 accounts for spin orientations. This definition is applicable to electrons in the atomic cell as they can be ascribed to single-particle states corresponding to a universal (density-temperature-dependent potential) $U\left(\overrightarrow{r}\right)$ (Mermin’s extension of Homberg-Kohn theorem). Inserting the $\delta$-function $\delta (ϵ-{ϵ}_n)$, Eq.(28) is rewritten as

$$\rho \left(\overrightarrow{r}\right)={\int }_{-\infty }^{\infty }dϵ\sum _{n-all}2{\psi }_n^{\ast }\left(\overrightarrow{r}\right){\psi }_n\left(\overrightarrow{r}\right)\delta (ϵ-{ϵ}_n){\rho }_{\mu }(\beta ,ϵ).$$

Next, use the Fourier representation $\delta \left(x\right)={\left(2\pi \right)}^{-1}{\int }_{-\infty }^{\infty }exp\left(itx\right)dt$ (with $i=\sqrt{-1}$) to obtain

$$\rho \left(\overrightarrow{r}\right)={\int }_{-\infty }^{\infty }dϵ{\left(2\pi \right)}^{-1}{\int }_{-\infty }^{\infty }dtC(\overrightarrow{r},it)exp\left(itϵ\right){\rho }_{\mu }(\beta ,ϵ).$$

Here, a new density $C(\overrightarrow{r},it)={\sum }_{n-all}\left[2{\psi }_n^{\ast }\left(\overrightarrow{r}\right){\psi }_n\left(\overrightarrow{r}\right)\right]exp(-it{ϵ}_n)$, called the Bloch density, is introduced. A nice feature of this expression is that temperature dependent part occurs separately. Now, on writing $[dϵexp\left(itϵ\right)]={\left(it\right)}^{-1}[dexp\left(itϵ\right)]$ and integrating by parts yields

$$\rho \left(\overrightarrow{r}\right)={\left(2\pi i\right)}^{-1}{\int }_{-\infty }^{\infty }(dt/t)C(\overrightarrow{r},it){\int }_{-\infty }^{\infty }dϵexp\left(itϵ\right)\beta {[2cosh(\beta (ϵ-\mu )/2)]}^{-2}.$$

The expression $d\left\{{[1+exp\left(x\right)]}^{-1}\right\}=-{[2cosh(x/2)]}^{-2}dx$ is employed above. Using the Fourier transform ${\int }_{-\infty }^{\infty }dxexp\left(i\omega x\right){[cosh\left(x\right)]}^{-2}=\pi \omega /sinh(\pi \omega /2)$, above equation reduces to

$$\rho \left(\overrightarrow{r}\right)={\left(2\pi i\right)}^{-1}{\int }_{-\infty }^{\infty }(dt/t)C(\overrightarrow{r},it)exp\left(it\mu \right){\displaystyle \frac{(\pi t/\beta )}{sinh(\pi t/\beta )}}.$$

(29)

Here, temperature T appears only in the factor $\left[\right(\pi t/\beta )/sinh(\pi t/\beta \left)\right]$ which approaches 1 when $T\to 0$. Next, the expression in Eq.(29) is evaluated semi-classically so that the lowest order term is the density in TF model. Higher order terms follow from a series expansion in powers of ℏ.

8.2.1. Bloch Density

To derive the semi-classical approximation, the Bloch-density is written as [29]

$$C(\overrightarrow{r},\tau )=\sum _{n-all}2{\psi }_n^{\ast }\left(\overrightarrow{r}\right)exp(-\tau {ϵ}_n){\psi }_n\left(\overrightarrow{r}\right)=\sum _{n-all}2{\psi }_n^{\ast }\left(\overrightarrow{r}\right)exp(-\tau H){\psi }_n\left(\overrightarrow{r}\right)$$

(30)

where $\tau$ is just a parameter, $H=({\hslash }^2/2m){\nabla }^2+U\left(\overrightarrow{r}\right)$ is the Hamiltonian operator, and it obeys the equation $H{\psi }_n={ϵ}_n{\psi }_n$. The last term in Eq.(30) is the trace of the matrix with elements ${\rho }_{n,m}=[2{\psi }_n^{\ast }\left(\overrightarrow{r}\right)exp(-\tau H){\psi }_m\left(\overrightarrow{r}\right)]$. Now, trace of a matrix is unaltered if its elements are computed with any complete set of orthogonal functions. So, ${\psi }_n\left(\overrightarrow{r}\right)$ is replaceable with free-particle wave function $exp(i\overrightarrow{p}·\overrightarrow{r}/\hslash )$ and the sum is changeable to integration over $\overrightarrow{p}$. This is a crucial point in the development of semi-classical approximations. With these alterations, $C(\overrightarrow{r},\tau )$ is expressed as [29]

$$C(\overrightarrow{r},\tau )=[2/{(2\pi \hslash )}^3]d\overrightarrow{p}exp(-i\overrightarrow{p}·\overrightarrow{r}/\hslash )exp(-\tau H)exp(i\overrightarrow{p}·\overrightarrow{r}/\hslash )$$

(31)

The factor in front provides the correct number of states in $\left[d\overrightarrow{r}d\overrightarrow{p}\right]$, although $d\overrightarrow{r}$ is not written explicitly. Denoting $I=exp(-i\overrightarrow{p}·\overrightarrow{r}/\hslash )exp(-\tau H)exp(i\overrightarrow{p}·\overrightarrow{r}/\hslash )$, it is easily verified that $\partial I/\partial \tau =-exp(-i\overrightarrow{p}·\overrightarrow{r}/\hslash )Hexp(i\overrightarrow{p}·\overrightarrow{r}/\hslash )I$. Further, $I\to 1$ as $\tau \to 0$. Then, operating with ${\nabla }^2$ in H yields [30]

$$\partial I/\partial \tau =-(p^2/2m+U)I+({\hslash }^2/2m)\{(2i/\hslash )\overrightarrow{p}·\nabla I+{\nabla }^{2}I\}$$

Note that $\chi =exp\left[\tau (p^2/2m+U)\right]I$ satisfies the equation

$$\partial \chi /\partial \tau =(1/2m)exp\left[\tau (p^2/2m+U)\right]\{(2i\hslash )\overrightarrow{p}·\nabla +{\hslash }^2{\nabla }^2\}exp[-\tau (p^2/2m+U)]\chi$$

and its right-hand-side $\to 0$ as $\hslash \to 0$. So $\chi$ may be expanded as $\chi =1+\hslash {\chi }_1+{\hslash }^2{\chi }_2+\cdots$. Then, it is easily verified that [30]

$$\begin{array}{ccc}\hfill \partial {\chi }_1/\partial \tau & =& -(i\tau /m)\overrightarrow{p}·\nabla U\hfill \\ \hfill \partial {\chi }_2/\partial \tau & =& (1/2m)\{-\left(2i\tau \right)(\overrightarrow{p}·\nabla U){\chi }_1+2i\overrightarrow{p}·\nabla {\chi }_1-\tau {\nabla }^{2}U+{\tau }^2{(\nabla U)}^2\}\hfill \end{array}$$

Solution of these equations give ${\chi }_1=-(i{\tau }^2/2m)\overrightarrow{p}·\nabla U$ and

$${\chi }_2=({\tau }^2/8m){(\overrightarrow{p}·\nabla U)}^2+({\tau }^2/6m^2)\sum _l\sum _k{p}_l{p}_k{\partial }^{2}U/\partial x^l\partial x_k+({\tau }^2/6m){(\nabla U)}^2-({\tau }^2/4m){\nabla }^{2}U$$

Substitution of I in to Eq.(31) provides $C(\overrightarrow{r},\tau )$ as [30]

$$C(\overrightarrow{r},\tau )=C_n[{\tau }^{-3/2}\sqrt{\pi }/2]exp(-\tau U)[1+{\tau }^3({\hslash }^2/24m){(\nabla U)}^2-{\tau }^2({\hslash }^2/12m){\nabla }^{2}U]$$

(32)

Note that ${\chi }_1$ does not contribute to the above as it is proportional to $\overrightarrow{p}$. The constant $C_n=(\sqrt{2}/{\pi }^2){(m/{\hslash }^2)}^{3/2}$ (see Eq.(2)) is also incorporated. The factor ${\tau }^{-3/2}(\sqrt{\pi }/2)$ arises due to integration over p over 0 to ∞.

8.2.2. Corrected Density

With the use of Eq.(32), it is now easy to obtain quantum-corrected expressions for density and other thermodynamic quantities. For instance, Eq.(29) readily provides $\rho \left(\overrightarrow{r}\right)$ as

$$\begin{array}{ccc}\hfill \rho \left(\overrightarrow{r}\right)& =& C_n(\sqrt{\pi }/2){\left(2\pi \right)}^{-1}{\int }_{-\infty }^{\infty }dt/{\left(it\right)}^{5/2}exp\left[it(\mu -U)\right]\hfill \\ & \times & [1+({\hslash }^2/24m)\{{\left(it\right)}^3{(\nabla U)}^2-2{\left(it\right)}^2{\nabla }^{2}U\}]{\displaystyle \frac{(\pi t/\beta )}{sinh(\pi t/\beta )}}.\hfill \end{array}$$

(33)

The steps leading to Eq.(33) show that the first term 1 inside the square bracket (devoid of $\nabla U$ and ${\nabla }^{2}U$) must provide the density $n_{tf}$ in TF model (for the potential U), and is given by

$$\begin{array}{ccc}\hfill n_{tf}(\overrightarrow{r},\mu )& =& C_n(\sqrt{\pi }/2){\left(2\pi \right)}^{-1}{\int }_{-\infty }^{\infty }dt{\left(it\right)}^{-5/2}exp\left[it(\mu -U)\right]{\displaystyle \frac{(\pi t/\beta )}{sinh(\pi t/\beta )}}\hfill \\ & =& C_n{\beta }^{-3/2}{I}_{1/2}\left[{\eta }_{tf}\right],{\eta }_{tf}=\beta (\mu -U).\hfill \end{array}$$

Using this definition, $\rho \left(\overrightarrow{r}\right)$ is expressed as

$$\rho \left(\overrightarrow{r}\right)=n_{tf}(\overrightarrow{r},\mu )+({\hslash }^2/24m)[n_{tf}^{‴}(\overrightarrow{r},\mu ){(\nabla U)}^2-2n_{tf}^{″}(\overrightarrow{r},\mu )\left({\nabla }^{2}U\right)]$$

(34)

Here, primes on $n_{tf}(\overrightarrow{r},\mu )$ denotes derivative with respect to $\mu$. This generalizes Schwinger’s expression [13] for zero-temperature.

8.2.3. Corrected Free Energy Density

As the free particle entropy is given by $S=-k_B{\sum }_{n-all}\{{\rho }_{n}ln{\rho }_n+(1-{\rho }_n)ln(1-{\rho }_n)\}$, where ${\rho }_n={\{1+exp\left[\beta ({ϵ}_n-\mu )\right]\}}^{-1}$ is the occupation probability, the free energy is $F={\sum }_{n-all}{ϵ}_n{\rho }_n-TS$. In fact, substitution of ${\rho }_n$ and a bit of simplification yield $F=\mu {\sum }_{n-all}{\rho }_n-{\beta }^{-1}{\sum }_{n-all}ln\{1+exp\left[\beta (\mu -{ϵ}_n)\right]\}$. Note the change of sign in the exponential term. This is convertible to the expression for free energy, given earlier (see Section 8.1) in this appendix, on replacing the sum by integration over $\overrightarrow{p}$, followed by a partial integration. However, the quantum-corrected free energy density is expressed as

$$f\left(\overrightarrow{r}\right)=\mu \sum _{n-all}2{\psi }_n^{\ast }\left(\overrightarrow{r}\right){\psi }_n\left(\overrightarrow{r}\right){\rho }_n-{\beta }^{-1}\sum _{n-all}2{\psi }_n^{\ast }\left(\overrightarrow{r}\right){\psi }_n\left(\overrightarrow{r}\right)ln\{1+exp\left[\beta (\mu -{ϵ}_n)\right]$$

(35)

Next, introducing the delta function $\delta (ϵ-{ϵ}_n)$ and free particle states $exp(i\overrightarrow{p}·\overrightarrow{r}/\hslash )$ to compute the trace (see Eq.(29) ) yield

$$f\left(\overrightarrow{r}\right)=\mu \rho \left(\overrightarrow{r}\right)-{\beta }^{-1}{\int }_{-\infty }^{\infty }dϵ{\left(2\pi \right)}^{-1}{\int }_{-\infty }^{\infty }dtC(\overrightarrow{r},it)exp\left(itϵ\right)ln\{1+exp\left[\beta (\mu -ϵ)\right].$$

(36)

The first sum is replaced with $\rho \left(\overrightarrow{r}\right)$, which is already obtained in Eq.(34). Two integration by parts converts the logarithmic term to $\left(\beta \right){[2cosh(\beta (ϵ-\mu )/2)]}^{-2}$. Then, a Fourier transformation of the cosh term provide

$$f\left(\overrightarrow{r}\right)=\mu \rho \left(\overrightarrow{r}\right)-{\left(2\pi \right)}^{-1}{\int }_{-\infty }^{\infty }\{dt/{\left(it\right)}^2\}C(\overrightarrow{r},it)exp\left(it\mu \right){\displaystyle \frac{(\pi t/\beta )}{sinh(\pi t/\beta )}}.$$

(37)

On introducing $C(\overrightarrow{r},it)$ yields

$$\begin{array}{ccc}\hfill f\left(\overrightarrow{r}\right)& =& \mu \rho \left(\overrightarrow{r}\right)-C_n(\sqrt{\pi }/2){\left(2\pi \right)}^{-1}{\int }_{-\infty }^{\infty }dt/{\left(it\right)}^{7/2}exp\left[it(\mu -U)\right]\hfill \\ & \times & [1+({\hslash }^2/24m)\{{\left(it\right)}^3{(\nabla U)}^2-2{\left(it\right)}^2{\nabla }^{2}U\}]{\displaystyle \frac{(\pi t/\beta )}{sinh(\pi t/\beta )}}.\hfill \end{array}$$

As done earlier in the case of $\rho \left(\overrightarrow{r}\right)$ , the terms devoid of $\nabla U$ and ${\nabla }^{2}U$ is identified as proportional to TF-kinetic energy $e_{tf}$, for the specified potential U, which is given by

$$\begin{array}{ccc}\hfill e_{tf}(\overrightarrow{r},\mu )& =& (3/2)C_n(\sqrt{\pi }/2){\left(2\pi \right)}^{-1}{\int }_{-\infty }^{\infty }dt{\left(it\right)}^{-7/2}exp\left[it(\mu -U)\right]{\displaystyle \frac{(\pi t/\beta )}{sinh(\pi t/\beta )}}\hfill \\ & =& C_n{\beta }^{-5/2}{I}_{3/2}\left[{\eta }_{tf}\right],{\eta }_{tf}=\beta (\mu -U).\hfill \end{array}$$

Then, $f\left(\overrightarrow{r}\right)$ is expressed as

$$f\left(\overrightarrow{r}\right)=\mu \rho \left(\overrightarrow{r}\right)-(2/3)e_{tf}(\overrightarrow{r},\mu )+{\displaystyle \frac{2}{3}}{\displaystyle \frac{{\hslash }^2}{24m}}[e_{tf}^{‴}\left(\mu \right){(\nabla U)}^2-2e_{tf}^{″}\left(\mu \right)\left({\nabla }^{2}U\right)]$$

(38)

The primes on $e_{tf}$ indicates derivative with $\mu$ here also. The chemical potential $\mu$ and energy parameter ${\eta }_{tf}$ corresponds to TF model, where as the free energy density is

$$f_{tf}\left(\overrightarrow{r}\right)=\mu n_{tf}(\overrightarrow{r},\mu )-(2/3)e_{tf}(\overrightarrow{r},\mu )=C_n{\beta }^{-5/2}\{{\eta }_{tf}{I}_{1/2}\left({\eta }_{tf}\right)-(2/3)I_{3/2}\left({\eta }_{tf}\right)\}+\rho \left(\overrightarrow{r}\right)n_{tf}\left(\overrightarrow{r}\right).$$

This form is already derived earlier in Section 8.1 .

8.2.4. QSM-Free Energy Density

For simplifying the expression for $f\left(\overrightarrow{r}\right)$, and for developing QSM further, density $n\left(\overrightarrow{r}\right)\equiv \rho \left(\overrightarrow{r}\right)$, including quantum corrections, is re-defined as

$$n\left(\overrightarrow{r}\right)=C_n{\beta }^{-3/2}{I}_{1/2}\left[{\eta }^{\ast }\right].$$

(39)

The new energy parameter (to be determined) is ${\eta }^{\ast }={\eta }_{tf}+{\eta }_c$, with ${\eta }_c$ representing a small correction ($\sim {\hslash }^2$) due to quantum effects. Employing the approximation $I_{3/2}\left[{\eta }_{tf}\right]=I_{3/2}[{\eta }^{\ast }-{\eta }_c]\approx I_{3/2}\left[{\eta }^{\ast }\right]+(3/2)I_{1/2}\left[{\eta }^{\ast }\right](-{\eta }_c)$, TF-free energy density reduces to $f_{tf}=C_n{\beta }^{-5/2}\{{\eta }^{\ast }{I}_{1/2}\left[{\eta }^{\ast }\right]-(2/3)I_{3/2}\left[{\eta }^{\ast }\right]\}+n\left(\overrightarrow{r}\right)U\left(\overrightarrow{r}\right)$. Furthermore, it is easily found that $e_{tf}^{‴}\sim {\beta }^3{I}_{3/2}^{‴}=(-3/8){\beta }^3{I}_{-3/2}$ as well as $e_{tf}^{″}\sim {\beta }^2{I}_{3/2}^{″}=(3/4){\beta }^2{I}_{-1/2}$. The primes on $I_{3/2}$ denote derivatives with ${\eta }^{\ast }$ and hence the factor ${\beta }^3$. With these changes, Eq.(38) reduces to

$$f\left(\overrightarrow{r}\right)=f_{tf}\left(\overrightarrow{r}\right)+{\displaystyle \frac{2}{3}}{\displaystyle \frac{{\hslash }^2}{24m}}{C}_n{\beta }^{-3/2}[-{\displaystyle \frac{3}{8}}{I}_{-3/2}{(\beta \nabla U)}^2-{\displaystyle \frac{3}{2}}{I}_{-1/2}\left(\beta {\nabla }^{2}U\right)]$$

(40)

The energy parameter in this expression is taken ${\eta }^{\ast }$, in lieu of ${\eta }_{tf}$, to the same order $\left({\hslash }^2\right)$ of accuracy. The gradient terms in U are now expressed in terms of $\nabla n$. Now, Eq.(39) yields $\nabla n=C_n{\beta }^{-3/2}{I}_{1/2}^{\prime }\nabla {\eta }^{\ast }=C_n{\beta }^{-3/2}(1/2)I_{-1/2}[-\beta \nabla U]$, where it is assumed that $\nabla {\eta }^{\ast }\sim \nabla {\eta }_{tf}=-\beta \nabla U$. In a similar manner ${\nabla }^{2}n=C_n{\beta }^{-3/2}(-1/4)I_{-3/2}{(\beta \nabla U)}^2+C_n{\beta }^{-3/2}(1/2)I_{-1/2}(-\beta {\nabla }^{2}U)$. These relations provide

$$\begin{array}{ccc}\hfill C_n{\beta }^{-3/2}{I}_{-3/2}{(\beta \nabla U)}^2& =& 4\left\{{\displaystyle \frac{I_{1/2}{I}_{-3/2}}{I_{-1/2}^2}}\right\}{\displaystyle \frac{{(\nabla n)}^2}{n}}\hfill \\ \hfill -{\displaystyle \frac{1}{2}}{C}_n{\beta }^{-3/2}{I}_{-1/2}\left(\beta {\nabla }^{2}U\right)& =& \left\{{\displaystyle \frac{I_{1/2}{I}_{-3/2}}{I_{-1/2}^2}}\right\}{\displaystyle \frac{{(\nabla n)}^2}{n}}+{\nabla }^{2}n\hfill \end{array}$$

Finally, $f\left(\overrightarrow{r}\right)$ is given by

$$f\left(\overrightarrow{r}\right)=f_{tf}\left(\overrightarrow{r}\right)+[\sigma \left({\eta }^{\ast }\right){\displaystyle \frac{{(\nabla n)}^2}{n}}+{\displaystyle \frac{{\hslash }^2}{12m}}{\nabla }^{2}n\left)\right],\sigma \left({\eta }^{\ast }\right)=\{{\displaystyle \frac{I_{1/2}{I}_{-3/2}}{I_{-1/2}^2}}\}{\displaystyle \frac{{\hslash }^2}{24m}}$$

(41)

The term ${\nabla }^{2}n$ does not contribute to total free energy as $\nabla n$ vanishes on the cell boundary and has finite derivative at the origin. Thus the gradient term $f_{grad}=\sigma \left({\eta }^{\ast }\right)[{(\nabla n)}^2/n]$ and the definition of $\sigma$ in Eq.(1) follows.

8.3. Exchange Free-Energy

Exchange contribution arises as a pure quantum effect from the symmetry requirement of wave function with respect to exchange of electron’s positions. The two-electron wave function (with spin and position) has to be anti-symmetric on exchange. The electrostatic interaction energy, as given in Eq.(1) (last term), is based on considering elections as particles obeying classical statistics, and the factor (1/2) in front is only because the term inside the integral is for a pair of electrons. In the free-gas model (which is used to compute the exchange effect), the two-electron wave function, with one at ${\overrightarrow{r}}_1$ and the other at ${\overrightarrow{r}}_2$, is ${\psi }_2={\mathcal{V}}^{-1}exp[i{\overrightarrow{k}}_1·{\overrightarrow{r}}_1]exp[i{\overrightarrow{k}}_2·{\overrightarrow{r}}_2]$ where ${\overrightarrow{k}}_1$ and ${\overrightarrow{k}}_2$ represent their momenta ($\overrightarrow{p}=\hslash \overrightarrow{k}$), respectively. Box-normalization of a free-particle state ${\mathcal{V}}^{-1/2}exp[i{\overrightarrow{k}}_1·{\overrightarrow{r}}_1]$ is implied here. The classical energy $e^2/|{\overrightarrow{r}}_1-{\overrightarrow{r}}_1|$ of a pair of electrons now is generalized as ${ϵ}_s=\int d{\overrightarrow{r}}_1\int d{\overrightarrow{r}}_2{\psi }_2^{\ast }[e^2/|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2]{\psi }_2$. However, the total wave function should be considered as ${\sigma }_A{\psi }_2$ where ${\sigma }_A$ is the (normalized) anti-symmetric spin component. There is also the possibility of wave function ${\sigma }_S{\psi }_{2A}$ (with symmetric spin component ${\sigma }_S$ ) where the anti-symmetric space-part is ${\psi }_{2A}=2^{-1/2}{\mathcal{V}}^{-1}\{exp[i{\overrightarrow{k}}_1·{\overrightarrow{r}}_1]exp[i{\overrightarrow{k}}_2·{\overrightarrow{r}}_2]-exp[i{\overrightarrow{k}}_1·{\overrightarrow{r}}_2]exp[i{\overrightarrow{k}}_2·{\overrightarrow{r}}_1]\}$. The factor $2^{-1/2}$ in front properly normalize ${\psi }_{2A}$. Energy of a pair in the latter state is ${ϵ}_a=\int d{\overrightarrow{r}}_1\int d{\overrightarrow{r}}_2{\psi }_{2A}^{\ast }[e^2/|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2\left|\right]{\psi }_{2A}$. Difference of the energy of the two states (in which the spin states are, respectively, ${\sigma }_A$ and ${\sigma }_S$) is called the exchange energy (per pair)

$$\begin{array}{c}\hfill {ϵ}_{ex}=-(1/2){\mathcal{V}}^{-2}\int d{\overrightarrow{r}}_1\int d{\overrightarrow{r}}_2[e^2/|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2]\{exp[-i{\overrightarrow{k}}_1·({\overrightarrow{r}}_1-{\overrightarrow{r}}_2)+i{\overrightarrow{k}}_2·({\overrightarrow{r}}_1-{\overrightarrow{r}}_2)]+C.C.\}\\ \hfill =-\int d{\overrightarrow{r}}_1\int d{\overrightarrow{r}}_2[e^2/|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2]RealPart\left\{G_{\overrightarrow{k_1}}^{\ast }({\overrightarrow{r}}_1-{\overrightarrow{r}}_2)G_{\overrightarrow{k_2}}({\overrightarrow{r}}_1-{\overrightarrow{r}}_2)\right\}\end{array}$$

where $G_{\overrightarrow{k}}\left(\overrightarrow{r}\right)={\mathcal{V}}^{-1}exp[i\overrightarrow{k}·\overrightarrow{r}]$ and $C.C.$ denotes complex conjugate. It is necessary to take ${ϵ}_{ex}={ϵ}_a-{ϵ}_s$ so that it decreases the (classical) repulsive energy of a pair, as required for fermions. To go over to quantum statistics in phase space, classical occupation probability (at zero-temperature), viz. $(d\overrightarrow{r}/\mathcal{V})$ is replaced as $[2/{(2\pi \hslash )}^3]d\overrightarrow{r}d\overrightarrow{p}=[2/{\left(2\pi \right)}^3]d\overrightarrow{r}d\overrightarrow{k}$. Further, the Fermi-Dirac function $f_k\left({\eta }_f\right)$, for finite-temperature occupation probability at momentum $\overrightarrow{k}$ is introduced. Thus, the definition of $G_{\overrightarrow{k}}\left(\overrightarrow{r}\right)$ is changed to $G_{\overrightarrow{k}}\left(\overrightarrow{r}\right)=[2/{\left(2\pi \right)}^3]exp[i\overrightarrow{k}·\overrightarrow{r}]f_k\left({\eta }_f\right)$, where ${\eta }_f=\beta \mu$ is the energy parameter in free-gas model. Now, note that the factor [$Realpart\left\{G_{\overrightarrow{k_1}}^{\ast }{G}_{\overrightarrow{k_2}}\right\}]$ corresponds to 2 pairs of electrons; two electrons at $(\overrightarrow{r_1},\overrightarrow{k_1})$ and other two at $(\overrightarrow{r_2},\overrightarrow{k_2})$. Therefore, the factor [$(1/4)Realpart\left\{G_{\overrightarrow{k_1}}^{\ast }{G}_{\overrightarrow{k_2}}\right\}]$ should be used for the exchange energy per electron. Integration over $\overrightarrow{k}$ changes $G_{\overrightarrow{k}}$ to

$$\tilde{G}(\overrightarrow{r},{\eta }_f)=[2/{\left(2\pi \right)}^3]\int d\overrightarrow{k}{f}_k\left({\eta }_f\right)exp[i\overrightarrow{k}·\overrightarrow{r}],f_k\left({\eta }_f\right)={[exp(\beta {\hslash }^2{k}^2/2m-{\eta }_f)+1]}^{-1}$$

which is real as $f_k\left({\eta }_f\right)$ is a function of $|\overrightarrow{k}|$. Now, the factor $(1/4){[\tilde{G}\left(({\overrightarrow{r}}_1-{\overrightarrow{r}}_2)\right),{\eta }_f]}^2$ corresponds to density of electrons $n\left(\overrightarrow{r_1}\right)$ at $\overrightarrow{r_1}$ and $n\left(\overrightarrow{r_2}\right)$ at $\overrightarrow{r_2}$. Therefore, the total exchange energy is [31]

$$F_{ex}\left({\eta }_f\right)=-(1/4)\int d{\overrightarrow{r}}_1\int d{\overrightarrow{r}}_2[e^2/|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2\left|\right]{\left[\tilde{G}({\overrightarrow{r}}_1-{\overrightarrow{r}}_2,{\eta }_f)\right]}^2$$

(42)

Change of the integral over ${\overrightarrow{r}}_2$ to integral over the relative co-ordinate $\overrightarrow{r}={\overrightarrow{r}}_1-{\overrightarrow{r}}_2$ shows that $F_{ex}\left({\eta }_f\right)=-(1/4)\int d{\overrightarrow{r}}_1\int d\overrightarrow{r}[e^2/|\overrightarrow{r}\left|\right]{\tilde{G}}^2(r,{\eta }_f)$. Integration of the angle variables yield $\tilde{G}(r,{\eta }_f)=[1/{\pi }^2]{\int }_0^{\infty }dk{k}^2{f}_k\left({\eta }_f\right)[sin\left(kr\right)/\left(kr\right)]$. It is possible to simplify the expression for $F_{ex}$ further. To that end, its derivative with ${\eta }_f$ is obtained as

$$F_{ex}^{\prime }\left({\eta }_f\right)=-(1/4)\int d{\overrightarrow{r}}_1\int d\overrightarrow{r}[e^2/r]\tilde{G}(r,{\eta }_f)[2\partial \tilde{G}(r,{\eta }_f)/\partial {\eta }_f]$$

This is expressed in terms of $I_{1/2}^{\prime }\left({\eta }_f\right)$ using an elegant method [31]. Note that with $\nu =(2m/\beta {\hslash }^2)$, $\partial f_k/\partial {\eta }_f=-(\nu /2k)[\partial f_k/\partial k]$. Substitution of this followed by an integration by parts over k yields $\partial \tilde{G}/\partial {\eta }_f=[1/{\pi }^2](\nu /2){\int }_0^{\infty }dk{f}_k[cos\left(kr\right)]$. Therefore, integration over the angle variable yields $F_{ex}^{\prime }\left({\eta }_f\right)=-[e^2/4{\pi }^4]\left(4\pi \right)\nu \int d{\overrightarrow{r}}_1{\int }_0^{\infty }dr[-\partial g\left(r\right)/\partial r]g\left(r\right)$ where $g\left(r\right)={\int }_0^{\infty }dk{f}_k[cos\left(kr\right)]$. Integral over r is evaluated as $(1/2){\left[g\left(0\right)\right]}^2$. Therefore, $F_{ex}^{\prime }\left({\eta }_f\right)=-[e^2/{\pi }^3](\nu /2)\int d{\overrightarrow{r}}_1{\left[g\left(0\right)\right]}^2$. Using the boundary condition $f_k\left({\eta }_f\right)\to 0$ as ${\eta }_f\to -\infty$, integration over ${\eta }_f$ yields

$$F_{ex}\left({\eta }_f\right)=-(e^2/{\pi }^3)(\nu /2)\int d{\overrightarrow{r}}_1{\int }_{-\infty }^{{\eta }_f}d\eta [{\int }_0^{\infty }dk{f}_k\left(\eta \right){]}^2$$

where the the definition of $g\left(0\right)$ is used. Now, it is easily verified that ${\int }_0^{\infty }dk{f}_k\left({\eta }_f\right)=\sqrt{\nu }(1/2)I_{-1/2}\left({\eta }_f\right)=\sqrt{\nu }{I}_{1/2}^{\prime }\left({\eta }_f\right)$. Thus, substitution of $\nu$ provides the final expression, viz.

$$F_{ex}\left({\eta }_f\right)=-(2/{\pi }^3)e^2{(m/{\hslash }^2)}^2{\beta }^{-2}\int d{\overrightarrow{r}}_1{\int }_{-\infty }^{{\eta }_f}d\eta {\left[I_{1/2}^{\prime }\left(\eta \right)\right]}^2.$$

The same expression with general ${\eta }^{\ast }$, in lieu of ${\eta }_f$, is used in Eq.(1).

8.4. Exchange-Correlation Free-Energy

The formula for exchange free energy derived in the previous section assumes free-electron wave functions. However, effects of electron-electron repulsion (on wave functions), within the uniform electron gas model, is needed in the free energy functional of QSM. The exchange effect when the wave functions also account for repulsive interactions is termed as exchange-correlation free-energy $F_{xc}(n,T)$. Quantum ab initio simulations are useful to compute this quantity. These use the standard random sampling techniques, however, the functional that is minimized is based on wave functions obeying symmetry requirements. Alternatively, the techniques based on path-integral approach compute averages like $\int d{\overrightarrow{r}}^N<R^{N}exp[-\beta H]R^N>$ (where $|R^N>$ denotes n-particle wave function dependent on co-ordinates ${\overrightarrow{r}}^N$) to obtain quantum mechanical partition function and thermodynamic properties.

First of all, note that the (configuration) free energy, of N particles in volume V, due to a inter-particle potential $\lambda u\left(r\right)$ is given by $exp[-\beta F\left(\lambda \right)]=Z_p=\int d{\overrightarrow{r}}^{N}exp[-\beta {\sum }_{j<i}\lambda u\left(\right|{\overrightarrow{r}}_i-{\overrightarrow{r}}_j\left|\right)\left)\right]$, where $\int d{\overrightarrow{r}}^N[\cdots ]$ denotes N-dimensional configuration integral. This is the definition of free-energy $F\left(\lambda \right)$ in terms of the partition function $Z_p$. Differentiation by $\lambda$ gives $\partial F/\partial \lambda =[N(N-1)/2]\int \int d{\overrightarrow{r}}_{1}d{\overrightarrow{r}}_{2}u\left(\right|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2\left|\right){\rho }_2\left(\right|{\overrightarrow{r}}_1-{\overrightarrow{r}}_2|,\lambda )$ where ${\rho }_2(r,\lambda )$ is the distribution function for two-particles separated by r. The pair-distribution function is $g(r,\lambda )$, which tends to 1 and $r\to \infty$, is related to ${\rho }_2$ as ${\rho }_2(r,\lambda )=(1/V^2)g(r,\lambda )$. On changing integration over ${\overrightarrow{r}}_2$ to $\overrightarrow{r}={\overrightarrow{r}}_1-{\overrightarrow{r}}_2$ and taking the thermodynamic limit yields $\partial F/\partial \lambda =N(n/2)\int d\overrightarrow{r}u\left(r\right)g(r,\lambda )$. Further, integration over $\lambda$ gives free energy per electron as $F/N=(n/2){\int }_0^{1}d\lambda \int d\overrightarrow{r}u\left(r\right)g(r,\lambda )$. This formula is applicable to the UEG model if $u\left(r\right)=e^2/r$, the Coulomb energy, and $g(r,\lambda )$ is obtained incorporating quantum effects. On subtracting the classical interaction energy per particle, viz. $(1/2)\int d\overrightarrow{r}u\left(r\right)$, the remaining exchange-correlation free energy per particle is obtained as $f_{xc}^{\ast }=(n/2){\int }_0^{1}d\lambda \int d\overrightarrow{r}u\left(r\right)h(r,\lambda )$. Here, $h(r,\lambda )=\left[g\right(r,\lambda )-1]$ is known as the pair correlation function. This is called the coupling-parameter integration for $f_{xc}^{\ast }$. It is the quantity $e_{in}(n,\lambda )=(n/2)\int d\overrightarrow{r}u\left(r\right)h(r,\lambda )$ that is computed via quantum ab initio simulations.

The integration variable $\lambda$ is replaceable with inter-particle distance $r_s$ defined as $(4\pi /3)r_s^3=n^{-1}$. Now, using $u=e^2/r$, first rewrite free energy as $f_{xc}^{\ast }={\int }_0^{1}d\lambda e_{in}\left(r_s\lambda \right)$ where $e_{in}(r_s,\lambda )=(3/2)(1/r_s^3)e^2\int drrh(r,\lambda )$. As $h(r,\lambda )$ is a functional of $\lambda u\left(r\right)=\lambda e^2/r$, it follows that $h(r,1)=h(\lambda r,\lambda )$. Multiplying and dividing by ${\lambda }^2$ yields $e_{in}(r_s,\lambda )=(1/{\lambda }^2)e_{in}(r_s,1)$ and so $f_{xc}^{\ast }={\int }_0^1(d\lambda /{\lambda }^2)e_{in}(r_s,1)$. A change of integration variable $\lambda \to r_s^{\prime }=\lambda r_s$ yields $f_{xc}^{\ast }=(1/r_s^2){\int }_0^{r_s}d{r}_s^{\prime }\left[r_s^{\prime }{e}_{in}(r_s^{\prime },1)\right]$. Differentiation with $r_s$ yields the relation $2f_{xc}^{\ast }+r_s\partial f_{xc}^{\ast }/\partial r_s=e_{in}(r_s,1)$. The quantum ab initio simulations generate $e_{in}(r_s,1)$ for a series of $r_s$ thereby determining $f_{xc}^{\ast }$ using the above formula. Using extensive data sets of $e_{in}(r_s^{\prime },1)$, generated for the UEG [19], exchange-correlation free energy per electron is fitted as

$$f_{xc}^{\ast }=-(1/r_s){\displaystyle \frac{[G_1+G_2{\sqrt{r}}_s+G_3{r}_s]}{[1+G_4{\sqrt{r}}_s+G_5{r}_s]}}.$$

(43)

The parameters $G_1,G_2,G_2,G_3,G_{4}and{G}_5$ are functions of $\theta =T/T_F$ where the Fermi-temperature $T_F$ is defined as $k_B{T}_F=({\hslash }^2/2m){\left(3{\pi }^{2}n\right)}^{2/3}$. Note that the right-hand-side is the Fermi-energy.

$$\begin{array}{ccc}\hfill G_1& =& 0.610887tanh\left[{\theta }^{-1}\right]{\displaystyle \frac{(0.75+3.04363{\theta }^2-0.09227{\theta }^3+1.7035{\theta }^4)}{(1+8.31051{\theta }^2+5.1105{\theta }^4)}}\hfill \\ \hfill G_2& =& tanh\left[{\theta }^{-1/2}\right]{\displaystyle \frac{(0.3436902+7.8215953{\theta }^2+0.3004839{\theta }^4)}{(1+15.8443467{\theta }^2+0.706281{\theta }^4)}}\hfill \\ \hfill G_4& =& tanh\left[{\theta }^{-1/2}\right]{\displaystyle \frac{(0.72700876+2.38264734{\theta }^2+0.30221237{\theta }^4)}{(1+4.39347718{\theta }^2+0.729951339{\theta }^4)}}\hfill \\ \hfill G_5& =& tanh\left[{\theta }^{-1}\right]{\displaystyle \frac{(0.25388214+0.815795138{\theta }^2+0.064684441{\theta }^4)}{(1+15.09846204{\theta }^2+0.230761357{\theta }^4)}}\hfill \\ \hfill G_3& =& G_5(0.8759442-0.2301308435exp[-1/\theta ]);\hfill \end{array}$$

(44)

The exchange-correlation energy (per electron) $e_{xc}^{\ast }=[\delta \left(\beta f^{\ast }\right)/\delta \beta ]$ is to be obtained numerically. Similar is the case for the effective potential energy (per electron) $\delta \left[n{f}_{xc}^{\ast }\right]/\delta n$.

The approach based on the liquid state theory also provides good results. If $h_p$ and $h_a$ denote the correlation functions for parallel and anti-parallel spin orientations, then $h(r,\lambda )=(1/2)[h_p(r,\lambda )+h_a(r,\lambda )]$. Perrot and Dharma-wardana [20] developed a mapping of the UEG model to a classical two-component fluid model (based on Ornstein-Zernike equation) for computing $h(r,\lambda )$. Extensive database, as a function of n and T, generated using this method is then converted to fitting formulas for $f_{xc}^{\ast }$ as well as the exchange-correlation energy $e_{xc}^{\ast }$ (provided as supplementary material). In a recent paper, Faussurier [21] reported the formula of Perrot and Dharma-wardana for $e_{xc}^{\ast }$. To provide this formula, first of all define three functions, viz. $T_1=0.610887tanh\left[{\theta }^{-1}\right]$, $T_2={[(e^2/k_B)/\left(r_s{T}_F\right)]}^{1/2}tanh\left[{\theta }^{-1/2}\right]$ and $T_3=[(e^2/k_B)/\left(r_s{T}_F\right)]tanh\left[{\theta }^{-1}\right]$. Five fitting functions are introduced as

$$\begin{array}{ccc}\hfill A& =& T_1(0.75+3.04363{\theta }^2-0.09227{\theta }^3+1.7035{\theta }^4)/(1+8.31051{\theta }^2+5.1105{\theta }^4)\hfill \\ \hfill B& =& T_2(0.341308+12.070873{\theta }^2+1.148889{\theta }^4)/(1+10.495346{\theta }^2+1.326623{\theta }^4)\hfill \\ \hfill D& =& T_2(0.614925+16.996055{\theta }^2+1.489056{\theta }^4)/(1+10.109350{\theta }^2+1.221840{\theta }^4)\hfill \\ \hfill E& =& T_3(0.539409+2.522206{\theta }^2+0.178484{\theta }^4)/(1+2.555501{\theta }^2+0.146319{\theta }^4)\hfill \\ \hfill C& =& E[0.872496+0.025248exp\{-1/\theta \left\}\right]\hfill \end{array}$$

With these parameters, energy per electron is expressed as

$$e_{xc}^{\ast }=-(1/r_s)(A+B+C)/(1+D+E)$$

(45)

Ichimaru et al show that this expression is analytically convertible to the free energy per electron [32]. With the definition of additional functions (see Faussurier’s paper [21] also)

$$\begin{array}{ccc}\hfill S& =& \sqrt{4E-D^2},D_1=C/E,Y=B-D_{1}D,D_2=2Y/E,Z=A-D_1\hfill \\ \hfill D_3& =& [Z-DY/E]ln[E+D+1]/E,D_4=2[DZ+Y\{2-D^2/E\}]/\left(ES\right)\hfill \\ \hfill D_5& =& arctan\left[\right(2E+D)/S]-arctan[D/S]\hfill \end{array}$$

the, free energy (per electron) $f_{xc}^{\ast }$ is expressed as

$$f_{xc}^{\ast }=-(1/r_s)(D_1+D_2+D_3-D_4{D}_5)$$

(46)

Note that $ln[\cdots ]$ in $D_3$ denotes natural logarithm. An alternate, but equivalent, fit ( using liquid state theory results) for $f_{xc}^{\ast }$ is provided by Perrot and Dharma-wardana [20].

8.5. Stationary Property of Free Energy

A procedure, introduced by Englert and Schwinger [14], to show the stationary property of $F(n,T)$ is extended to finite temperature in this section. The stationary property is assumed (necessary condition) to derive the Euler-Lagrange equations, but the reverse provides the sufficient condition. Further, this demonstration illustrates Homberg-Kohn-Mermin theorem for the specific functional that is postulated in QSM.

Let the triplet $(n,V,{\eta }^{\ast })$ be perturbed as $n=\rho +\chi$, $V=U+\psi$ and ${\eta }^{\ast }={\eta }^{†}+\xi$, where $(\chi ,\psi ,\xi )$ are first order $O\left(1\right)$ quantities. It is intended to show that $F(n,T)=F(\rho ,T)+O\left(2\right)$ if the triplet $(\rho ,U,{\eta }^{†})$ are chosen appropriately. As V and U need to obey the Poisson’s equations, with densities n and $\rho$, respectively, it follows that ${\nabla }^2\psi =\left(4\pi e\right)\chi$. Now, Taylor’s expansion of $I_{1/2}\left({\eta }^{\ast }\right)$ around ${\eta }^{†}$ shows that $\chi =C_n{\beta }^{-3/2}{I}_{1/2}^{\prime }\left({\eta }^{†}\right)\xi +O\left(2\right)$. Thus, only one of parameters in $(\chi ,\psi ,\xi )$ is independent. Similarly, Taylor’s expansion of $[-(2/3)I_{3/2}\left({\eta }^{\ast }\right)+{\eta }^{\ast }{I}_{1/2}\left({\eta }^{\ast }\right)]$ together with the relation between $\chi$ and $\xi$ readily yields $f_{kin}\left({\eta }^{\ast }\right)=f_{kin}\left({\eta }^{†}\right)+{\beta }^{-1}\chi {\eta }^{†}+O\left(2\right)$. Thus, the increment in $f_{kin}$, in changing from n to $\rho$, turns out to be $\Delta f_{kin}={\beta }^{-1}\chi {\eta }^{†}+O\left(2\right)$.

Using Taylor’s expansion, the incremental change in gradient term is expressed as $\Delta f_{grad}={\left[(\partial /\partial n)f_{grad}\right]}_{\rho }\chi +O\left(2\right)$. Similarly, the exchange-correlation term contributes the incremental change $\Delta f_{xc}={\left[(\partial /\partial n)f_{xc}\right]}_{\rho }\chi +O\left(2\right).$ Both these together is expressed as $\Delta f_{gxc}={\beta }^{-1}{\eta }_{gxc}\chi +O\left(2\right)$. (The term ${\eta }_{gxc}$, originating from exchange and gradient terms, is also made use of in Section 3). Thus, the sum of the changes is $[\Delta f_{kin}+\Delta f_{gxc}]={\left(4\pi e\beta \right)}^{-1}[{\eta }^{†}+{\eta }_{gxc}]{\nabla }^2\psi +O\left(2\right)$. The relation between $\chi$ and $\psi$ is also employed in this step.

Now, let it be assumed that $[{\eta }^{†}+{\eta }_{gxc}]=\beta [eU+\mu ]$, where $\mu$ corresponds to the $(\rho ,U,{\eta }^{†})$-system. This relation is identical to the Euler-Lagrange equation given in Eq.(9), however, the triplet $(\rho ,U,{\eta }^{†})$ occur in lieu of $(n,V,{\eta }^{\ast })$. Substitution of $[{\eta }^{†}+{\eta }_{gxc}]$ then yield $[\Delta f_{kin}+\Delta f_{gxc}]={\left(4\pi e\right)}^{-1}[eU+\mu ]{\nabla }^2\psi +O\left(2\right)$. The term involving $\mu$ is dropped as it does not contribute to the integral over the atomic cell. This is so because $\int {\nabla }^2\psi d\overrightarrow{r}=\left(4\pi e\right)\int \chi d\overrightarrow{r}=0$ since total charge $Z=\int nd\overrightarrow{r}=\int \rho d\overrightarrow{r}$ is the same. The expression for $\int d\overrightarrow{r}[\Delta f_{kin}+\Delta f_{gxc}]={\left(4\pi \right)}^{-1}\int d\overrightarrow{r}\left[U{\nabla }^2\psi \right]+O\left(2\right)$ completely cancels with the change in electrostatic interaction energy as shown below. However, Englert and Schwinger [14] chose $\psi$ so that this terms cancels the gradient term. This is readily seen to be possible after two partial integration as surface terms vanish. However, this approach equates an O(1) and O(0) quantities which is undesirable.

The interaction free energy is $\int d\overrightarrow{r}{f}_{int}\left(n\right)=-e\int d\overrightarrow{r}Vn-(e^2/2)\int d\overrightarrow{r}\int d{\overrightarrow{r}}^{\prime }n\left(\overrightarrow{r}\right)n\left({\overrightarrow{r}}^{\prime }\right)/|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|$, as provided by Eq.(1 and Eq.(7). Then, the increment is $\int d\overrightarrow{r}\Delta f_{int}=-e\int d\overrightarrow{r}[U\chi +\psi \rho ]-e^2\int d\overrightarrow{r}\int d{\overrightarrow{r}}^{\prime }[\rho \left(\overrightarrow{r}\right)\chi \left({\overrightarrow{r}}^{\prime }\right)/|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|+O\left(2\right)$, where symmetry of $|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|$ with respect to interchange of $\overrightarrow{r}$ and $\overrightarrow{r^{\prime }}$ is employed. Now, the relation $\psi =-e\int d{\overrightarrow{r}}^{\prime }\chi \left({\overrightarrow{r}}^{\prime }\right)/|\overrightarrow{r}-{\overrightarrow{r}}^{\prime }|$, which follows from the Poisson’s equation for $\psi$, shows that $\int d\overrightarrow{r}\Delta f_{int}=-{\left(4\pi \right)}^{-1}\int d\overrightarrow{r}U{\nabla }^2\psi +O\left(2\right)$ which cancels $\int d\overrightarrow{r}[\Delta f_{kin}+\Delta f_{eg}]$. Thus, the sufficient condition, $[{\eta }^{†}+{\eta }_{gxc}]=\beta [eU+\mu ]$, (i.e. the Euler-Lagrange equation) is shown to vanish the first order derivative terms, and hence, the stationary property of $F[n,T]$ in QSM. No appeal to TF model is needed, as employed in the earlier work (of Englert and Schwinger) for zero-temperature [14].

8.6. Strongly Bound Electrons

The nearly free electron concept embedded in the kinetic energy of QSM is inadequate for strongly bound electrons. For the cold isolated atom, Schwinger [12] re-derived Scott’s correction term $\left[(1/2)Z^2(e^2/a_0)\right]$, where $a_0={\hslash }^2/m{e}^2$ is the Bohr radius, to TF binding energy [$-0.7687Z^{7/3}(e^2/a_0)]$. The correction arises because of inadequate accounting of kinetic and electrostatic energies of strongly bound electrons in TF model. The derivation given below is more complete, and provides the basis for extension to finite-temperatures.

The free energy functional in Eq.(1) should be corrected with the term ${\int }_{r_0}d\overrightarrow{r}[f_{qm}\left(\overrightarrow{r}\right)-f_{tf}\left(\overrightarrow{r}\right)]$ which amounts to replacing free energy (kinetic and electrostatic) of strongly bound electrons in TF model with a quantum mechanical (QM) expression ( see Section 3 ) . The latter accounts for discrete structure of electron states within a sphere of radius $r_0$. When minimizing Eq.(1) with respect to density $n\left(\overrightarrow{r}\right)$, under the constraint ${\int }_{R}d\overrightarrow{r}\delta n\left(\overrightarrow{r}\right)=0$, it is necessary to insist that the number of states within $r_0$ is the same in both (TF and QM) models. Since entropy is directly related to number of states, it may be assumed to be the same in both models. Then, the correction term is expressed in terms of energy difference, viz., $\Delta E={\int }_{r_0}d\overrightarrow{r}[e_{qm}\left(\overrightarrow{r}\right)-e_{tf}\left(\overrightarrow{r}\right)]$ where electrons with energy $<-ϵ$ are included in each term in the integral. Here, $ϵ$ is a positive number, like binding energy of lowest bound state. The conservation condition of number of states provide the equation to fix the arbitrary value of $r_0$ which is related to $-ϵ$. For these electrons, the potential energy is nearly Coulomb, and is (approximately) given by $u_c\left(r\right)=-Z{e}^2/r+{ϵ}_0$, where ${ϵ}_0>0$ is the contribution of electrons outside radius $r_0$.

In the TF model, the number of these states is $n_{tf}=[2/{(2\pi \hslash )}^3]{\int }_{r_0}d\overrightarrow{r}{\int }_{p_m\left(r\right)}d\overrightarrow{p}$ where the momentum-integral is restricted by $p_m\left(r\right)$ because $p^2/2m+u_c\left(r\right)<-ϵ$. This inequality shows that $r_0$ is fixed as $u_c\left(r_0\right)=-ϵ$, which yields $r_0=Z{e}^2/{ϵ}^{\prime }$ with ${ϵ}^{\prime }=ϵ+{ϵ}_0$. Further, it also gives the maximum momentum $p_m\left(r\right)={\left[2m{ϵ}^{\prime }(r_0/r-1)\right]}^{1/2}$. These relations, together with the value of the integral ${\int }_0^{1}dx{x}^{1/2}{(1-x)}^{3/2}=\pi /16$, readily give the result $n_{tf}=(2/3)n^{\prime 3}$ where $n^{\prime }=Z{(e^2/2a_0{ϵ}^{\prime })}^{1/2}=Z{(E_0/2{ϵ}^{\prime })}^{1/2}$. Here, $E_0=e^2/a_0=27.22\left(eV\right)$ is twice the ground state binding energy of H in non-relativistic QM [12]. The cumulative number of states via QM model is ${\sum }_{i=1}^k\left(2i^2\right)=2[k(k+1)(2k+1)/6]$ for the Coulomb potential $u_c\left(r\right)$, where k is an integer (the principal quantum number), and the factor 2 is for spin degeneracy. As this is more than $(2/3)k^3$, the parameter $n^{\prime }$ should be chosen [12] between k and $k+1$. Therefore, $n^{\prime }$ equals $1.442,2.466$ and 3.475 for $k=1,2$ and 3, respectively, and $n^{\prime }\to k+1/2$ for large k. Choice of a specific value of k (say, k=1 for the ground state) fixes $n^{\prime }$ and hence also ${ϵ}^{\prime }=Z^2(E_0/2n^{\prime 2})$ and $r_0$.

Kinetic and electrostatic energy of $n^{\prime }$ electrons in TF model is easily obtained as ${\int }_{r_0}d\overrightarrow{r}{e}_{tf}=[2/{(2\pi \hslash )}^3]{\int }_{r_0}d\overrightarrow{r}{\int }_{p_m\left(r\right)}d\overrightarrow{p}\{p^2/2m+u_c\left(r\right)\}=-\left[Z^2{E}_0\right]n^{\prime }+(2/3)n^{\prime 3}{ϵ}_0$. Values of the integrals ${\int }_0^{1}dx{x}^{-1/2}{(1-x)}^{5/2}=5\pi /16$ and ${\int }_0^{1}dx{x}^{-1/2}{(1-x)}^{3/2}=6\pi /16$ are also used here. As explained earlier (see Section 3), contributions from the gradient and exchange-correlation terms need not be considered. The energy of electrons in k states via QM is $-{\sum }_{i=1}^k\left(2i^2\right)\{\left[Z^2{E}_0\right]/\left(2i^2\right)+{ϵ}_0\}=-\left[Z^2{E}_0\right]k+(2/3)n^{\prime 3}{ϵ}_0$. This is so because in the potential $u_c\left(r\right)$, the $i^{th}$ electron state has energy $-\left[Z^2{E}_0\right]/\left(2i^2\right)+{ϵ}_0$ which should be less than $-ϵ$, that is, $-\left[Z^2{E}_0\right]/\left(2i^2\right)<-{ϵ}^{\prime }$. The good approximation $n^{\prime }=k+1/2$ provides, for the cold free atom, ${\int }_{r_0}d\overrightarrow{r}[e_{qm}\left(\overrightarrow{r}\right)-e_{tf}\left(\overrightarrow{r}\right)]=(1/2)Z^2{E}_0$, which is Scott’s correction. Using the semi-classical WKB method of March [27], for computing energy of bound electrons in Coulomb potential, Barnes[26] also derived the same result.

For finite temperature, number of electrons must be computed using the Fermi-Dirac occupation probability $\{exp\beta [p^2/2m+u_c\left(r\right)-\mu ]+1{\}}^{-1}$. As this probability is required for $0\le r\le r_0$, where the potential energy dominates, $\mu$ is neglected and $u_c\left(r\right)$ is taken as the Coulomb term $-Z{e}^2/r$. Then, the number of electrons and energy in TF model are expressed as $N_{tf}=(2/3)n^{\prime 3}{F}_1\left({\beta }^{\prime }\right)$ and $E_{tf}=-\left[Z^2{E}_0\right]n^{\prime }{F}_2\left({\beta }^{\prime }\right)+N_{tf}{ϵ}_0$. The temperature-dependent factors (which are normalized to 1 as $T\to 0$) are easily obtained as

$$\begin{array}{ccc}\hfill F_1\left({\beta }^{\prime }\right)& =& (16/\pi ){\int }_0^1{x}^{2}dx{\int }_0^{(1-x)/x}dy(3/2)y^{1/2}f(y,x)\hfill \\ \hfill F_2\left({\beta }^{\prime }\right)& =& (16/\pi ){\int }_0^1{x}^{2}dx{\int }_0^{(1-x)/x}dy[(3/2)y^{1/2}(1/x)-(5/2)y^{3/2}]f(y,x)\hfill \end{array}$$

Here, $f(y,x)=[exp\beta {ϵ}^{\prime }\{y-1/x\}+1{]}^{-1}$. Note that $\{y-1/x\}\le -1$ ( within the integrals ), and so $f(y,x)$ approaches 1 as $T\to 0$. Further, a very plausible choice $k=2$ (i.e. ten states) shows that ${ϵ}^{\prime }=Z^2{E}_0/12.164$ so that $f(y,x)\sim 1$ even for high temperatures. This is also evident from the corresponding quantities of the QM model. These are expressed as $N_{qm}=2[k(k+1)(2k+1)/6]G_1\left(\beta \right)$ and $E_{qm}=-Z^2{E}_{0}k{G}_2\left(\beta \right)+N_{qm}{ϵ}_0$, with the normalized temperature-dependent factors

$$\begin{array}{ccc}\hfill G_1\left(\beta \right)& =& \left[\sum _{i=1}^k\left(2i^2\right){]}^{-1}\sum _{i=1}^k\left(2i^2\right)\right\{exp\beta \{-Z^2{E}_0/\left(2i^2\right)\}+1{\}}^{-1}\hfill \\ \hfill G_2\left(\beta \right)& =& {\left[k\right]}^{-1}\sum _{i=1}^k\{exp\beta \{-Z^2{E}_0/\left(2i^2\right)\}+1{\}}^{-1}\hfill \end{array}$$

The factor $\beta Z^2{E}_0$ shows that $G_1$ and $G_2$ are close to 1 even for high values of $k_{B}T$. For all numerical applications, the fitted function $\Delta E=0.46Z^2(e^2/a_0){\{exp(-2.5\ast \beta Z^2{E}_0)+1\}}^{-1}$ provide accurate results ( error $\le 4\%$) for all temperatures (see Figure 1).

8.7. Solution via Nyström’s Method

The Nyström’s method to solve Eq.(25), and the source function in Eq.(27), converts them to a matrix equation [25]. To this end, the integral term and the source function are suitably evaluated by a quadrature formula. In this work, $\varphi \left(x\right)$ (as well as $\tilde{\varphi }\left(x\right)and{\tilde{\eta }}^{\ast }\left(x\right))$ are approximated using piece-wise quadratic polynomials in the sub-intervals covering [0, 1]. For this purpose, the interval [$x_1$, 1] is divided in to N (even) sub-intervals defined with $N+1$ (odd) mesh points defined as $x_j={[j/(N+1)]}^2(j=1,2,3,\cdots N+1)$. This choice of $x_j$ yields more mesh points near $x=0$ where $\varphi \left(x\right)$ varies more steeply. There are only $N+1$ unknowns ${\varphi }_j=\varphi \left(x_j\right)$ to be determined as $\varphi \left(0\right)$ is known. In the sub-interval ${\Delta }_j=x_{j-1}\le x\le x_{j+1}$ ($j=2,4,\cdots N$), the unknown function $\varphi \left(x\right)$ is represented as a quadratic polynomial $\varphi \left(x\right)={\varphi }_{j-1}{f}_j^{j-1}\left(x\right)+{\varphi }_j{f}_j^j\left(x\right)+{\varphi }_{j+1}{f}_j^{j+1}\left(x\right)$. The Lagrange polynomials are defined as $f_j^k\left(x\right)={\prod }_{i=j-1,i\ne k}^{j+1}[(x-x_i)/(x_k-x_i)]$ (for $j-1\le k\le j+1$) so that these are normalized to unity at $x_k$ and pass through others. Similar representation is employed for $\tilde{\varphi }\left(x\right)$ and ${\tilde{\eta }}^{\ast }\left(x\right)$. Substitution of $\varphi \left(x\right)$ into Eq.(25) yields Nyström’s interpolation formula

$$\begin{array}{ccc}\hfill \varphi \left(x\right)& +& \mathcal{A}\sum _{j=1}^{N+1}{W}_j\left(x\right){\varphi }_j-{\varphi }_{N+1}x=S^{\ast }\left(x\right),\hfill \end{array}$$

(47)

$$\begin{array}{ccc}\hfill S^{\ast }\left(x\right)& =& S\left(x\right)-{\mathcal{A}}_0^{x_1}\mathcal{K}(x,y)Q\left({\tilde{\eta }}^{\ast }\right)\varphi \left(y\right)dy.\hfill \end{array}$$

(48)

The weight functions $W_j\left(x\right)$ (for $j=1,2,\cdots N+1$) are given by

$$\begin{array}{ccc}\hfill W_1\left(x\right)& =& \underset{x_1}{\overset{x_3}{ʃ}}\mathcal{K}(x,y)Q\left({\tilde{\eta }}^{\ast }\right)f_2^1\left(y\right)dy,W_{N+1}\left(x\right){=}_{x_{N-1}}^{x_{N+1}}\mathcal{K}(x,y)Q\left({\tilde{\eta }}^{\ast }\right)f_N^{N+1}\left(y\right)dy,\hfill \\ \hfill W_j\left(x\right)& =& \underset{x_{j-1}}{\overset{x_{j+1}}{ʃ}}\mathcal{K}(x,y)Q\left({\tilde{\eta }}^{\ast }\right)f_j^j\left(y\right)dy,j=2,4,6,\cdots N,\hfill \\ \hfill W_j\left(x\right)& =& \underset{x_{j-2}}{\overset{x_j}{ʃ}}\mathcal{K}(x,y)Q\left({\tilde{\eta }}^{\ast }\right)f_{j-1}^j\left(y\right)dy{+}_{x_j}^{x_{j+2}}\mathcal{K}(x,y)Q\left({\tilde{\eta }}^{\ast }\right)f_{j+1}^j\left(y\right)dy,j=3,5,7,\cdots N-1\hfill \end{array}$$

(49)

For expressing $S^{\ast }\left(x\right)$, additional weight functions $W_j^{\ast }\left(x\right)$$(j=2,4,6,\cdots N)$ are defined as

$$W_j^{\ast }\left(x\right){=}_{x_{j-1}}^{x_{j+1}}\mathcal{K}(x,y)y{I}_{1/2}\left({\tilde{\eta }}^{\ast }\right)dy,j=2,4,\cdots N.$$

(50)

Then, the source term $S^{\ast }\left(x\right)$ is expressed as

$$S^{\ast }\left(x\right)={\varphi }_0+\mathcal{A}\sum _{j=1}^{N+1}{W}_j\left(x\right){\tilde{\varphi }}_j-\mathcal{A}\sum _{j=2★★}^N{W}_j^{\ast }\left(x\right)-{\mathcal{A}}_0^{x_1}\mathcal{K}(x,y)y{I}_{1/2}\left({\eta }^{\ast }\right)dy$$

(51)

Note that the second sum (with $★★$ on lower limit), involving $W_j^{\ast }$, is only over even values of j. The last integral term is obtained by adding up all contributions in the interval $[0,x_1]$, and also using the approximation $I_{1/2}\left[{\eta }^{\ast }\right]\sim I_{1/2}\left[{\tilde{\eta }}^{\ast }\right]+(1/2)I_{-1/2}\left[{\tilde{\eta }}^{\ast }\right]({\eta }^{\ast }-{\tilde{\eta }}^{\ast })$.

Expressing ${\eta }^{\ast }\left(x\right)=\zeta \left(x\right)/x$ and expanding $I_{-1/2}[\zeta /x]\sim {[\zeta /x]}^{1/2}$ near $x\sim 0$, it is easy to show that $\zeta \left(x\right)\sim \varphi \left(0\right)+\mathcal{B}\varphi {\left(0\right)}^{1/2}{x}^{1/2}+({\varphi }^{\prime }\left(0\right)+{\mathcal{B}}^2/2)x$. The definition of $\mathcal{B}$ involving $C_{ex}^{†}$ is used here to account for the exchange contribution. Thus, $\zeta \left(y\right)$ in $(0<y<x_1)$ is approximated by the quadratic polynomial $\zeta \left(y\right)=\varphi \left(0\right)[1-{(y/x_1)}^2]+\mathcal{B}\varphi {\left(0\right)}^{1/2}{y}^{1/2}[1-{(y/x_1)}^{3/2}]+\zeta \left(x_1\right){(y/x_1)}^2+[{\varphi }^{\prime }\left(0\right)+(1/2){\mathcal{B}}^2](x_1-y)(y/x_1)$ for all computations. This polynomial passes through $\zeta \left(0\right)=\varphi \left(0\right)$ and $\zeta \left(x_1\right)$ and has the correct variation near the origin. The slope ${\varphi }^{\prime }\left(0\right)$ is given by

$${\varphi }^{\prime }\left(0\right)={\varphi }^{\prime }\left(1\right)-{\mathcal{A}}_0^{1}y{I}_{1/2}\left[{\eta }^{\ast }\left(y\right)\right]dy={\varphi }_{N+1}-{\mathcal{A}}_0^{x_1}y{I}_{1/2}\left[{\eta }^{\ast }\left(y\right)\right]dy-\mathcal{A}\sum _{j=2★★}^N{W}_j^{\ast }\left(0\right).$$

(52)

Differentiation of Eq.(23) provides the first equality. Currently available values of $\zeta \left(x_1\right)=x_1{\eta }^{\ast }\left(x_1\right)$ and ${\varphi }^{\prime }\left(0\right)$ are employed in Newton’s iterations. Furthermore, all the integrals are evaluated using low order (say 6 or 8) Gauss quadrature formula. It is to be noted that the asymptotic term $\left\{y(2/3){[\varphi \left(0\right)/y]}^{3/2}\right\}$ should be subtracted in the integral in Eq.(52) (before applying Gauss quadrature formula) and its contribution $\left[(4/3)\varphi \left(0\right)x_1^{1/2}\right]$ added separately. The steps involved in obtaining Eq.(47) are the same as in deriving three-point Simpson’s integration formula. Hence the local truncation error in the method is proportional to ${\Delta }_N^4$. In the collocation rule in Nyström’s method, Eq.(47) is enforced at the discrete set of values $x_k(k=1,2,\cdots N+1)$, which yields a matrix equation $\mathbf{M}\overrightarrow{\varphi }={\mathbf{S}}^{\ast }$ of order $(N+1)$ for the vector $\overrightarrow{\varphi }$. The coefficient matrix is $\mathbf{M}=\mathbf{I}+\mathcal{A}\mathbf{W}-\mathbf{I}\mathbf{x}$. The components ${\mathbf{W}}_{k,j}=W_j\left(x_k\right)$ and those of of the vectors ${\mathbf{S}}^{\ast }$ and $\mathbf{x}$ are $S_j^{\ast }$ and $x_j$, respectively.