Analytical Representations of Thermodynamic Functions of Thomas-Fermi Model

1. Introduction

The Thomas-Fermi (TF) model is the simplest of all orbital-free density functional theories for determining the electron-components of thermodynamic properties of materials [1,2]. Its finite-temperature generalization, together with the spherical atomic cell concept, is commonly used for computing electron-properties of hot dense matter [3]. The TF model leads to a nonlinear boundary value problem that is to be solved using numerical techniques [4]. Important improvements to the original algorithm to solve the TF equation have also been reported [5]. Addition of the quantum mechanical exchange correction to the TF model leads to Thomas-Fermi-Dirac (TFD) model, while the Thomas-Fermi-Dirac-Weizsäcker (TFDW) model is obtained when electron-density gradient corrections are accounted [6]. The quantum statistical model (QSM) belongs to the same class although its development yielded the correct coefficient of the gradient term in the energy functinal [7]. Numerical computations of electron-component of equation of state (EOS) using the QSM is somewhat involved. Moreover, there is no scaling of thermodynamic properties with atomic (Z) and mass (A) numbers [8,9]. Applications of TF model indeed demonstrate the need to account for the electron binding effects even in obtaining thermal part of electronic properties [10]. Available analytical representations [11,12] of tabulated numerical results are not sufficiently accurate in the wide density-temperature plane required for applications. Perturbation solutions starting from the hot-curve (fully ionized state) do not explicitly account for the nuclear potential in the results [13]. The excellent compilation of approximate analytical solutions of TF equation [14] reinforces the need for efficient numerical schemes.

The well known method to solve the TF equation is to rewrite it as a Voltera integral equation and thereafter formulate a discrete version and then apply the shooting method to obtain a solution [4]. However, numerical difficulties arise in the shooting method for lower temperature range, and special techniques are needed to overcome them [5]. An alternate method is to apply Newton’s method to solve the non-linear TF equation. A discrete version of the resulting (linear second order) differential equation, obtained via the finite difference method (FDM), and the double sweep method provide an efficient algorithm, as the coefficient matrix is tri-diagonal [15]. However, the second order derivative of the TF-function does not exist at the origin and is very large in its neighborhood where the FDM is not applicable. One of the aims of this paper is to develop a method with forth order accuracy to solve the linear equations which result from Newton’s method. Then, it is pointed out that Brachman’s differential equation [16] is employable as it provides directly the thermal component of electron energy. Furthermore, with the use of very accurate rational function approximations of Fermi-Dirac integrals [17], the present paper offers a simple algorithm that works in the entire scaled temperature-density plane.

The paper is organized as follows. The relevant formulas related to the TF model are collected in thomas-fermi-model, while the numerical scheme is detailed numerical-scheme. Analytical representations of thermodynamic functions, which include electron energy, pressure, ionization, Fermi energy and initial slope of TF-function, are provided and discussed in analytical-representation. Accuracy of these representations is demonstrated, in applications, via computation of electron-Hugoniot and the Hugoniot of Cu and comparison with experimental (or theoretical) data up to about 20.4 TPa pressure. Finally, summary provides a summary of the work.

2. Thomas-Fermi Model

The finite-temperature TF model follows from the Poisson equation for electrostatic potential in an atomic cell, where the electron density is given by the Fermi-Dirac distribution [3]. In scaled variables, it yields a nonlinear second order differential equation [4]

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

(1)

Here, the dimensionless radial co-ordinate $x=r/R$, where R denotes the spherical cell radius. The TF-function $\varphi \left(x\right)$ and dimensionless parameter $\mathcal{A}$ are defined as

$$\begin{array}{ccc}\hfill \varphi \left(x\right)& =& (r/R)[{\mathcal{E}}_{\mathcal{F}}+eV\left(r\right)]/\left(k_{B}T\right),\hfill \\ \hfill \mathcal{A}& =& {(R/C)}^2,C=h^{3/2}/\left[4\pi e{\left(2m\right)}^{3/4}{\left(k_{B}T\right)}^{1/4}\right],\hfill \end{array}$$

where $V\left(r\right)$ is the electrostatic potential, ${\mathcal{E}}_{\mathcal{F}}$ the Fermi energy, $k_B$ Boltzmann’s constant and T temperature. Note that $e$ denotes the magnitude of electronic charge so that $-eV$ is the electron potential energy. Remaining symbols are electron mass $m$ and Planck’s constant h. The Fermi-Dirac function of order $n$ is given by

$$I_n\left(\eta \right)=\underset{0}{\overset{\infty }{\int }}dy{y}^n{[exp(y-\eta )+1]}^{-1}.$$

(2)

Accurate rational function approximations [17] for these integrals (and their inverse functions) are available for $n=-{\displaystyle \frac{1}{2}},{\displaystyle \frac{1}{2}},{\displaystyle \frac{3}{2}}$ and ${\displaystyle \frac{5}{2}}$. The limiting condition that $V\left(r\right)$ must approach the nuclear potential $Ze/r$ as $r\to 0$ yields one boundary condition $\varphi \left(0\right)=Z{e}^2/\left(R{k}_{B}T\right)$. The other follows from the condition that the atomic sphere is electrically neutral, which implies that potential gradient $V^{\prime }\left(R\right)=0$ at R. This fact, together with the choice of potential value $V\left(R\right)=0$ at R (which is arbitrary) provides the condition ${\varphi }^{\prime }\left(1\right)=\varphi \left(1\right)$. Now, note that $\varphi \left(1\right)={\mathcal{E}}_{\mathcal{F}}/k_{B}T$ which is also implicitly defined through the relation $Z^{*}/\mathcal{V}=(4\pi /h^3){\left(2m{k}_{B}T\right)}^{3/2}{I}_{1/2}[{\mathcal{E}}_{\mathcal{F}}/k_{B}T]$, where $Z^{*}$ is number of free electrons at temperature T and $\mathcal{V}$ the cell volume.

2.1. Thermodynamic Properties

Detailed derivations of all thermodynamic properties within the TF model are well known [4]. For instance, electron pressure is given by

$$P={\displaystyle \frac{2}{9}}{\displaystyle \frac{1}{\mathcal{V}}}{\displaystyle \frac{\mathcal{A}R}{e^2}}{\left(k_{B}T\right)}^2{I}_{3/2}\left[\varphi \left(1\right)\right].$$

(3)

Only $\varphi \left(1\right)$ is needed for computing P. Alternatively, if P is known, ${\mathcal{E}}_{\mathcal{F}}/k_{B}T$ and thereafter $Z^{*}$ can be obtained. Kinetic energy of electrons (per atom) needs $\varphi \left(x\right)$ over the full interval, and is expressed in the integral

$$E_{kin}={\displaystyle \frac{\mathcal{A}R}{e^2}}{\left(k_{B}T\right)}^2\underset{0}{\overset{1}{\int }}dx{x}^2{I}_{3/2}\left[{\displaystyle \frac{\varphi \left(x\right)}{x}}\right].$$

(4)

Similarly, the potential energy of electrons is given by

$$E_{pot}={\displaystyle \frac{2}{3}}{\displaystyle \frac{\mathcal{A}R}{e^2}}{\left(k_{B}T\right)}^2\{I_{3/2}\left[\varphi \left(1\right)\right]-3\underset{0}{\overset{1}{\int }}dx{x}^2{I}_{3/2}[{\displaystyle \frac{\varphi \left(x\right)}{x}}\left]\right\}.$$

(5)

Note that the last three equations (3 - 5) yield the virial theorem expressed as $2E_{kin}+E_{pot}=3P\mathcal{V}$. A partial integration in Eq.(5) together with the TF equation and the relation $I_{3/2}^{\prime }(\varphi /x)=(3/2)I_{1/2}(\varphi /x)({\varphi }^{\prime }/x-\varphi /x^2)$ yields the alternate expression

$$E_{pot}=-{\displaystyle \frac{\mathcal{A}R}{e^2}}{\left(k_{B}T\right)}^2\underset{0}{\overset{1}{\int }}dxx{I}_{1/2}\left[{\displaystyle \frac{\varphi \left(x\right)}{x}}\right]\mathcal{A}\underset{x}{\overset{1}{\int }}dy{y}^2{I}_{1/2}\left[{\displaystyle \frac{\varphi \left(y\right)}{y}}\right].$$

Asymptotic forms of $I_n[\varphi \left(x\right)/x]$ should be used near $x=0$ in evaluating these integrals. An alternate way to compute total energy of electrons $E=E_{kin}+E_{pot}$ (per atom) is the following.

The free energy F and thereafter the entropy S are obtained [16] by integrating the Gibbs-Helmholtz equation $E=-T^2\partial (F/T)/\partial T$. Entropy so obtained is given by

$$TS={\displaystyle \frac{5}{3}}{E}_{kin}+2E_{ee}+E_{eN}-Z{\mathcal{E}}_{\mathcal{F}},$$

(6)

where $E_{ee}$ and $E_{eN}$ are, respectively, the electron-electron and electron-nuclear components of potential energy $E_{pot}$. Using the definition of $E_{eN}$ and Poisson equation for V, it is readily found taht $E_{eN}=Z[k_{B}T{\varphi }^{\prime }\left(0\right)-{\mathcal{E}}_{\mathcal{F}}]$, where ${\varphi }^{\prime }\left(0\right)$ denotes the slope of TF function (with respect to x) at the origin [15]. Use of this expression, together with the virial theorem, reduces Eq.(6) to the expression $TS=-P\mathcal{V}+(7/3)E-Z{k}_{B}T{\varphi }^{\prime }\left(0\right)$. Finally, the thermodynamic relation $T{(\partial S/\partial T)}_V=T{(\partial E/\partial T)}_V$ yields Brachman’s differential equation for energy [16,18]

$$T{\displaystyle \frac{\partial E}{\partial T}}-{\displaystyle \frac{7}{4}}E={\displaystyle \frac{3}{4}}{T}^2{\displaystyle \frac{\partial }{\partial T}}\left[{\displaystyle \frac{P\mathcal{V}+Z{k}_{B}T{\varphi }^{\prime }\left(0\right)}{T}}\right],$$

(7)

where the partial derivatives with T, like $(\partial E/\partial T)$, are taken at constant volume. This first order differential equation can be integrated from a suitable starting value of T. Note that only ${\varphi }^{\prime }\left(0\right)$, and not the full TF function $\varphi \left(x\right)$, is needed in Brachman’s equation, however, the differential equation needs to be integrated.

In the high temperature limit when electron density is uniform, $eV=(Z{e}^2/R)(1/x+x^2/2-3/2)$, and consequently, $\varphi \left(x\right)=x({\mathcal{E}}_F/k_{B}T)+\varphi \left(0\right)(1+x^3/2-3x/2)$ and $k_{B}T{\varphi }^{\prime }\left(0\right)={\mathcal{E}}_F-(3/2)Z{e}^2/R$. Fermi energy reduces to ${\mathcal{E}}_{\mathcal{F}}=-k_{B}Tlog[8\pi /3]-(3/2)k_{B}Tlog[2\pi m{k}_{B}T{R}^2/\left(h^2{Z}^{2/3}\right]$ in the same limit. This limiting form yields (see Eq.(7)) the ideal gas limit $E=(3/2)Z{k}_{B}T$, as $(P\mathcal{V}/T)$ becomes independent of T in this limit. In the zero-temperature limit [18], Eq.(7) reduces to $E_0=(3/7)[P_0\mathcal{V}+(Z^2{e}^2/\mu ){\varphi }_{0in}^{\prime }]$, where ${\varphi }_{0in}^{\prime }$ is the initial slope of (the zero-temperature) TF-function, $\mu =a_0(9{\pi }^2/128)Z^{-1/3}$ is the TF scaling length, and $a_0={\overline{h}}^2/\left(m{e}^2\right)$ the Bohr radius. Substituting $E=E_0+E_{th}$, and integrating the resulting equation from T to ∞, the thermal energy $E_{th}$ of electrons is obtained as

$$E_{th}(\mathcal{V},T)=(3/4)g(\mathcal{V},T)-(9/16)T^{7/4}\underset{T}{\overset{\infty }{\int }}{t}^{-11/4}g(\mathcal{V},t)dt,$$

(8)

where $g=(P-P_0)\mathcal{V}+Z{k}_{B}T{\varphi }^{\prime }\left(0\right)-(Z^2{e}^2/\mu ){\varphi }_{0in}^{\prime }$. It is necessary to integrate from T to ∞ for the sake of numerical stability (see Eq.(20)). Further, the zero-temperature terms are subtracted in $g(\mathcal{V},t)$ so that the integral does not diverge as $T\to 0$. Note that the expression for $E_{th}$ is consistent with the high-temperature limits $P_{th}=P-P_0=Z{k}_{B}T/\mathcal{V}$ and $E_{th}=(3/2)Z{k}_{B}T$ for thermal pressure and energy of electrons in the cell. The term $-(3/2)log\left[k_{B}T\right]$, arising from ${\varphi }^{\prime }\left(0\right)$ in this limit, cancels out. Very accurate numerical approximations for zero-temperature pressure $P_0$ and slope ${\varphi }_{0in}^{\prime }$ versus cell radius are available [21,22] (see analytical-representation).

3. Numerical Scheme

To develop a numerical method, the TF equation is converted to an integral equation. On integrating from x to 1, and thereafter from 0 to x, Eq.(1) reduces to

$$\varphi \left(x\right)+\mathcal{A}\underset{0}{\overset{1}{\int }}\mathcal{K}(x,y)y{I}_{1/2}\left[{\displaystyle \frac{\varphi \left(y\right)}{y}}\right]dy-\varphi \left(1\right)x=\varphi \left(0\right),$$

(9)

where the kernel is $\mathcal{K}(x,y)=x$ for $y\ge x$ and $\mathcal{K}(x,y)=y$ for $y\le x$. The boundary conditions on $\varphi \left(x\right)$ are incorporated in Eq.(9). To employ Newton’s method, first the non-linear term is approximated using the expansion $I_{1/2}\left[\varphi \right]\sim I_{1/2}\left[\tilde{\varphi }\right]+(1/2)I_{-1/2}\left[\tilde{\varphi }\right](\varphi -\tilde{\varphi })$, which follows from Taylor’s expansion of $I_{1/2}\left[\varphi \right]$ around the function $\tilde{\varphi }$. The recurrence relation $I_n^{\prime }\left(x\right)=n{I}_{n-1}\left(x\right)$ is also used. Substitution into Eq.(9) yield linear equation

$$\varphi \left(x\right)+{\displaystyle \frac{\mathcal{A}}{2}}\underset{0}{\overset{1}{\int }}\mathcal{K}(x,y)I_{-1/2}\left[{\displaystyle \frac{\tilde{\varphi }\left(y\right)}{y}}\right]\varphi \left(y\right)dy-\varphi \left(1\right)x=S\left(x\right),$$

(10)

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

$$S\left(x\right)=\varphi \left(0\right)+\mathcal{A}\underset{0}{\overset{1}{\int }}\mathcal{K}(x,y)\left\{{\displaystyle \frac{1}{2}}{I}_{-1/2}\right[{\displaystyle \frac{\tilde{\varphi }\left(y\right)}{y}}]\tilde{\varphi }\left(y\right)-y{I}_{1/2}[{\displaystyle \frac{\tilde{\varphi }\left(y\right)}{y}}\left]\right\}dy,$$

(11)

The integrals in these equations exist for $x\ge 0$ because $I_{-1/2}[\tilde{\varphi }/y]\sim \tilde{\varphi }\left(0\right)/y^{1/2}$ and $I_{1/2}[\tilde{\varphi }/y]\sim \tilde{\varphi }\left(0\right)/y^{3/2}$ as $y\to 0$. Several numerical schemes exist [19] for solving the linear Fredholm integral equations (FIE) as in Eq.(10). Thereafter, a few iterations, together with the replacement $\tilde{\varphi }\left(y\right)\leftarrow \varphi \left(y\right)$, provides the required solution, if Newton’s method converges.

In the Nyström’s method to solve FIE, which is employed to convert the FIE to a matrix equation [20], the integral term is evaluated using a suitable quadrature. In the procedure adapted here, $\varphi \left(x\right)$ (and also $\tilde{\varphi }\left(x\right)$) is approximated using piece-wise quadratic polynomials in the sub-intervals in [0, 1]. To this end, the interval [$x_1$, 1] is divided in to N (even) sub-intervals defined with $N+1$ (odd) mesh points $x_j={[j/(N+1)]}^2(j=1,2,3,\cdots N+1)$. This choice of $x_j$ provides 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=r_{j-1}\le r\le r_{j+1}$ ($j=2,4,\cdots N$), the unknown function $\varphi \left(x\right)$ is represented as as a quadratic polynomial

$$\begin{array}{ccc}\hfill \varphi \left(x\right)& =& {\varphi }_{j-1}{f}_j^L\left(x\right)+{\varphi }_j{f}_j^M\left(x\right)+{\varphi }_{j+1}{f}_j^R\left(x\right),\hfill \\ \hfill f_j^L\left(x\right)& =& (x_j-x)(x_{j+1}-x)/\left[(x_j-x_{j-1})(x_{j+1}-x_{j-1})\right],\hfill \\ \hfill f_j^M\left(x\right)& =& (x-x_{j-1})(x_{j+1}-x)/\left[(x_j-x_{j-1})(x_{j+1}-x_j)\right],\hfill \\ \hfill f_j^R\left(x\right)& =& (x-x_{j-1})(x-x_j)/\left[(x_{j+1}-x_{j-1})(x_{j+1}-x_j)\right].\hfill \end{array}$$

(12)

Substitution of $\varphi \left(x\right)$ into Eq.(10) yields Nyström’s interpolation formula

$$\varphi \left(x\right)+{\displaystyle \frac{\mathcal{A}}{2}}\sum _{j=1}^{N+1}{W}_j\left(x\right){\varphi }_j-{\varphi }_{N+1}x=S^{*}\left(x\right),$$

(13)

where the modified source function $S^{*}\left(x\right)$ is given by

$$S^{*}\left(x\right)=S\left(x\right)-{\displaystyle \frac{\mathcal{A}}{2}}\underset{0}{\overset{x_1}{\int }}\mathcal{K}(x,y)I_{-1/2}\left[{\displaystyle \frac{\tilde{\varphi }\left(y\right)}{y}}\right]\varphi \left(y\right)dy.$$

(14)

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}{\int }}\mathcal{K}(x,y)I_{-1/2}\left[{\displaystyle \frac{\tilde{\varphi }}{y}}\right]f_2^L\left(y\right)dy,W_{N+1}\left(x\right)=\underset{x_{N-1}}{\overset{x_{N+1}}{\int }}\mathcal{K}(x,y)I_{-1/2}\left[{\displaystyle \frac{\tilde{\varphi }}{y}}\right]f_N^R\left(y\right)dy,\hfill \end{array}$$

$$\begin{array}{ccc}\hfill W_j\left(x\right)& =& \underset{x_{j-1}}{\overset{x_{j+1}}{\int }}\mathcal{K}(x,y)I_{-1/2}\left[{\displaystyle \frac{\tilde{\varphi }}{y}}\right]f_j^M\left(y\right)dy,j=2,4,6,\cdots N,\hfill \\ \hfill W_j\left(x\right)& =& \underset{x_{j-2}}{\overset{x_j}{\int }}\mathcal{K}(x,y)I_{-1/2}\left[{\displaystyle \frac{\tilde{\varphi }}{y}}\right]f_{j-1}^R\left(y\right)dy+\underset{x_j}{\overset{x_{j+2}}{\int }}\mathcal{K}(x,y)I_{-1/2}\left[{\displaystyle \frac{\tilde{\varphi }}{y}}\right]f_{j+1}^L\left(y\right)dy,j=3,5,7,\cdots N-1\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill W_j^{*}\left(x\right)& =& \underset{x_{j-1}}{\overset{x_{j+1}}{\int }}\mathcal{K}(x,y)y{I}_{1/2}\left[{\displaystyle \frac{\tilde{\varphi }}{y}}\right]dy,j=2,4,\cdots N.\hfill \end{array}$$

(16)

An additional set of functions $W_j^{*}\left(x\right)$$(j=2,4,6,\cdots N)$ are also introduced for expressing the source term $S^{*}\left(x\right)$ as

$$S^{*}\left(x\right)={\varphi }_0+{\displaystyle \frac{\mathcal{A}}{2}}\sum _{j=1}^{N+1}{W}_j\left(x\right){\tilde{\varphi }}_j-\mathcal{A}\sum _{j=2★}^N{W}_j^{*}\left(x\right)-\mathcal{A}\underset{0}{\overset{x_1}{\int }}\mathcal{K}(x,y)y{I}_{1/2}\left[{\displaystyle \frac{\varphi \left(y\right)}{y}}\right]dy$$

(17)

Note that the second sum (with ★ on lower limit) involving $W_j^{*}$ is only over even values of j. The last integral term is obtained using the approximation $I_{1/2}\left[\varphi \right]\sim I_{1/2}\left[\tilde{\varphi }\right]+(1/2)I_{-1/2}\left[\tilde{\varphi }\right](\varphi -\tilde{\varphi })$. In all computations, $\varphi \left(y\right)$ in $(0<y<x_1)$ is approximated by the quadratic polynomial $\varphi \left(y\right)=\varphi \left(0\right)[1-{(x/x_1)}^2]+{\varphi }_1{(x/x_1)}^2+{\varphi }^{\prime }\left(0\right)(x_1-x)(x/x_1)$. This polynomial passes through $\varphi \left(0\right)$ and ${\varphi }_1$ and has the slope ${\varphi }^{\prime }\left(0\right)$ at the origin. The slope at the origin given by

$$\begin{array}{ccc}\hfill {\varphi }^{\prime }\left(0\right)& =& {\varphi }^{\prime }\left(1\right)-\mathcal{A}\underset{0}{\overset{1}{\int }}y{I}_{1/2}\left[{\displaystyle \frac{\varphi \left(y\right)}{y}}\right]dy\hfill \\ & =& {\varphi }_{N+1}-\mathcal{A}\underset{0}{\overset{x_1}{\int }}y{I}_{1/2}\left[{\displaystyle \frac{\varphi \left(y\right)}{y}}\right]dy-\mathcal{A}\sum _{j=2★}^N{W}_j^{*}\left(0\right),\hfill \end{array}$$

(18)

The first line follows by integrating Eq.(1) over [0,1]. Latest available values of ${\varphi }_1=\varphi \left(x_1\right)$ and ${\varphi }^{\prime }\left(0\right)$ are employed in computer implementation. Rigorously, a term proportional to $y^{3/2}$ must be added in $\varphi \left(y\right)$ near $y\sim 0$ to account for the divergence of ${\varphi }^{\prime \prime }\left(0\right)$. However, this is not considered here since the dominant term contributing to the integral (in $[0,x_1]$) is $\varphi \left(0\right)/y^{1/2}$. All the integrals are readily evaluated using low order (say 6 or 8) Gauss quadrature. However, the asymptotic term $[y(2/3)\varphi \left(0\right)/y^{3/2}]$ should be subtracted within the integral in Eq.(18) (before applying Gauss quadrature formula) and its contribution $\left[(4/3)\varphi \left(0\right)x_1^{1/2}\right]$ added separately. Similar technique may be used in obtaining $S^{*}\left(x_1\right)$ from Eq.(17) for better accuracy. The approximation leading to Eq.(13) is the same as in obtaining Simpson’s integration formula and hence the local truncation error is proportional to ${\Delta }_N^4$. The collocation rule in Nyström’s method amounts to enforcing Eq.(13) 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}}^{*}$ of order $(N+1)$ for the vector $\overrightarrow{\varphi }$. The coefficient matrix is $\mathbf{M}=\mathbf{I}+{\displaystyle \frac{\mathcal{A}}{\mathbf{2}}}\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}}^{*}$ and $\mathbf{x}$ are $S_j^{*}$ and $x_j$, respectively. Elements in a particular column of the lower triangular part of $\mathbf{M}$ are the same. The upper triangular part can be written as $\mathbf{I}\mathbf{x}\mathbf{U}$ with $\mathbf{U}$ sharing the same property.

The number of Newton’s iterations $N_{it}$ needed for convergence to a specified accuracy is indeed depended on the starting guess solution. For illustration, Figure 1 shows $N_{it}$ needed versus $T{Z}^{-4/3}$ to reduce the Euclidean norm between successive iterates to ${10}^{-10}$. The starting guess solution for each temperature corresponds to uniform electron density and $N=20$ are the number of sub-intervals used. The scaled cell radius in this case is $R{Z}^{1/3}=4.339A^o$, however, similar results are obtained for expanded as well as compressed configurations. Even for very low temperatures $T{Z}^{-4/3}\sim {10}^{-4}$ eV, Newton’s method converges in about 10 iterations to the specified accuracy. The inset in this figure is the converged profile $\varphi \left(x\right)$ for reduced temperature $T{Z}^{-4/3}=1.122\times {10}^{-4}$ eV, which justify piece-wise quadratic polynomials used in the method. It is also found that the spatial profile of $\varphi \left(x\right)$ converges with about 10 sub-intervals N.

3.1. Thermodynamic Quantities

All the properties described earlier (see thermodynamic-properties) is computed using the converged solution vector $\varphi$. For instance, pressure is obtained directly from Eq.(3) using the converged value of ${\varphi }_N$. Next, the kinetic energy is obtained as

$$E_{kin}={\displaystyle \frac{\mathcal{A}R}{e^2}}{\left(k_{B}T\right)}^2\left\{\underset{0}{\overset{x_1}{\int }}{y}^2{I}_{3/2}\right[{\displaystyle \frac{\varphi \left(y\right)}{y}}]dy+\sum _{j=2★}^N\underset{x_{j-1}}{\overset{x_{j+1}}{\int }}{y}^2{I}_{3/2}[{\displaystyle \frac{\varphi \left(y\right)}{y}}\left]dy\right\}`$$

(19)

Note that the second sum is only over even values of j. The contribution of the asymptotic term in $y^2{I}_{3/2}(\varphi /y)\sim (2/5)\varphi {\left(0\right)}^{5/2}/y^{1/2}$ to the first integral should be accounted as explained earlier in the case of ${\varphi }^{\prime }\left(0\right)$. Further, $E_{pot}$ follows from the virial theorem. Alternatively, the thermal energy $E_{th,n}$ at temperature $T_n$ is obtained directly via a recursion relation

$$E_{th,n}=(T_n/T_{n+1}{)}^{7/4}\{E_{th,n+1}-(3/4)g_{n+1}\}+(3/4)g_n-(9/16)T_n^{7/4}\underset{T_n}{\overset{T_{n+1}}{\int }}{t}^{-11/4}g\left(t\right)dt.$$

(20)

This follows directly from Eq.(8), and the dependence of $E_{th,n}$ and $g_n$ on cell volume $\mathcal{V}$ is suppressed. For sufficiently large n, the starting value of thermal energy $E_{th,n+1}=(3/2)Z{k}_B{T}_{n+1}$. An interpolation function of $g(\mathcal{V},t)$ is first constructed from its tabulation versus T and then used in this relation. The last integral in Eq.(20) is computed with three-point Simpson’s rule. The slope ${\varphi }^{\prime }\left(0\right)$ is given by Eq.(18) (also see analytical-representation).

4. Analytical Representations

Using the results of extensive calculations, covering the temperature range ${10}^{-4}\le T{Z}^{-4/3}\le {10}^3$ (eV) and cell radius range $0.93\le R{Z}^{1/3}\le 20.1$ (A^o), analytical representations for thermal energy $E_{th}(R,T)$, thermal pressure $P_{th}(R,T)$, thermal ionization $Z^{*}(R,T)$, Fermi energy ${\mathcal{E}}_F(R,T)$ and initial slope ${\varphi }^{\prime }\left(0\right)$ of TF-function are derived. As explained earlier, while ionization provides Fermi energy, the initial slope gives entropy and free energy.

1. Thermal Energy

Thermal energy is represented in the form

$$E_{th}(R,T)=Z^{7/3}\{w_1{\overline{T}}^2+w_2{\overline{T}}^3\}\{1+w_3{\overline{T}}^{w_4}+{\overline{T}}^2{\}}^{-1}.$$

(21)

The factor $Z^{7/3}$ is the scale parameter for energy in TF-model and $\overline{T}=T/Z^{4/3}$ (eV) is the scaled temperature. The parameters $w_n$ (for $n=1,3,4$) are functions of the scaled cell radius $\overline{R}=R{Z}^{1/3}$ (A^o). The expression for $E_{th}$ is proportional to the asymptotic forms $\sim {\overline{T}}^2$ and $\sim \overline{T}$ in the limits $\overline{T}\to 0$ and $\overline{T}\to \infty$, respectively. Therefore, it is evident that $w_2=3/2$. The parameter $w_1$ is already determined by Gilvarry [21,22] using solutions of temperature-perturbed TF equation. It’s fit as a function of $x_b=\overline{R}/{\mu }^{*}$ is

$$w_1=P_0\left(\overline{R}\right){\displaystyle \frac{15}{2}}(\sigma +2\tau +3\omega )\xi {\left[Z^{4/3}\right]}^2.$$

Here, $P_0\left(\overline{R}\right)$ is the pressure at zero-temperature and $\sigma ={\chi }_b/{\varphi }_b$ , $\omega ={\chi }_i^{\prime }/\left(x_b^{1/2}{\varphi }_b^{5/2}\right)$ and $\tau ={(x_b/{\varphi }_b)}^2$. The factor ${\left[Z^{4/3}\right]}^2$ arises because $w_1{\overline{T}}^2$ should be proportional to $T^2$. The dimensionless parameter $x_b=R{Z}^{1/3}/{\mu }^{*}$ is proportional scaled cell-radius as ${\mu }^{*}=a_0{(9{\pi }^2/128)}^{1/3}$ ($a_0$ the Bohr radius) and $\xi =({\pi }^2/8){({\mu }^{*}{Z}^{-4/3}/e^2)}^2$ are constants. The zero-temperature pressure $P_0\left(\overline{R}\right)$, the boundary value of zero-temperature TF-function ${\varphi }_b$ and the parameters ${\chi }_b$ and ${\chi }_i^{\prime }$ are fitted functions of $x_b$ given by

$$\begin{array}{ccc}\hfill P_0& =& \left\{Z^{10/3}{e}^2\right\}{\left\{10\pi {\mu }^{*4}\right\}}^{-1}{\{0.48075x_b^2+0.06934x_b^3+0.0033704x_b^4+0.0097x_b^{7/2}\}}^{-5/2}\hfill \\ \hfill {\varphi }_b& =& {\{0.48075x_b^{2/2}+0.06934x_b^{4/2}+0.0097x_b^{5/2}+0.0033704x_b^{6/2}\}}^{-1}\hfill \\ \hfill {\chi }_b& =& -0.3205x_b^3-0.04021x_b^4-0.002519x_b^5\hfill \\ \hfill {\chi }_i^{\prime }& =& {\{-0.005805x_b^{-0.228}-0.3755x_b^{-1}-3.12x_b^{-2}\}}^{-1}\hfill \end{array}$$

The remaining parameters $w_3\left(\overline{R}\right)$ and $w_4\left(\overline{R}\right)$ are fitted to the rational expressions

$$\begin{array}{ccc}\hfill w_3& =& \{-9.46531+2.1505\overline{R}-0.33598{\overline{R}}^2+0.00357102{\overline{R}}^3-4.5534*{10}^{-5}{\overline{R}}^4+\hfill \\ & +& 1.9721*{10}^{-7}{\overline{R}}^5\}{\{1-0.829696\overline{R}+0.00159347{\overline{R}}^2\}}^{-1}\hfill \\ \hfill w_4& =& \{0.952738+0.01242\overline{R}+0.00373893{\overline{R}}^2\}{\{1+0.0171359\overline{R}+0.00734102{\overline{R}}^2\}}^{-1}\hfill \end{array}$$

As application of these formulas, plots of scaled energy $E_{th}/Z^{7/3}$ (eV/atom) versus scaled temperature $\overline{T}=T{Z}^{-4/3}$ (eV) are shown in Figure 2 for three values of scaled cell radius $\overline{R}=R{Z}^{1/3}\left(A^o\right)$. Here, the symbols are numerical values and lines represent analytical representations. It is important to note that all the three curves merge together at high temperatures thereby indicating approach to ideal gas limit.

2. Thermal Pressure

The analytical representation for thermal pressure is chosen as

$$P_{th}(R,T)=Z^{10/3}\{p_1{\overline{T}}^2+p_2{\overline{T}}^3\}\{1+p_3{\overline{T}}^{p_4}+{\overline{T}}^2{\}}^{-1}.$$

(22)

Here the factor $Z^{10/3}$ is the scale parameter for pressure. Now, all the parameters $p_n$ ($1\le n\le 4$) are functions of the scaled cell radius $\overline{R}=R{Z}^{1/3}$. Since the expression for $P_{th}$ is proportional to the asymptotic form $\sim \overline{T}$ in the limit $\overline{T}\to \infty$, the parameter $p_2={\left(Z\mathcal{V}\right)}^{-1}$. The parameter $p_1$ (already determined by Gilvarry [21,22] using the temperature-perturbed TF equation) is given by

$$p_1=P_0\left(\overline{R}\right){\displaystyle \frac{5}{2}}(\sigma +2\tau )\xi {\left[Z^{4/3}\right]}^2$$

The factor ${\left[Z^{4/3}\right]}^2$ arises in this expression because $p_1{\overline{T}}^2$ should be proportional to $T^2$. Other quantities $\sigma$, $\tau$ and $\xi$ are already given above. The parameters $p_3$ and $p_4$ are represented as

$$\begin{array}{ccc}\hfill p_3& =& \{-12.795+3.98404\overline{R}-0.659487{\overline{R}}^2+0.000101473{\overline{R}}^3-0.0000148037{\overline{R}}^4\}\hfill \\ & \times & {\{1-0.828017\overline{R}-0.00295839{\overline{R}}^2\}}^{-1}\hfill \\ \hfill p_4& =& \{0.735021-0.0368263\overline{R}+0.0144916{\overline{R}}^2-0.0000886233{\overline{R}}^3+5.74717*{10}^{-7}{\overline{R}}^4\}\hfill \\ & \times & {\{1-0.172122\overline{R}+0.0251956{\overline{R}}^2\}}^{-1}.\hfill \end{array}$$

Plots of $P_{th}/Z^{10/3}$ (GPa) versus $\overline{T}=T/Z^{4/3}$ (eV) are shown in Figure 3 for same values of scaled cell radius $\overline{R}=R{Z}^{1/3}\left(A^o\right)$. Here also the symbols denote numerical values while the lines depict analytical representations. The linear variation at high temperature in all the three curves indicates approach to ideal gas limit. Unlike the fits available in the literature [11,12], the expressions in Eq.(21) and 22 have the correct dependence on cell radius $\overline{R}$, and are applicable in compressed as well as expanded volume regions (see applications).

4. Thermal Ionization

As the numerical procedure gives $\varphi \left(1\right)={\mathcal{E}}_F/k_{B}T$, thermal ionization is directly given in terms of $I_{1/2}({\mathcal{E}}_F/k_{B}T)$ (see thomas-fermi-model). A numerical fit to $Z^{*}$ as a function of density and temperature is given by More [23]. However, it is found that the relative error in this formula is as much as 30 percent for low temperature and expanded volume states. The data generated now is expressed in rational function representation as

$$Z^{*}(R,\overline{T})=Z\{z_0+z_1{\overline{T}}^{z_4}+z_2{\overline{T}}^{2.15}\}\{1.+z_3\overline{T}+z_2{\overline{T}}^{2.15}{\}}^{-1}.$$

(23)

This form approaches approaches Z as $\overline{T}\left(eV\right)\to \infty$, and the parameter $z_0$ provides ionization at zero-temperature. Theoretically, the value of $z_0$ in the limit $\overline{R}\to 0$ should approach 1; the numerical value obtained $\sim 1.05635$ is close to this theoretical limit. Plots of $Z^{*}$ versus $\overline{T}$ are shown in Figure 4 for same values (0.9348, 4.339 and 20.14 $A^o$) of scaled cell radius $\overline{R}=R{Z}^{1/3}\left(A^o\right)$ considered earlier. Again, symbols are numerical values and lines analytical representations. The present expression has about $\sim 4$ percent accuracy in the entire plane.

Dependence of all the five parameters $z_n$$(0\le n\le 4)$ on $\overline{R}$ is obtained as

$$\begin{array}{ccc}\hfill z_0& =& \{1.05635+0.0348859\overline{R}-0.00050309{\overline{R}}^2\}\hfill \\ & \times & {\{1+0.639018\overline{R}+0.158541{\overline{R}}^2+0.0143746{\overline{R}}^3\}}^{-1}\hfill \\ \hfill z_1& =& \{-0.00685217+0.0666106\overline{R}+0.00659339{\overline{R}}^2\}\hfill \\ & \times & {\{1-0.182044\overline{R}+0.0162883{\overline{R}}^2-0.000677451{\overline{R}}^3+0.0000117727{\overline{R}}^4\}}^{-1}\hfill \\ \hfill z2& =& \{-0.0932445+0.150902\overline{R}-0.0691882{\overline{R}}^2+0.0129664{\overline{R}}^3-0.000409478{\overline{R}}^4\}\hfill \\ & \times & {\{1-0.151435\overline{R}+0.0146767{\overline{R}}^2-0.000672239{\overline{R}}^3+0.0000111289{\overline{R}}^4\}}^{-1}\hfill \\ \hfill z3& =& \{-0.560091+0.955685\overline{R}-0.365586{\overline{R}}^2+0.0766116{\overline{R}}^3-0.00267292{\overline{R}}^4\}\hfill \\ & \times & {\{1-0.105577\overline{R}+0.0144308{\overline{R}}^2-0.000610433{\overline{R}}^3+5.68036*{10}^{-}6{\overline{R}}^4\}}^{-1}\hfill \\ \hfill z4& =& \{0.709056+1.05963\overline{R}-0.0861352{\overline{R}}^2+0.00187143{\overline{R}}^3\}\hfill \\ & \times & {\{1+0.756378\overline{R}-0.064099{\overline{R}}^2+0.00143063{\overline{R}}^3\}}^{-1}\hfill \end{array}$$

3. Fermi Energy

Although the solution procedure of TF model directly gives ${\mathcal{E}}_F/k_{B}T$, it can also be obtained in terms of ionization using the inverse of $I_{1/2}$ integral (see thomas-fermi-model). Since very accurate representation of this inverse function is available [17], it is unnecessary to represent ${\mathcal{E}}_F$ directly. However, for the sake of completeness, scaled values ${\mathcal{E}}_F{Z}^{-4/3}$ versus scaled temperature $\overline{T}=T/Z^{4/3}$ (eV) are shown in Figure 5 for the three values (0.9348, 4.339 and 20.14 $A^o$) of scaled radius $\overline{R}=R{Z}^{1/3}$$\left(A^o\right)$.

4. Initial Slope

The initial slope of the Thomas-Fermi function ${\varphi }^{\prime }\left(0\right)$ is needed in calculating entropy (and free energy) as it is defined by $TS=-P\mathcal{V}+(7/3)E-Z{k}_{B}T{\varphi }^{\prime }\left(0\right)$. To this end, the parameter $\psi =k_B\overline{T}{\varphi }^{\prime }\left(0\right)$ is introduced so that entropy is expressed in terms of thermal quantities as $TS=-P_{th}\mathcal{V}+(7/3)E_{th}-Z{\psi }_{th}$. The representation chosen $\psi$ is

$$\psi (\overline{R},\overline{T})={\psi }_0\left(\overline{R}\right){\psi }_{fit}(\overline{R},\overline{T})+{\psi }_{uni}(\overline{R},\overline{T}).$$

(24)

Here, ${\psi }_{uni}={\mathcal{E}}_F/Z^{4/3}-(3/2)e^2/\overline{R}$ is the limiting form of $\psi$ corresponding to uniform distribution of electrons. This follows readily from ${\varphi }_{uni}\left(x\right)=x({\mathcal{E}}_F/k_{B}T)+\varphi \left(0\right)(1+x^3/2-3x/2)$. Further, ${\psi }_0\left(\overline{R}\right)$ is the limiting value of $\psi$ as $\overline{T}\to 0$. This separation is useful as ${\mathcal{E}}_F/Z^{4/3}\to -log[8\pi /3]k_B\overline{T}-(3/2)k_B\overline{T}log[2\pi m{k}_B\overline{T}{\overline{R}}^2/h^2]$ as $\overline{T}\to \infty$. The representation for ${\psi }_{fit}$, where ${\mathcal{E}}_F$ is computed using ionization $Z^{*}(\overline{R},\overline{T})$, is given by

$${\psi }_{fit}=a_0\{1+exp\left[a_1\{{log}_{10}\left(\overline{T}\right)-a_2\}\right]{\}}^{-0.45}.$$

(25)

where $\overline{T}$ and $\overline{R}$ are in eV and A^o, respectively. The parameters $a_n$ ($0\le n\le 2$) and the zero-temperature slope ${\psi }_0\left(\overline{R}\right)$ are found to be

$$\begin{array}{ccc}\hfill a_0& =& \{0.559395+0.683896\overline{R}\}{\{1+0.675501\overline{R}\}}^{-1}\hfill \\ \hfill a_1& =& \{2.06764+1.48356\overline{R}-0.000678249{\overline{R}}^2\}{\{1+1.15387\overline{R}\}}^{-1}\hfill \\ \hfill a_2& =& \{6.00953+1.17532\overline{R}-0.0142208{\overline{R}}^2\}{\{1+1.12646\overline{R}\}}^{-1}\hfill \\ \hfill {\psi }_0& =& \{2.19714-107.115\overline{R}-98.7829{\overline{R}}^2-104.134{\overline{R}}^3-40.1999{\overline{R}}^4\}\hfill \\ & *& {\{1+1.52848\overline{R}+2.17814{\overline{R}}^2+2.11488{\overline{R}}^3+0.824127{\overline{R}}^4\}}^{-1}\hfill \end{array}$$

The variation of ${\psi }_{th}/\overline{T}=(\psi -\psi \left(0\right))/T{Z}^{-4/3}$ versus $T{Z}^{-4/3}$ is shown in Figure 6 for three different values (0.9348, 4.339 and 20.14 $A^o$) of $\overline{R}=R{Z}^{1/3}$. The chosen functions represent the data accurately. For the sake of completeness, the dependence of ${\psi }_0$ versus $\overline{R}$ is shown in the inset in this figure.

The definition of $\psi$ and Eq.(8) shows that ${\psi }_0$ is the same as $(e^2/{\mu }^{*}){\varphi }_{0in}^{\prime }$ . The fit for ${\varphi }_{0in}^{\prime }$ (initial slope of zero-temperature TF-function) provided by Gilvarry [21,22] is

$${\varphi }_{0in}^{\prime }={\varphi }_{0in\infty }^{\prime }+{\{\sum _{j=1}^3{b}_j{x}_b^{j+1}+b_4{x}_b^{4.5058}+\sum _{j=5}^7{b}_j{x}_b^j+b_8{x}_b^{7.772}\}}^{-1},$$

with $b_1=0.48075,b_2=0.41602,b_3=0.1052,b_4=0.032589,b_5=.010578,b_6=0.0027271,b_7=0.0000415,b_8=0.000001786$. Note that ${\varphi }_{0in\infty }^{\prime }=-1.58807103$ is the slope for isolated atom. This yields $E_{\infty }=(3/7)Z^{7/3}(e^2/{\mu }^{*}){\varphi }_{0in\infty }^{\prime }=-0.7687451Z^{7/3}(e^2/a_0)$ as the energy of the isolated atom at zero-temperature. Numerical results show that the two representations for ${\psi }_0$ agree well.

5. Applications

A three-component description of the EOS [24,25] employs the zero-temperature isotherms $(P_c,E_c)$, the ionic components of thermal energy and pressure $(P_{it},E_{it})$ and the electron components $(P_{th},E_{th})$. Usually, the electron components are obtained as interpolated data from large tables [24]. For applications that do not involve very high temperatures, empirical fits [11,12] of the data published by Later [4] are also employed. These fit are interpolation functions between low-temperature electron properties (as in Fermi-gas model) and the ideal gas model. The representations provided in this work (see analytical-representation) offers an alternate procedure which provides results close to that obtainable from data tables.

5.1. Electron Hugoniot

The principal Hugoniot of only electrons within the Fermi-gas model is well known [26]. Considering a point nuclear charge $Ze$ at the cell center, together with the analytical representations of $P_0+P_{th}$ and $E_0+E_{th}$, the electron-Hugoniot is easily computed. In the results shown in Figure 7, pressure ratio $P/P_0$ is shown versus the compression factor $V_0/V$. Curve-1 corresponding to TF model shows that the maximum compression obtainable is $\sim 5$. This figure also contains the results obtained using the empirical fits [11,12] mentioned above. The maximum compression of $\sim 4$ is found with the fits because of the absence of any binding effect of electrons in the cell. The inset in this figure shows reduced specific heat $(C_V/k_B)$ (per electron) for cell radius $\overline{R}=4.339A^o$. The peak found in $C_V$ around $\overline{T}\sim 3.5\times {10}^4$ K is due to electron binding to the nucleus.

5.2. Hugoniot of Cu

For application to Cu, Vinet’s zero-temperature isotherm [27,28] given by $P_v=3B_0{(V/V_{c0})}^{-2}(1-(V/Vc0))exp[(1-(V/V_{c0}))(B_0^{\prime }-1)3/2]$ is employed. The corresponding analytical expression for energy $E_v$ follows from the relation $d{E}_v/dV=-P_v$ and $E_v\left(V_{c0}\right)=E_{coh}$. Here $B_0$, $B_0^{\prime }$ and $E_{coh}$ are, respectively, the bulk modulus, its pressure derivative and cohesive energy (per gram) at volume $V_{c0}$. This volume is determined such that total pressure (including those of ions and electrons) is ${10}^{-4}$ GPa at ambient temperature $T_0$ (300 K) and volume $V_0$. These formulas are suitably corrected using QSM [7,8] at very high pressures. Pressure within QSM is analytically expressed as $P_q\left(V\right)=({\hslash }^2/2m)[(2/5){\left(3{\pi }^2\right)}^{2/3}{N}_s^{5/3}-(2/a_B)(13/36){(3/\pi )}^{1/3}{N}_s^{4/3}]$ (ℏ is reduced Planck’s constant, m electron-mass and $a_B$ the Bohr radius). The quantity $N_s\left(\mathcal{V}\right)$ (electron density at the atomic cell surface) is given by $N_s=(Z/V)exp[-\alpha (R/a_B)-\beta {(R/a_B)}^2]$ where $\alpha =0.1935Z^{[0.495-0.039\left({log}_{10}Z\right)]}$ and $\beta =0.068+0.078\left({log}_{10}Z\right)-0.086{\left({log}_{10}Z\right)}^2$ (R denotes cell radius). Energy in QSM also follows from the relation $d{E}_q/dV=-P_q$. An interpolation procedure [29] to smoothly go over from Vinet’s model to QSM uses an adjustable volume $V_m$ such that the corrected pressure $P_c=P_v$ and $E_c=E_v$ for $V\ge V_M$. These for smaller volumes are obtained by imposing continuity of $P_c\left(V\right)$ and its first two derivatives at $V_M$. Note that the definition ensures that $E_c\left(V\right)$ is continuous at $V_M$ [30].

Ionic thermal contributions are obtained using Johnson’s model [31] wherein the thermal ionic energy (per gram) and pressure are expressed as $E_{it}(V,T)=E_D+N{k}_{B}T(E_0+E_{\psi })$ and $P_{it}(V,T)=P_D+(1/V)(2{\Gamma }_s-2/3)N{k}_{B}T(E_0+E_{\psi })$. Here, Debye model $(P_D,E_D)$ provides the EOS in solid phase. The components $E_0$ and $E_{\psi }$ (non-zero for $T\ge T_M$, the melting temperature) accounts for changes in energy due to melting, heat of fusion and subsequent approach to ideal gas behavior [25,31]. A three-term fit [32] to the Grüneisen parameters ${\Gamma }_s\left(V\right)$ is employed in the present application.

The principal Hugoniot and zero-temperature isotherm are compared with available experimental (also theoretical) data [33] in Figure 8. The Hugoniot agrees very well over an extended range of pressures up to $20.4$ TPa. The zero-temperature isotherm is compared over a larger range up to 100 TPa pressure in the insert. Improvements to $P_c$ due to the corrections using QSM (with $V_M=V_0/2$) are also shown here. To check the ionic thermal model, the constant-pressure specific heat $C_P$ in the solid phase at $V_0$ is compared with reported results [33] in Figure 9. This figure also depicts $C_V$ by the dashed line. Results agree well in spite of the simplicity of the ionic model. The sharp increase in $C_P$ and $C_V$ at $T_M=1360$ K occurs due to the addition of heat of fusion. The room-temperature (300 K) isotherm shown in the inset agrees excellently.

6. Summary

An efficient numerical algorithm based on Newton’s method to solve the TF equation is developed in this paper. The basic idea in the method, which rewrites the TF equation as a non-linear Fredholm integral equation, circumvents some of the pitfalls in earlier work using FDM. Results presented show that Newton’s iterations converge within few iterations in the whole temperature-density plane. Furthermore, the algorithm needs small number of sub-intervals for discrete approximations. The use of Brachman’s equation for direct calculation of electron-thermal energy is elaborated. Analytical representations of thermodynamic quantities are obtained from extensive tables of data obtained with the present method. These include thermal energy, pressure, ionization, Fermi energy and initial slope of TF-function (which is needed for entropy and free energy). Accuracy of the representations is ascertained via computation of the electron Hugoniot and principal Hugoniot of Cu and comparison with experimental (or theoretical) data. The present paper significantly extends to all temperatures the earlier work (restricted to low temperatures) of Gilvarry on the TF model.