The free electron model is a simple model for the electronic properties of metals. It is based on the particle-in-a-box problem familiar from an introductory quantum-mechanics course except we add many particles to the box. The electrons are 'free' because we neglect the Coulomb interactions between them. Below we will calculate thermodynamic properties of this system on noninteracting fermions such as the internal energy and the specific heat. When this model is applied to metals, the fermion density and the mass of the noninteracting fermions are treated as fit parameters. You might think that the electron density is easily calculated because we know the atomic density of metals like copper or aluminum and we know how many electrons each atom has. However, the core electrons typically do not contribute to the thermodynamic properties. The reason for this will become clear when we examine the electron band model for solids.

Consider $N$ noninteracting fermions confined to a box with dimensions $L\times L\times L$. The potential energy in the box is zero so the Schrödinger equation for this case just includes the kinetic energy terms,

$$\left(- \sum \limits_i^N \frac{\hbar^2}{2m}\nabla_i^2\right) \Psi(\vec{r}_1,\vec{r}_2,\cdots,\vec{r}_N) = E \Psi(\vec{r}_1,\vec{r}_2,\cdots,\vec{r}_N).$$

Here $i$ sums over the particles. This many-particle equation separates into $N$ single-particle equations and each of the single-particle equations is the same, $$- \frac{\hbar^2}{2m}\nabla_i^2 \psi_i(\vec{r}_i) = E_i \psi_i(\vec{r}_i).$$

The many-particle wavefunction is a product of the single-particle spin orbitals and the total energy is the sum of the single-particle energies,

$$\Psi(\vec{r}_1,\vec{r}_2,\cdots,\vec{r}_N)= \mathcal{A}\left(\psi_1(\vec{r}_1)\uparrow\psi_1(\vec{r}_1)\downarrow\psi_2(\vec{r}_2)\uparrow\psi_2(\vec{r}_2)\downarrow\cdots \right),$$ $$E= E_1\uparrow + E_1\downarrow + E_2\uparrow + E_2\downarrow + \cdots.$$

Here $\mathcal{A}$ is the antisymmetrizing operator that ensures that if particles $j$ and $k$ are exchanged, the total wavefunction changes sign.

Next, we need to solve the single particle Schrödinger equation. The particle-in-a-box problem is usually solved with fixed boundary conditions of $\psi =0$ at the boundaries. Here we will use periodic boundary conditions in each coordinate, 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 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 and it 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(k)dk.$$

Here we have defined

$$D(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 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(k)dk\qquad\rightarrow\qquad D(E) = D(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. 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.

$$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$. 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 and performing 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 probablility 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 $$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}^{-3}$$ $$h \approx nE_F + \frac{\pi^2}{3}\frac{n}{E_F}(k_BT)^2 \quad\text{ J m}^{-3}$$ $$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.

 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