The free electron gas is based on the following assumptions: no electron-electron interactions, the electrons experience no potential due to the crystal, they are confined to a box. In addition, we apply a constant homogeneous external magnetic field.
First, we restrict the discussion to the coupling of the electronic spin to the magnetic field, i.e. we ignore the angular moment. The energy one-particle eigenvalues are \[\begin{align*} \varepsilon_{\sigma}(\vv k) &= \frac{\hbar^{2} \vv k^{2}}{2m} +\sigma b\;,\\ b &= \frac{\mu_{B} g_{e} }{2}B\;.\end{align*}\] with the quantized wave vectors \(\vv k\). The mean occupation of the one-particle orbitals is given by the Fermi-Dirac distribution \[\begin{align*} n_{F}(\varepsilon_{\sigma}(\vv k)|T,\mu)\end{align*}\] The mean total number of electrons in the Zeeman level represented by \(\sigma\) is \[\begin{align*} N_{\sigma} &=\sum_{\vv k} n_{F}(\varepsilon_{\sigma}(\vv k)|T,\mu) = \int d\varepsilon \rho(\varepsilon) n_{F}(\varepsilon +\sigma b|T,\mu)\;.\end{align*}\] As derived in the 3D dos is
\[\begin{align} \rho(\varepsilon)&= D \sqrt{\varepsilon}\;,\\ \text{with}\quad D&=\frac{V m^{3/2}}{\hbar^{3} \pi^{2}\sqrt{2}}\;.\end{align}\]
Then \[\begin{align*} N_{\sigma} &=D\; \int_{0}^{\infty} d\varepsilon \sqrt{\varepsilon} \; n_{F}(\varepsilon +\sigma b|T,\mu)\;.\end{align*}\] For small magnetic field we can use a Taylor expansion in \(b\) abound \(b=0\) \[\begin{align*} N_{\sigma} &=D\; \int_{0}^{\infty} d\varepsilon \sqrt{\varepsilon} \; \bigg(n_{F}(\varepsilon|T,\mu) +\sigma b \;\frac{\partial }{\partial \varepsilon}n_{F}(\varepsilon|T,\mu) +{\cal O}(b^{2})\bigg)\\ &=D\; \int_{0}^{\infty} d\varepsilon \sqrt{\varepsilon} \; n_{F}(\varepsilon|T,\mu) +\sigma b D \int_{0}^{\infty} d\varepsilon \sqrt{\varepsilon} \;n'_{F}(\varepsilon|T,\mu) +{\cal O}(b^{2})\end{align*}\] The total particle number follows as \[\begin{align} N &= N_{+}+N_{-} = 2 D \int_{0}^{\infty} d\varepsilon \sqrt{\varepsilon} \;n_{F} (\varepsilon|T,\mu) + {\cal O}(b^{2})\notag\\ &= 2 D \int_{0}^{\infty} d\varepsilon \sqrt{\varepsilon} \;n_{F} (\varepsilon|T,\mu) + {\cal O}(b^{2})\;.\label{eq:N:total}\end{align}\] For small \(b\) and low temperature we can use the Sommerfeld expansion, outlined in ,
\[\begin{align} I&=\int_{0}^{\infty} f(\varepsilon) n_{F}(\varepsilon|\mu,T) d\varepsilon\nonumber\\ &=\int_{0}^{\mu} f(\varepsilon)d\varepsilon + 2 \sum_{n=1}^{\text{odd}}\big(1-\frac{1}{2^{n}}\big)\;\zeta(n+1) \big(k_{B}T\big)^{n+1}\;f^{(n)}(\mu)\;,\end{align}\]
which reads to leading order \[\begin{align*} I &= \int_{0}^{\mu} f(\varepsilon) d\varepsilon + \frac{\pi^{2}}{6}\;\big(k_{B}T\big)^{2} \;f'(\mu)+ {\cal O}(\big(\frac{k_{B}T}{\mu}\big)^{4})\;.\end{align*}\] For the total particle number we obtain \[\begin{align*} N&=2 D \bigg(\int_{0}^{\mu}\sqrt{\varepsilon} +\frac{\pi^{2}}{6}\frac{d}{d\mu}\sqrt{\mu}\bigg( k_{B}T \bigg)^{2}+\ldots\bigg)\\ &=2 D \bigg(\frac{2}{3}\mu^{3/2} +\frac{\pi^{2}}{12}\mu^{-1/2}\big( k_{B}T \big)^{2}+\ldots \bigg)\;.\end{align*}\] \[\begin{align} \label{eq:spin:para:aux2} \frac{3 N}{4 D}&= \mu^{3/2}\bigg(1 +\frac{\pi^{2}}{8}\bigg(\frac{ k_{B}T }{\mu}\bigg)^{2}+\ldots \bigg)\end{align}\]
For \(T=0\) the chemical potential is equivalent to the Fermi energy \(\varepsilon_{F}\) and we have \[\begin{align} \label{eq:spin:para:aux1} \frac{3 N}{4 D} &= \varepsilon_{F}^{3/2}\;,\end{align}\] or rather \[\begin{align*} \varepsilon_{F} &= \bigg(\frac{3 N \hbar^{3} (2\pi)^{2} }{4 V(2m)^{3/2}} \bigg)^{2/3} = \bigg(\frac{3 N \pi^{2} }{ V} \bigg)^{2/3}\;\frac{\hbar^{2}}{2 m}\;.\end{align*}\]
\[\begin{align} \label{eq:spin:para:EF} \varepsilon_{F}&= \bigg(3 \pi^{2} n\bigg)^{2/3} \frac{\hbar^{2}}{2m} \;.\end{align}\]
This defines the Fermi wave number \(k_{F}\), through \[\begin{align*} \varepsilon_{F} & = \bigg(3 \pi^{2} n\bigg)^{2/3}\frac{\hbar^{2}}{2m} = \frac{\hbar^{2} k_{F}^{2} }{2m}\\ k_{F}&=\bigg(3 \pi^{2} n\bigg)^{1/3} = \bigg(\frac{3 \pi^{2}}{v}\bigg)^{1/3 }\;,\end{align*}\] where \(v\) is the average volume per electron. Hence \[\begin{align*} k_{F} &\propto \frac{1}{r}\;,\end{align*}\] where \(r\) is the mean distance between the electrons. Inserting \(\eqref{eq:spin:para:aux1}\) in \(\eqref{eq:spin:para:aux2}\) yields apart from higher order terms \[\begin{align*} \varepsilon_{F} &= \mu\bigg(1 +\frac{\pi^{2}}{8}\bigg(\frac{ k_{B}T }{\mu}\bigg)^{2}\bigg)^{2/3}\\ \mu &= \varepsilon_{F}\bigg(1 +\frac{\pi^{2}}{8}\bigg(\frac{ k_{B}T }{\mu}\bigg)^{2}\bigg)^{-2/3}\;.\end{align*}\] We can solve this equation iteratively, starting with \(\mu=\varepsilon_{F}\). The first iteration yields \[\begin{align*} \mu &=\varepsilon_{F}\bigg(1 +\frac{\pi^{2}}{8}\bigg(\frac{ k_{B}T }{\varepsilon_{F}}\bigg)^{2}\bigg)^{-2/3}\;.\end{align*}\] For low temperatures, \(\frac{ k_{B}T }{\varepsilon_{F}}<1\) we can as well write \[\begin{align*} \mu &=\varepsilon_{F}\bigg(1 -\frac{\pi^{2}}{8}\frac{2}{3}\bigg(\frac{ k_{B}T }{\varepsilon_{F}}\bigg)^{2} +{\cal O}\bigg( \big( \frac{k_{B}T}{\varepsilon_{F}} \big)^{4} \bigg)\bigg)\;.\end{align*}\] Further iterations do not change the second order term and we generally have
\[\begin{align} \mu &=\varepsilon_{F}\bigg(1 -\frac{\pi^{2}}{12}\bigg(\frac{ k_{B}T }{\varepsilon_{F}}\bigg)^{2} \bigg)\;+{\cal O}\bigg( \big( \frac{k_{B}T}{\varepsilon_{F}} \big)^{4} \bigg)\end{align}\]
The partition function for non-interacting particles has already been derived previously. For one-particle energies \(\varepsilon_{\nu}\) it was \[\begin{align*} \ln(Z) &= \sum_{\nu} \ln\bigg( 1+ e^{-\beta(\varepsilon_{\nu}-\mu)}\bigg)\;.\end{align*}\] In the present case the index \(\nu\) stands for the wave vector \(\vv k\) and the spin direction. Hence \[\begin{align} \label{eq:para:pauli:ln:Z} \ln(Z) %&= \sum_{\sigma} \sum_{\vv k} \ln\bigg( 1+ e^{-\beta(\varepsilon(\vv k)+\sigma b-\mu)}\bigg)\\ &= \sum_{\sigma} \int_{0}^{\infty}\;\rho(\varepsilon)\; \ln\bigg( 1+ e^{-\beta(\varepsilon+\sigma b-\mu)}\bigg)\;.\end{align}\]
The corresponding grand potential reads \[\begin{align*} \Omega(T,\vv B) &= -k_{B}T \ln(Z) =-k_{B}T \sum_{\sigma} \int_{0}^{\infty}\;\rho(\varepsilon)\; \ln\bigg( 1+ e^{-\beta(\varepsilon+\sigma b-\mu)}\bigg)\;.\end{align*}\] The magnetization in z-direction is obtain via \[\begin{align*} M &= - \pder{\Omega}{B}{T} \\ &=k_{B}T\;\sum_{\sigma}\big( -\frac{\beta\sigma \mu_{B}g_{e}}{2} \big) \int_{0}^{\infty} d\varepsilon\;\rho(\varepsilon)\;\frac{e^{-\beta(\varepsilon+\sigma b-\mu)}}{1+e^{-\beta(\varepsilon+\sigma b-\mu)}}\\ &=-\frac{\mu_{B}g_{e}}{2} \;\sum_{\sigma}\sigma \;\; \int_{0}^{\infty} d\varepsilon\;\rho(\varepsilon)\;n_{F}(\varepsilon+\sigma b|T,\mu)\\ &= -\frac{\mu_{B}g_{e}}{\hbar} \;\underbrace{ \frac{\hbar(N_{+}-N_{-})}{2} }_{\color{blue} = \langle S_\text{total}^{z} \rangle}\end{align*}\] In agreement with \(\eqref{eq:magmo:spin}\). Again we assume that \(b\) is small and employ a Taylor expansion. \[\begin{align*} M &= -\frac{\mu_{B}g_{e} }{2} \sum_{\sigma} \sigma \int_{0}^{\infty}d\varepsilon \bigg(\rho(\varepsilon) n_{F}(\varepsilon|T,\mu) + \sigma b n'_{F}(\varepsilon|T,\mu)+{O}(b^{2})\bigg)\\ &= -\frac{\mu_{B}g_{e} }{2} 2b \int_{0}^{\infty}d\varepsilon \rho(\varepsilon) n'_{F}(\varepsilon|T,\mu)+{O}(b^{2})\\ &= -\frac{(\mu_{B}g_{e})^{2} }{2} B\int_{0}^{\infty}d\varepsilon \rho(\varepsilon) n'_{F}(\varepsilon|T,\mu)+{O}(b^{2})\;.\end{align*}\] Then the susceptibility reads \[\begin{align} \chi_{T} &= \mu_{0}\pder{M}{B}{T}\bigg|_{B=0}\notag\\ &=-\mu_{0}\frac{(\mu_{B}g_{e})^{2} }{2} \int_{0}^{\infty}d\varepsilon \rho(\varepsilon) n'_{F}(\varepsilon|T,\mu)\notag\;.\end{align}\] This can also be written as \[\begin{align*} \chi_{T} &= -\mu_{0} \mu_{B}^{2} \bigg(\frac{g_{e}}{2}\bigg)^{2} 2 \int d\varepsilon \rho(\varepsilon) n'_{F} (\varepsilon|T,\mu)\\ &= \mu_{0} \mu_{B}^{2} \bigg(\frac{g_{e}}{2}\bigg)^{2} 2 \frac{\partial }{\partial \mu} \int d\varepsilon \rho(\varepsilon) n_{F} (\varepsilon|T,\mu)\;.\end{align*}\] Comparison with \(\eqref{eq:N:total}\) yields \[\begin{align} \label{eq:Pauli:dN:dmu} \chi_{T} &= \mu_{0} \mu_{B}^{2} \bigg(\frac{g_{e}}{2}\bigg)^{2} \pder{N}{\mu}{T,B=0}\end{align}\]
We use the Sommerfeld expansion, derived in to expand the integral in powers of \(k_{B}T/\mu\). \[\begin{align*} \chi &=-\mu_{0}\frac{(\mu_{B}g_{e})^{2} }{2} \int_{0}^{\infty}d\varepsilon \rho(\varepsilon) n'_{F}(\varepsilon|T,\mu)\notag\\ &=\mu_{0}\frac{(\mu_{B}g_{e})^{2} }{2} \int_{0}^{\infty}d\varepsilon \rho'(\varepsilon) n_{F}(\varepsilon|T,\mu)\\ &=\mu_{0}\frac{(\mu_{B}g_{e})^{2} }{2} \bigg(\int_{0}^{\mu}d\varepsilon \rho'(\varepsilon) +\frac{\pi^{2}}{6}(k_{B}T)^{2}\rho''(\mu)+ {\cal O}(\big(\frac{k_{B}T}{\mu}\big)^{4}) \bigg)\\ &=\mu_{0}\frac{(\mu_{B}g_{e})^{2} }{2} \bigg( \rho(\mu) +\frac{\pi^{2}}{6}(k_{B}T)^{2}\rho''(\mu)+ {\cal O}(\big(\frac{k_{B}T}{\mu}\big)^{4}) \bigg)\;.\end{align*}\] The final result reads \[\begin{align} \label{eq:para:chi} \chi &=\mu_{0}\frac{(\mu_{B}g_{e})^{2} }{2} \rho(\mu)\bigg( 1 +\frac{\pi^{2}}{6}(k_{B}T)^{2}\frac{\rho''(\mu)}{\rho(\mu)}+\ldots \bigg)\end{align}\] Now we recall, that \(\rho(\varepsilon) = D \sqrt{\varepsilon}\), resulting in \(\rho''(\mu)/\rho(\mu)=-1/(4\mu^{2})\).
\[\begin{align} \chi_{P} &= \mu_{0}\frac{(\mu_{B}g_{e})^{2} }{2} \rho(\mu) \bigg( 1 - \frac{\pi^{2}}{24}\bigg(\frac{k_{B}T}{\mu}\bigg)^{2} +\ldots\bigg)\;.\end{align}\]
We see that indeed \(k_{B}T/\mu\) is the relevant small parameter. This describes the Pauli-spin-paramagnetism, which is almost temperature independent for low \(T\); in strong contrast to the Curie \(1/T\) behaviour. The reason for the discrepancy lies in the fermi-statistics. Only the spin of the electrons in the vicinity of \(\mu\) can contribute. The number of thermally excited electrons is \(k_{b} T \rho(\mu)\) which compensates the \(1/T\) behaviour.
So far we have only considered the spin degrees of freedom of the free electron gas. The orbital moments also contribute to the magnetization, which results in the
\[\begin{align} \chi_{L}&= -\frac{1}{3}\chi_{P}\end{align}\]
For the derivation see appendix [@app:free:electron:gas]