The valence electrons of a metal experience an approximately constant potential throughout the crystal, and the potential energy term can be set to zero, $V(x) = 0$. The Schrödinger equation for this case is,

$$- \frac{\hbar^2}{2m}\nabla^2 \psi(\vec{r}) = E \psi(\vec{r}).$$

To avoid the issues associated with the surface of the crystal we apply periodic boundary conditions i.e., $\psi(0,y,z) = \psi(L,y,z),\,\psi(x,0,z) = \psi(x,L,z),\,\psi(x,y,0) = \psi(x,y,L)$. The normalized solutions to the Schrödinger equation are,

$$\psi(\vec{r}) = \frac{1}{\sqrt{L^3}}e^{i\vec{k}\cdot\vec{r}},$$

where the periodic boundary conditions require that $\vec{k} = \left[\frac{2\pi n_x}{L},\frac{2\pi n_y}{L},\frac{2\pi n_z}{L}\right]$ and the $n$'s are positive or negative integers, i.e., $n_x,n_y,n_z = \cdots -2,\,-1,\,0,\,1,\,2,\,\cdots$.

$\psi$		$n_x=$1
	$x$

This figure shows the wave function $e^{i2\pi n_xx}$ in one dimension. Notice that $\psi^*\psi$ is a constant for every value of $n_x$. It is equally probable to find the particle everywhere.

The $\vec{k}$ vectors that satisfy the boundary conditions form a cubic grid, and each $\vec{k}$ vector is associated with a reciprocal space volume of $\left(\frac{2\pi}{L}\right)^3$. Thus, the density of electron states in reciprocal space is $D(\vec{k}) = \frac{2}{(2\pi)^3}$. There are twice as many electron states as $\vec{k}$-states because of spin.

Multiplying the real space volume $L^3$ times $D(\vec{k})$ times a volume of $k$-space tells you how many allowed electron states are in that $k$-space volume. Sometimes we only want to know how many electron states there are with values between $|k_1|$ and $|k_2|$ without regard to which direction the $\vec{k}$-vectors are pointing. That would be $D(\vec{k})$ times the volume between 2 spherical shells, one with radius $|k_1|$ and the other with radius $|k_2|$. The number of $\vec{k}$-points between $k$ and $k+dk$ is then

$$L^3D(\vec{k})4\pi k^2dk= N =\frac{L^3k^2}{\pi^2}dk = L^3D(|\vec{k}|) dk.$$

Here we have defined

$$D(|\vec{k}|) =\frac{k^2}{\pi^2}\quad\text{m}^{-2}$$

as the density of electrons states in $|\vec{k}|$ per unit volume. It grows like $k^2$ because the number of $\vec{k}$-points with the same magnitude $k$ lie on a spherical shell and the area of the spherical shell grows like $k^2$ as $k$ gets larger.

The energy of the solutions can be found by substituting them into the Schrödinger equation, which yields,

$$E = \frac{\hbar^2 k^2}{2m}\quad\text{J}.$$

This is called the energy-dispersion relation and is just the kinetic energy since the potential energy was set to zero. To determine how many states there are in an energy interval between $E$ and $E+dE$ we can equate,

$$D(E)dE = D(|\vec{k}|) dk\qquad\rightarrow\qquad D(E) = D(|\vec{k}|) \frac{dk}{dE}.$$

From the dispersion relation we obtain, $dE = \frac{\hbar^2 k}{m}dk$, and the energy density of electron states can be calculated to be,

$$D(E) = \frac{\sqrt{2}m^{3/2}}{\pi^2\hbar^3}\sqrt{E}\quad\text{J}^{-1}\,\text{m}^{-3}.$$

The electron density of states grows like $\sqrt{E}$, so there are more electron states at higher energies. The density of states was calculated for periodic boundary conditions, but it turns out to be the same if fixed boundary conditions are used. At low temperatures, only the states with the lowest energy levels will be occupied. The electron density $n$ is then equal to the sum of the occupied states,

$$n =\int\limits_{-\infty}^{E_F}D(E)dE\quad\text{m}^{-3}.$$

Here $E_F$ is the Fermi energy. In a metal, the Fermi energy is the energy of the highest occupied electron state.

For free electrons,

$$n = \frac{2\sqrt{2}m^{3/2}}{3\pi^2\hbar^3}E_F^{3/2}\quad\text{ m}^{-3}\qquad\rightarrow\qquad E_F = \frac{\hbar^2}{2m}\left(3\pi^2n\right)^{2/3}\quad\text{ J}.$$

An alternate derivation of this result is to say that the states in reciprocal space are occupied up to the Fermi wave vector $k_F$.

$$n = \int\limits_{-\infty}^{k_F}D(|\vec{k}|) dk.$$

The volume of reciprocal space that is occupied is a sphere with volume $4\pi k_F^3/3$. The volume of a single $k$ state is $(2\pi)^{-3}$ per unit volume and each $k$ state can hold two electrons (spin) so,

$$ n = \frac{2}{(2\pi)^3}\frac{4\pi k_F^3}{3}\quad\text{ m}^{-3}\qquad \rightarrow \qquad k_F = (3\pi^2 n)^{1/3}\quad\text{ m}^{-1},$$ $$E_F = \frac{\hbar^2k_F^2}{2m}\quad\text{J.}$$

The density of states takes a simple form when it is written in terms of the Fermi energy,

$$D(E) = \frac{3}{2}\frac{n\sqrt{E}}{E_F^{3/2}}\quad\text{J}^{-1}\text{ m}^{-3}.$$

The contribution of the electrons to the thermodynamic properties depends mostly on the states near the Fermi energy so it is useful to calculate the density of states at the Fermi energy and the derivative of the density of states at the Fermi energy, $$D(E_F) = \frac{3}{2}\frac{n}{E_F}\quad\text{ J}^{-1}\text{ m}^{-3},\qquad \frac{dD(E)}{dE}\Big{|}_{E_F} =D'(E_F) =\frac{3}{4}\frac{n}{E_F^2}\quad\text{ J}^{-2}\text{ m}^{-3}.$$

At finite temperatures, the occupation of the electron states is given by the Fermi function,

$$f(E) = \frac{1}{\exp{\left(\frac{E-\mu}{k_BT}\right)}+1}.$$

Here $\mu$ is the chemical potential. At $T=0$, $\mu= E_F$. The electron density at finite temperatures is then,

$$n =\int\limits_{-\infty}^{\infty}D(E)f(E)dE = \int\limits_{-\infty}^{\infty}\frac{D(E)}{\exp{\left(\frac{E-\mu}{k_BT}\right)}+1}dE.$$

This equation can't be solved analytically for $\mu$ but it is possible to determine $\mu$ numerically by guessing a chemical potential, performing the integral, and comparing the resulting electron density with the desired electron density. This process is repeated until the correct chemical potential is found. A suitable initial guess is $\mu = E_F$ because the temperature dependence of a typical metal is weak.

The Sommerfeld expansion can also be used to find the temperature dependence of the chemical potential at low temperatures. The Sommerfeld result is,

$$\mu(T) \approx E_F - \frac{\pi^2}{12}\frac{(k_BT)^2}{E_F}\quad\text{ J}.$$

The internal energy density $u$ is the energy $E$, times the number of electron states at that energy $D(E)$, times the probability that a state at that energy is occupied $f(E)$, summed over all energies,

$$u =\int\limits_{-\infty}^{\infty}ED(E)f(E)dE \approx \frac{3}{5}nE_F + \frac{\pi^2}{4}\frac{n}{E_F}(k_BT)^2\quad\text{ J m}^{-3}.$$

Once the internal energy density is known, other thermodynamic quantities can be calculated using the laws of thermodynamics. The specific heat at constant volume is the derivative of the internal energy density with respect to temperature,

$$c_v = \left(\frac{\partial u}{\partial T}\right)_V \approx \frac{\pi^2}{2}\frac{n}{E_F}k_B^2T\quad\text{ J K}^{-1}\text{ m}^{-3}.$$

The contribution of the electrons to the specific heat is linear in the temperature.

The entropy density $s$ can be calculated from the specific heat,

$$s =\int\frac{c_v}{T}dT \approx \frac{\pi^2}{2}\frac{n}{E_F}k_B^2T\quad\text{ J K}^{-1}\text{ m}^{-3}.$$

The Helmholtz free energy density is,

$$f = u -Ts \approx \frac{3}{5}nE_F - \frac{\pi^2}{4}\frac{n}{E_F}(k_BT)^2\quad\text{ J m}^{-3}.$$

The kinetic energy of the electrons will cause a pressure if they strike a boundary. This pressure is $P= -\left(\frac{\partial F}{\partial V}\right)_{N,T}$ where $F$ is the total Helmholtz free energy. Using the relations, $F=Vf$, $n = N/V$, and $E_F = \frac{\hbar^2}{2m}\left(\frac{3\pi^2 N}{V}\right)^{2/3}$, where $N$ is the total number of electrons and $V$ is the volume, the total free energy can be written,

$$F \approx \frac{3}{5}N\frac{\hbar^2(3\pi^2)^{2/3}}{2m}\left(\frac{N}{V}\right)^{2/3} - \frac{\pi^2N}{4}\frac{2m}{\hbar^2(3\pi^2)^{2/3}}\left(\frac{V}{N}\right)^{2/3}(k_BT)^2\quad\text{J}.$$ $$P = -\left(\frac{\partial F}{\partial V}\right)_{N,T} \approx \frac{2}{5}nE_F + \frac{\pi^2}{6}\frac{n}{E_F}(k_BT)^2\quad\text{ N m}^{-2}.$$

The bulk modulus $B$ tells us how much pressure we need to apply to change the volume of the system, $B = -V\left(\frac{\partial P}{\partial V}\right)_{N,T}$. Expressing $P$ in terms of $N$ and $V$ we have,

$$P \approx \frac{2}{5}\frac{\hbar^2(3\pi^2)^{2/3}}{2m}\left(\frac{ N}{V}\right)^{5/3} + \frac{\pi^2}{6}\frac{2m}{\hbar^2(3\pi^2)^{2/3}}\left(\frac{N}{V}\right)^{1/3}(k_BT)^2\quad\text{ N m}^{-2}.$$ $$B = -V\left(\frac{\partial P}{\partial V}\right)_{N,T} \approx \frac{2}{3}nE_F + \frac{\pi^2}{18}\frac{n}{E_F}(k_BT)^2\quad\text{ N m}^{-2}.$$

The bulk modulus and the specific heat are two experimentally accessible quantities that are often used to evaluate how well a metal can be described by the free-electron model.

The enthalpy density $h$ is,

$$h = u + P \approx nE_F + \frac{5\pi^2}{12}\frac{n}{E_F}(k_BT)^2 \quad\text{ J m}^{-3}.$$

The Gibbs energy density $g = u+P-Ts = h-Ts$ is,

$$g = h-Ts \approx nE_F - \frac{\pi^2}{12}\frac{n}{E_F}(k_BT)^2\quad\text{ J m}^{-3}.$$

Formulas for the free-electron model in one-, two-, and three-dimensions are summarized in the table below.

$\,\,$Single particle solutions $\,\,$to Schrödinger equation	$$\psi(\vec{r}) = \frac{1}{\sqrt{L}}e^{ikx}\quad \frac{1}{\sqrt{\text{m}}}$$	$$\psi(\vec{r}) = \frac{1}{\sqrt{L^2}}e^{i\vec{k}\cdot\vec{r}}\quad \frac{1}{\text{m}}$$	$$\psi(\vec{r}) = \frac{1}{\sqrt{L^3}}e^{i\vec{k}\cdot\vec{r}}\quad \frac{1}{\text{m}^{3/2}}$$
$\,\,$Allowed k values	$$k_x = \frac{2\pi n_x}{L}\quad\frac{1}{\text{m}}\\n_x= \cdots -2,\,-1,\,0,\,1,\,2,\,\cdots$$	$$\vec{k} = \left[\frac{2\pi n_x}{L},\frac{2\pi n_y}{L}\right]\quad\frac{1}{\text{m}}\\n_x,n_y = \cdots -2,\,-1,\,0,\,1,\,2,\,\cdots$$	$$\vec{k} = \left[\frac{2\pi n_x}{L},\frac{2\pi n_y}{L},\frac{2\pi n_z}{L}\right]\quad\frac{1}{\text{m}}\\n_x,n_y,n_z = \cdots -2,\,-1,\,0,\,1,\,2,\,\cdots$$
$\,\,$Density of electron states $\,\,$in reciprocal space	$$D(\vec{k}) = \frac{2}{2\pi}$$	$$D(\vec{k}) = \frac{2}{(2\pi)^2}$$	$$D(\vec{k}) = \frac{2}{(2\pi)^3}$$
$\,\,$Density of electron states $\,\,$in reciprocal space	$$D(|\vec{k}|) = \frac{2}{\pi}$$	$$D(|\vec{k}|) = \frac{k}{\pi}\quad\text{m}^{-1}$$	$$D(|\vec{k}|) = \frac{k^2}{\pi^2}\quad\text{m}^{-2}$$
$\,\,$Dispersion relation	$$E = \frac{\hbar^2 k^2}{2m}\quad\text{J}$$	$$E = \frac{\hbar^2 k^2}{2m}\quad\text{J}$$	$$E = \frac{\hbar^2 k^2}{2m}\quad\text{J}$$
$\,\,$Fermi wave vector  implicitly defined by$$\quad n = \int\limits_{-\infty}^{k_F}D(|\vec{k}|) dk$$	$$k_F = \frac{\pi n}{2}\quad\text{m}^{-1}$$	$$k_F = \sqrt{2\pi n}\quad\text{m}^{-1}$$	$$k_F = (3\pi^2 n)^{1/3}\quad\text{m}^{-1}$$
$\,\,$Fermi energy EF$$E_F = \frac{\hbar^2 k_F^2}{2m}$$	$$E_F = \frac{\hbar^2}{2m}\frac{\pi^2 n^2}{4}\quad\text{J}$$	$$E_F = \frac{\hbar^2}{2m}2\pi n\quad\text{J}$$	$$E_F = \frac{\hbar^2}{2m}(3\pi^2 n)^{2/3}\quad\text{J}$$
$\,\,$Density of states $D(E)$	$$D(E) = \frac{n}{2\sqrt{E_F}}\frac{1}{\sqrt{E}}\quad\text{J}^{-1}\text{ m}^{-1}$$	$$D(E) = \frac{n}{E_F}\quad\text{J}^{-1}\text{ m}^{-2}$$	$$D(E) = \frac{3}{2}\frac{n\sqrt{E}}{E_F^{3/2}}\quad\text{J}^{-1}\text{ m}^{-3}$$
$$D(E_F)$$	$$D(E_F) = \frac{n}{2E_F}\quad\text{J}^{-1}\text{ m}^{-1}$$	$$D(E_F) = \frac{n}{E_F}\quad\text{ J}^{-1}\text{ m}^{-2}$$	$$D(E_F) = \frac{3}{2}\frac{n}{E_F}\quad\text{ J}^{-1}\text{ m}^{-3}$$
$$D'(E_F)=\frac{dD(E)}{dE}\Big{|}_{E_F}$$	$$D'(E_F) =-\frac{1}{4}\frac{n}{E_F^2}\quad\text{J}^{-2}\text{ m}^{-1}$$	$$D'(E_F) =0\quad\text{J}^{-2}\text{ m}^{-2}$$	$$D'(E_F) =\frac{3}{4}\frac{n}{E_F^2}\quad\text{J}^{-2}\text{ m}^{-3}$$
$\,\,$Chemical potential μ  implicitly defined by$$\quad n = \int\limits_{-\infty}^{\infty}\frac{D(E)}{\exp{\left(\frac{E-\mu}{k_BT}\right)}+1}dE\quad$$	$$\mu(T) \approx E_F+\frac{\pi^2}{12}\frac{(k_BT)^2}{E_F}\quad\text{ J}$$	$$\mu(T) = k_BT\ln\left(\exp\left(\frac{E_F}{k_BT}\right) -1\right)\quad\text{ J}$$	$$\mu(T) \approx E_F - \frac{\pi^2}{12}\frac{(k_BT)^2}{E_F}\quad\text{ J}$$
$\,\,$Internal energy$$u =\int\limits_{-\infty}^{\infty}ED(E)f(E)dE$$	$$u \approx \frac{1}{3}nE_F +\frac{\pi^2}{12}\frac{n}{E_F}(k_BT)^2\quad\text{J m}^{-1}$$	$$u \approx \frac{1}{2}nE_F +\frac{\pi^2}{6}\frac{n}{E_F}(k_BT)^2\quad\text{J m}^{-2}$$	$$\quad u \approx \frac{3}{5}nE_F + \frac{\pi^2}{4}\frac{n}{E_F}(k_BT)^2\quad\text{J m}^{-3}\quad$$
$\,\,$Specific heat$$c_v = \left(\frac{\partial u}{\partial T}\right)_V $$	$$c_v \approx \frac{\pi^2}{6}\frac{n}{E_F}k_B^2T\quad\text{J K}^{-1}\text{ m}^{-1}$$	$$c_v \approx \frac{\pi^2}{3}\frac{n}{E_F}k_B^2T\quad\text{J K}^{-1}\text{ m}^{-2}$$	$$c_v \approx \frac{\pi^2}{2}\frac{n}{E_F}k_B^2T\quad\text{J K}^{-1}\text{ m}^{-3}$$
$\,\,$Entropy$$s =\int\frac{c_v}{T}dT$$	$$s \approx \frac{\pi^2}{6}\frac{n}{E_F}k_B^2T\quad\text{J K}^{-1}\text{ m}^{-1}$$	$$s \approx \frac{\pi^2}{3}\frac{n}{E_F}k_B^2T\quad\text{J K}^{-1}\text{ m}^{-2}$$	$$s \approx \frac{\pi^2}{2}\frac{n}{E_F}k_B^2T\quad\text{ J K}^{-1}\text{ m}^{-3}$$
$\,\,$Helmholtz free energy$$f = u -Ts$$	$$f \approx \frac{3}{5}nE_F - \frac{\pi^2}{12}\frac{n}{E_F}(k_BT)^2\quad\text{J m}^{-1}$$	$$f \approx \frac{3}{5}nE_F - \frac{\pi^2}{6}\frac{n}{E_F}(k_BT)^2\quad\text{J m}^{-2}$$	$$f \approx \frac{3}{5}nE_F - \frac{\pi^2}{4}\frac{n}{E_F}(k_BT)^2\quad\text{J m}^{-3}$$
$\,\,$Pressure$$P = -\left(\frac{\partial F}{\partial V}\right)_{N,T}$$	$$P \approx \frac{2}{3}nE_F + \frac{\pi^2}{6}\frac{n}{E_F}(k_BT)^2\quad\text{ N }$$	$$P \approx \frac{1}{2}nE_F + \frac{\pi^2}{6}\frac{n}{E_F}(k_BT)^2\quad\text{ N m}^{-1}$$	$$P \approx \frac{2}{5}nE_F + \frac{\pi^2}{6}\frac{n}{E_F}(k_BT)^2\quad\text{ N m}^{-2}$$
$\,\,$Bulk modulus$$B = -V\left(\frac{\partial P}{\partial V}\right)_{N,T}$$	$$B \approx 2nE_F - \frac{\pi^2}{6}\frac{n}{E_F}(k_BT)^2\quad\text{ N}$$	$$B \approx nE_F \quad\text{ N m}^{-1}$$	$$B \approx \frac{2}{3}nE_F + \frac{\pi^2}{18}\frac{n}{E_F}(k_BT)^2\quad\text{ N m}^{-2}$$
$\,\,$Enthalpy$$h = u + P$$	$$h \approx nE_F + \frac{\pi^2}{4}\frac{n}{E_F}(k_BT)^2 \quad\text{ J m}^{-1}$$	$$h \approx nE_F + \frac{\pi^2}{3}\frac{n}{E_F}(k_BT)^2 \quad\text{ J m}^{-2}$$	$$h \approx nE_F + \frac{5\pi^2}{12}\frac{n}{E_F}(k_BT)^2 \quad\text{ J m}^{-3}$$
$\,\,$Gibbs energy$$g = h-Ts$$	$$g \approx nE_F - \frac{\pi^2}{12}\frac{n}{E_F}(k_BT)^2\quad\text{ J m}^{-1}$$	$$g \approx nE_F \quad\text{ J m}^{-2}$$	$$g \approx nE_F - \frac{\pi^2}{12}\frac{n}{E_F}(k_BT)^2\quad\text{ J m}^{-3}$$

Below is an app that evaluates mathematical expressions that are typed into the textbox with the blue border. The formulas for the free-electron model have been entered.

---

  

<p>The mathematical functions that can be used are list below. Multiplication must be specified with a '*' symbol, 3*cos(x) not 3cos(x). Powers are specified with the 'pow' function: x³ is pow(x,3) or x**3 not x^3.</p> <table align='center'><tr><td><ul><li>abs(<i>x</i>) - absolute value</li><li>acos(<i>x</i>) - inverse cosine</li><li>acosh(<i>x</i>) - inverse hyperbolic cosine</li> <li>asin(<i>x</i>) - inverse sine</li><li>asinh(<i>x</i>) - inverse hyperbolic sine</li><li>atan(<i>x</i>) - inverse tangent</li><li>atanh(<i>x</i>) - inverse hyperbolic tangent</li><li>cos(<i>x</i>) - cosine</li><li>cosh(<i>x</i>) - hyperbolic cosine</li> <li>exp(<i>x</i>) - <i>e^t</i></li></ul></td> <td><ul><li>log(<i>x</i>) - natural logarithm</li><li>pi = 3.141592653589793</li><li>pow(<i>t,y</i>) - compute <i>t^y</i></li><li>round(<i>x</i>) - round to the nearest integer</li> <li>sin(<i>x</i>) - sine</li><li>sinh(<i>x</i>) - hyperbolic sine</li><li>sqrt(<i>x</i>) - square root</li><li>tan(<i>x</i>) - tangent</li><li>tanh(<i>x</i>) - hyperbolic tangent</li><li>H(<i>x</i>) - Heaviside function = 0 for <i>x</i> < 0, 1 for <i>x</i> > 0</li></ul></td></tr></table> <br /> <a href="https://lampz.tugraz.at/~hadley/physikm/script/intro/vectorscript.en.php?print">Vector functions</a><br /> <table align='center'><tr><td><ul><li>cross(A,B) - cross product</li><li>dot(A,B) - scalar product of two vectors </li><li>dot(a,B) - scalar a times vector B</li><li>length(A) - length of a vector</li></ul></td><td valign="top"><ul><li>unit(A) - unit vector in direction of A</li><li>vadd(A,B) - vector addition A + B</li><li>vsub(A,B) - vector subtraction A -B</li></ul></td></tr></table> <p>● print(x) - Text or variables can be printed to the Script Output box with the print(x) function. Text should be put in quotes, print("Print some text"). The value of variable x can be printed with print(x) and text can be mixed with variables by concatenating them together with a + sign, print("The length is "+x+" m."). </p>

me = 9.1093837015E-31; // [kg] electron mass kB = 1.380649E-23; // [J/K] Boltzmann's constant hbar = 1.054571817E-34; // [J s] reduced Planck's constant e = 1.602176634E-19; // [C] elementary charge pi = 3.141592653589793; n = 1E28; // electron density in [1/m^3] m = 1*me; // effective mass [kg] me = electron mass T = 300; // temperature [K] kF = pow(3*pi*pi*n,1/3); // Fermi wavenumber [1/m] EF = hbar*hbar*kF*kF/(2*m); //Fermi energy [J] EF_eV = EF/e; // Fermi energy [eV] mu = EF - (pi*pi/(6*EF))*pow(kB*T,2); //chemical potential [J] mu_eV = mu/e; //chemical potential [eV] u = 3*n*EF/5 + (pi*pi/4)*(n/EF)*pow(kB*T,2); // internal energy density J/m^3 cv = (pi*pi/2)*(n/EF)*(kB*kB*T); // specific heat J/(K m^3) s = (pi*pi/2)*(n/EF)*(kB*kB*T); // entropy density J/(K m^3) f = 3*n*EF/5 - (pi*pi/4)*(n/EF)*pow(kB*T,2); // Helmholtz free energy density J/m^3 P = 2*n*EF/5 + (pi*pi/6)*(n/EF)*pow(kB*T,2); // Pressure J/m^3 = N/m^2 B = 2*n*EF/3 + (pi*pi/18)*(n/EF)*pow(kB*T,2); // Bulk modulus J/m^3 = N/m^2 h = n*EF + (5*pi*pi/12)*(n/EF)*pow(kB*T,2); // enthalpy density J/m^3 g = n*EF - (pi*pi/12)*(n/EF)*pow(kB*T,2); // Gibbs energy density J/m^3
Script Output

print("3D free electron model"); me = 9.1093837015E-31; // [kg] electron mass kB = 1.380649E-23; // [J/K] Boltzmann's constant hbar = 1.054571817E-34; // [J s] reduced Planck's constant e = 1.602176634E-19; // [C] elementary charge pi = 3.141592653589793; n = 1E28; // electron density in [1/m^3] m = 1*me; // effective mass [kg] me = electron mass T = 300; // temperature [K] kF = pow(3*pi*pi*n,1/3); // Fermi wavenumber [1/m] EF = hbar*hbar*kF*kF/(2*m); //Fermi energy [J] EF_eV = EF/e; // Fermi energy [eV] mu = EF - (pi*pi/(6*EF))*pow(kB*T,2); //chemical potential [J] mu_eV = mu/e; //chemical potential [eV] u = 3*n*EF/5 + (pi*pi/4)*(n/EF)*pow(kB*T,2); // internal energy density J/m^3 cv = (pi*pi/2)*(n/EF)*(kB*kB*T); // specific heat J/(K m^3) s = (pi*pi/2)*(n/EF)*(kB*kB*T); // entropy density J/(K m^3) f = 3*n*EF/5 - (pi*pi/4)*(n/EF)*pow(kB*T,2); // Helmholtz free energy density J/m^3 P = 2*n*EF/5 + (pi*pi/6)*(n/EF)*pow(kB*T,2); // Pressure J/m^3 = N/m^2 B = 2*n*EF/3 + (pi*pi/18)*(n/EF)*pow(kB*T,2); // Bulk modulus J/m^3 = N/m^2 h = n*EF + (5*pi*pi/12)*(n/EF)*pow(kB*T,2); // enthalpy density J/m^3 g = n*EF - (pi*pi/12)*(n/EF)*pow(kB*T,2); // Gibbs energy density J/m^3 print("This code estimates some properties of lithium using the free electron model"); print("define some constants"); me = 9.1093837015E-31; // [kg] electron mass kB = 1.380649E-23; // [J/K] Boltzmann's constant hbar = 1.054571817E-34; // [J s] reduced Planck's constant e = 1.602176634E-19; // [C] elementary charge pi = 3.141592653589793; u = 1.6605402E-27; // [kg] atomic mass constant print("calculate the electron density"); m_Li = 6.941*u; // mass of a lithium atom rho = 0.534; // g/cm^3 density of lithium n_Li = rho*1000/m_Li; //atomic density in atoms/m^3 n = n_Li; // [1/m^3] electron density. Lithium has a valence of 1 print("Hall measurements support this value of the electron density see Phys. Rev. vol. 185, p. 876 (1969)"); m = 2.18*me; // effective mass [kg] from specific heat measurements (Kittel) T = 300; // temperature [K] kF = pow(3*pi*pi*n,1/3); // Fermi wavenumber [1/m] EF = hbar*hbar*kF*kF/(2*m); //Fermi energy [J] EF_eV = EF/e; // Fermi energy [eV] mu = EF - (pi*pi/(6*EF))*pow(kB*T,2); //chemical potential [J] mu_eV = mu/e; //chemical potential [eV] u0 = 3*n*EF/5; // internal energy density @ T=0 J/m^3 u = 3*n*EF/5 + (pi*pi/4)*(n/EF)*pow(kB*T,2); // internal energy density J/m^3 u_T = u - u0; // internal energy density J/m^3 print("u_T = u - u0 @ T= 300 K = 4177893 [J/m^3] from http://lampz.tugraz.at/~hadley/ss1/materials/thermo/dos2u.html"); cv = (pi*pi/2)*(n/EF)*(kB*kB*T); // specific heat J/(K m^3) print("cv = 27862 [J/(K m^3)] from http://lampz.tugraz.at/~hadley/ss1/materials/thermo/dos2cv.html"); s = (pi*pi/2)*(n/EF)*(kB*kB*T); // entropy density J/(K m^3) f = 3*n*EF/5 - (pi*pi/4)*(n/EF)*pow(kB*T,2); // Helmholtz free energy density J/m^3 P = 2*n*EF/5 + (pi*pi/6)*(n/EF)*pow(kB*T,2); // Pressure J/m^3 = N/m^2 B = 2*n*EF/3 + (pi*pi/18)*(n/EF)*pow(kB*T,2); // Bulk modulus J/m^3 = N/m^2 h = n*EF + (5*pi*pi/12)*(n/EF)*pow(kB*T,2); // enthalpy density J/m^3 g = n*EF - (pi*pi/12)*(n/EF)*pow(kB*T,2); // Gibbs energy density J/m^3 print("2D free electron model"); me = 9.1093837015E-31; // [kg] electron mass kB = 1.380649E-23; // [J/K] Boltzmann's constant hbar = 1.054571817E-34; // [J s] reduced Planck's constant e = 1.602176634E-19; // [C] elementary charge pi = 3.141592653589793; n = 1E19; // electron density in [1/m^2] m = 1*me; // effective mass [kg] me = electron mass T = 300; // temperature [K] kF = sqrt(2*pi*n); // Fermi wavenumber [1/m] EF = hbar*hbar*kF*kF/(2*m); //Fermi energy [J] EF_eV = EF/e; // Fermi energy [eV] mu = kB*T*log(exp(EF/(kB*T))-1); //chemical potential [J] mu_eV = mu/e; //chemical potential [eV] u = n*EF/2 + (pi*pi/6)*(n/EF)*pow(kB*T,2); // internal energy density J/m^2 cv = (pi*pi/3)*(n/EF)*(kB*kB*T); // specific heat J/(K m^2) s = (pi*pi/3)*(n/EF)*(kB*kB*T); // entropy density J/(K m^2) f = 3*n*EF/5 - (pi*pi/6)*(n/EF)*pow(kB*T,2); // Helmholtz free energy density J/m^2 P = n*EF/2 + (pi*pi/6)*(n/EF)*pow(kB*T,2); // Pressure J/m^2 = N/m B = n*EF; // Bulk modulus J/m^2 = N/m h = n*EF + (pi*pi/3)*(n/EF)*pow(kB*T,2); // enthalpy density J/m^2 g = n*EF; // Gibbs energy density J/m^2 print("1D free electron model"); me = 9.1093837015E-31; // [kg] electron mass kB = 1.380649E-23; // [J/K] Boltzmann's constant hbar = 1.054571817E-34; // [J s] reduced Planck's constant e = 1.602176634E-19; // [C] elementary charge pi = 3.141592653589793; n = 1E9; // electron density in [1/m] m = 1*me; // effective mass [kg] me = electron mass T = 300; // temperature [K] kF = pi*n/2; // Fermi wavenumber [1/m] EF = hbar*hbar*kF*kF/(2*m); //Fermi energy [J] EF_eV = EF/e; // Fermi energy [eV] mu = EF + (pi*pi/(12*EF))*pow(kB*T,2); //chemical potential [J] mu_eV = mu/e; //chemical potential [eV] u = n*EF/3 + (pi*pi/12)*(n/EF)*pow(kB*T,2); // internal energy density J/m cv = (pi*pi/6)*(n/EF)*(kB*kB*T); // specific heat J/(K m) s = (pi*pi/6)*(n/EF)*(kB*kB*T); // entropy density J/(K m) f = 3*n*EF/5 - (pi*pi/12)*(n/EF)*pow(kB*T,2); // Helmholtz free energy density J/m P = 2*n*EF/3 + (pi*pi/6)*(n/EF)*pow(kB*T,2); // Pressure J/m^2 = N B = 2*n*EF - (pi*pi/6)*(n/EF)*pow(kB*T,2); // Bulk modulus J/m = N h = n*EF + (pi*pi/4)*(n/EF)*pow(kB*T,2); // enthalpy density J/m g = n*EF - (pi*pi/12)*(n/EF)*pow(kB*T,2); // Gibbs energy density J/m print("This code estimates some properties of copper using the free electron model"); print("define some constants"); me = 9.1093837015E-31; // [kg] electron mass kB = 1.380649E-23; // [J/K] Boltzmann's constant hbar = 1.054571817E-34; // [J s] reduced Planck's constant e = 1.602176634E-19; // [C] elementary charge pi = 3.141592653589793; u = 1.6605402E-27; // [kg] atomic mass constant print("calculate the electron density"); m_Cu = 63.55*u; // mass of a copper atom rho = 8.933; // g/cm^3 density of copper n_Cu = rho*1000/m_Cu; //atomic density in atoms/m^3 n = n_Cu*1; // [1/m^3] electron density. Copper has 1 valence electron m = 1.38*me; // effective mass [kg] from specific heat measurements (Kittel) T = 300; // temperature [K] kF = pow(3*pi*pi*n,1/3); // Fermi wavenumber [1/m] EF = hbar*hbar*kF*kF/(2*m); //Fermi energy [J] EF_eV = EF/e; // Fermi energy [eV] mu = EF - (pi*pi/(6*EF))*pow(kB*T,2); //chemical potential [J] mu_eV = mu/e; //chemical potential [eV] u0 = 3*n*EF/5; // internal energy density @ T=0 J/m^3 u = 3*n*EF/5 + (pi*pi/4)*(n/EF)*pow(kB*T,2); // internal energy density J/m^3 u_T = u - u0; // internal energy density J/m^3 print("u_T = u - u0 @ T= 300 K = 4425074 [J/m^3] from http://lampz.tugraz.at/~hadley/ss1/materials/thermo/dos2u.html"); cv = (pi*pi/2)*(n/EF)*(kB*kB*T); // specific heat J/(K m^3) print("cv = 29048 [J/(K m^3)] from http://lampz.tugraz.at/~hadley/ss1/materials/thermo/dos2cv.html"); s = (pi*pi/2)*(n/EF)*(kB*kB*T); // entropy density J/(K m^3) f = 3*n*EF/5 - (pi*pi/4)*(n/EF)*pow(kB*T,2); // Helmholtz free energy density J/m^3 P = 2*n*EF/5 + (pi*pi/6)*(n/EF)*pow(kB*T,2); // Pressure J/m^3 = N/m^2 B = 2*n*EF/3 + (pi*pi/18)*(n/EF)*pow(kB*T,2); // Bulk modulus J/m^3 = N/m^2 h = n*EF + (5*pi*pi/12)*(n/EF)*pow(kB*T,2); // enthalpy density J/m^3 g = n*EF - (pi*pi/12)*(n/EF)*pow(kB*T,2); // Gibbs energy density J/m^3 print("This code estimates some properties of aluminium using the free electron model"); print("define some constants"); me = 9.1093837015E-31; // [kg] electron mass kB = 1.380649E-23; // [J/K] Boltzmann's constant hbar = 1.054571817E-34; // [J s] reduced Planck's constant e = 1.602176634E-19; // [C] elementary charge pi = 3.141592653589793; u = 1.6605402E-27; // [kg] atomic mass constant print("calculate the electron density"); m_Al = 26.982*u; // mass of an aluminium atom rho = 2.70; // g/cm^3 density of aluminium n_Al = rho*1000/m_Al; //atomic density in atoms/m^3 n = n_Al*3; // [1/m^3] electron density. Aluminium has 3 valence electrons m = 1.48*me; // effective mass [kg] from specific heat measurements (Kittel) T = 300; // temperature [K] kF = pow(3*pi*pi*n,1/3); // Fermi wavenumber [1/m] EF = hbar*hbar*kF*kF/(2*m); //Fermi energy [J] EF_eV = EF/e; // Fermi energy [eV] mu = EF - (pi*pi/(6*EF))*pow(kB*T,2); //chemical potential [J] mu_eV = mu/e; //chemical potential [eV] u0 = 3*n*EF/5; // internal energy density @ T=0 J/m^3 u = 3*n*EF/5 + (pi*pi/4)*(n/EF)*pow(kB*T,2); // internal energy density J/m^3 u_T = u - u0; // internal energy density J/m^3 print("u_T = u - u0 @ T= 300 K = 4308538 [J/m^3] from http://lampz.tugraz.at/~hadley/ss1/materials/thermo/dos2u.html"); cv = (pi*pi/2)*(n/EF)*(kB*kB*T); // specific heat J/(K m^3) print("cv = 28467 [J/(K m^3)] from http://lampz.tugraz.at/~hadley/ss1/materials/thermo/dos2cv.html"); s = (pi*pi/2)*(n/EF)*(kB*kB*T); // entropy density J/(K m^3) f = 3*n*EF/5 - (pi*pi/4)*(n/EF)*pow(kB*T,2); // Helmholtz free energy density J/m^3 P = 2*n*EF/5 + (pi*pi/6)*(n/EF)*pow(kB*T,2); // Pressure J/m^3 = N/m^2 B = 2*n*EF/3 + (pi*pi/18)*(n/EF)*pow(kB*T,2); // Bulk modulus J/m^3 = N/m^2 h = n*EF + (5*pi*pi/12)*(n/EF)*pow(kB*T,2); // enthalpy density J/m^3 g = n*EF - (pi*pi/12)*(n/EF)*pow(kB*T,2); // Gibbs energy density J/m^3 print("This code estimates some properties of gold using the free electron model"); print("define some constants"); me = 9.1093837015E-31; // [kg] electron mass kB = 1.380649E-23; // [J/K] Boltzmann's constant hbar = 1.054571817E-34; // [J s] reduced Planck's constant e = 1.602176634E-19; // [C] elementary charge pi = 3.141592653589793; u = 1.6605402E-27; // [kg] atomic mass constant print("calculate the electron density"); m_Au = 196.967*u; // mass of a gold atom rho = 19.282; // g/cm^3 density of gold n_Au = rho*1000/m_Au; //atomic density in atoms/m^3 n = n_Au*1; // [1/m^3] electron density. Gold has 1 valence electron m = 1.26*me; // effective mass [kg] from specific heat measurements (Kittel) T = 300; // temperature [K] kF = pow(3*pi*pi*n,1/3); // Fermi wavenumber [1/m] EF = hbar*hbar*kF*kF/(2*m); //Fermi energy [J] EF_eV = EF/e; // Fermi energy [eV] mu = EF - (pi*pi/(6*EF))*pow(kB*T,2); //chemical potential [J] mu_eV = mu/e; //chemical potential [eV] u0 = 3*n*EF/5; // internal energy density @ T=0 J/m^3 u = 3*n*EF/5 + (pi*pi/4)*(n/EF)*pow(kB*T,2); // internal energy density J/m^3 u_T = u - u0; // internal energy density J/m^3 print("u_T = u - u0 @ T= 300 K = 2951849 [J/m^3] from http://lampz.tugraz.at/~hadley/ss1/materials/thermo/dos2u.html"); cv = (pi*pi/2)*(n/EF)*(kB*kB*T); // specific heat J/(K m^3) print("cv = 19354 [J/(K m^3)] from http://lampz.tugraz.at/~hadley/ss1/materials/thermo/dos2cv.html"); s = (pi*pi/2)*(n/EF)*(kB*kB*T); // entropy density J/(K m^3) f = 3*n*EF/5 - (pi*pi/4)*(n/EF)*pow(kB*T,2); // Helmholtz free energy density J/m^3 P = 2*n*EF/5 + (pi*pi/6)*(n/EF)*pow(kB*T,2); // Pressure J/m^3 = N/m^2 B = 2*n*EF/3 + (pi*pi/18)*(n/EF)*pow(kB*T,2); // Bulk modulus J/m^3 = N/m^2 h = n*EF + (5*pi*pi/12)*(n/EF)*pow(kB*T,2); // enthalpy density J/m^3 g = n*EF - (pi*pi/12)*(n/EF)*pow(kB*T,2); // Gibbs energy density J/m^3