This plot shows band structure of graphite along some high symmetry lines using a DFT calcualtion.

In the following figure the resulting density of states (DOS) are plotted.

This plot shows the density of states of graphite as given in the DOSCAR file of the DFT calculation using VASP.

The Fermi energy is at $E_F = 5.04 \; eV$. Since the DOS are nearly symmetric around the Fermi energy, the chemical potential is almost constant with temperature and for further calculations the chemical potential is approximated by the Fermi energy. The calculated $E(\vec{k})$ dispersion relation from the EIGENVAL file is read in by the following Matlab function.

The resulting density of states are plotted in the following figure.

This plot shows the density of states of graphite derived from the tight binding model.

The resulting fit parameters are $p_0 = 3.27e+3 \; A$, $p_1 = 5.34e+5 \; AK$, $p_2 = 9.65e+6 \; AK^2$, $p_3 = 3.24e+8 \; AK^3$, $p_4 = -2.94e-3 \; 1/K$. In the following figure the tabulated values are plotted together with the assumed function.

This plot shows the tabulated values of the phonon contribution of the free path length and the fitted function.

As the assumed function describes the tabulated values quite well, this function is used interpolate between the tabulated values. Now the tensor elements can be calculated by using the chosen discretization of $b_1$ = (-20,...,20)/41, $b_2$ = (-20,...,20)/41, $b_3$ = (-5,...,6)/12 from the DFT calculation for the numerical integration. The tensor elements are then given by \begin{equation} \sigma_{22} = \frac{e^2}{\pi h^2} \frac{\sqrt{3} \pi}{c} \frac{1}{41 \cdot 41 \cdot 12} \sum_{bands,j} \sum_{b_1} \sum_{b_2} \sum_{b_3} \left(\frac{v(\vec{k})}{l_{def}} + \frac{v(\vec{k})}{l_{ph}(T)}\right)^{-1} \frac{\partial f_0}{\partial \mu} \left( \frac{\partial}{\partial b_2}E_j(\vec{b})\right)^2 \;, \end{equation} \begin{equation} \sigma_{33} = \frac{e^2}{\pi h^2} \frac{4 \pi c}{\sqrt{3} a^2} \frac{1}{41 \cdot 41 \cdot 12} \sum_{bands,j} \sum_{b_1} \sum_{b_2} \sum_{b_3} \left(\frac{v(\vec{k})}{l_{def}} + \frac{v(\vec{k})}{l_{ph}(T)}\right)^{-1} \frac{\partial f_0}{\partial \mu} \left(\frac{\partial}{\partial b_2}E_j(\vec{b})\right)^2 \;, \end{equation} with the velocity $v(\vec{k})$ approximated by \begin{equation} v(\vec{k}) = \frac{2 \pi}{h} \left|\vec{\nabla}_{\vec{k}} E(\vec{k}) \right| = \frac{2 \pi}{h} \sqrt{ \left(\frac{\partial E(\vec{k})}{\partial k_x} \right)^2 + \left(\frac{\partial E(\vec{k})}{\partial k_y} \right)^2 + \left(\frac{\partial E(\vec{k})}{\partial k_z} \right)^2} = \frac{1}{h} \sqrt{ a^2 \left(\frac{\partial E(\vec{b})}{\partial b_1} - \frac{1}{2} \frac{\partial E(\vec{b})}{\partial b_2}\right)^2 + \frac{3 a^2}{4} \left(\frac{\partial E(\vec{b})}{\partial b_2} \right)^2 + c^2 \left(\frac{\partial E(\vec{b})}{\partial b_3} \right)^2}\;, \end{equation} and the partial derivatives approximated by \begin{equation} \frac{\partial}{\partial b_1}E_j(b_1(i),b_2(j),b_3(m)) = \frac{E_j(b_1(i+1),b_2(j),b_3(m)) - E_j(b_1(i-1),b_2(j),b_3(m))}{2 \cdot 41} \;, \end{equation} \begin{equation} \frac{\partial}{\partial b_2}E_j(b_1(i),b_2(j),b_3(m)) = \frac{E_j(b_1(i),b_2(j+1),b_3(m)) - E_j(b_1(i),b_2(j-1),b_3(m))}{2 \cdot 41} \;, \end{equation} \begin{equation} \frac{\partial}{\partial b_3}E_j(b_1(i),b_2(j),b_3(m)) = \frac{E_j(b_1(i),b_2(j),b_3(m+1)) - E_j(b_1(i),b_2(j),b_3(m-1))}{2 \cdot 12} \;. \end{equation} The approximate gradient is calculated by the following Matlab function.

And the diagonal elemets of the conductivity tensor are calculated by the following Matlab function.

In this plot the calculated tensor elements are shown over a wide temperature range on a logarithmic scale. It can be seen that for low temperatures the conductivity goes to 0, as less and less states are contributing in the transport equation. The conductivity reaches its maximum between 500 K and 600 K and decays exponentially for higher temperature, due to the exponential law in the used approximation for the temperature dependent relaxation time. At room temperature (20°C) the calculated conductivity $\sigma_{22}$ is 4.1e+5 S/m (DFT) / 4.6e+5 S/m (tight binding), which fits in the order of magnitude to the value of 2e+5-3e+5 S/m given on Wikipedia , the calculated conductivity $\sigma_{33}$ is 7.7e+3 S/m (DFT) / 3.5e+3 S/m (tight binding), which is about a factor 10 higher than the value of 3.3e+2 S/m given on Wikipedia .

In this plot the the ratio $\sigma_{22}/\sigma_{33}$ is shown. It can be seen that the ratio increases for lower temperatures, since the in-plane conductivity $\sigma_{22}$ corresponds to the semimetal graphene, while the perpendicular conductivity $\sigma_{33}$ describes an insulator.

Here the calculated in-plane resistivity $\rho_{22}$ is plotted on the left hand side and compared with a graphic from the website of the European Carbon and Graphite Association on the right. When comparing the two figures, one can see that the shape of the curves is quiet similar and that the minimum is shifted by about 100 °C.