Transport Properties using the Boltzmann Equation

by Maria Troppenz & Santiago Rigamonti for exciting neon

(Jupyter notebook by Tejus Rohatgi, Mara Voiculescu & Martin Kuban)

---

Purpose: In this tutorial you will learn how to calculate electronic transport coefficients with exciting. The coefficients are obtained by solving the linearized Boltzmann equation in the constant relaxation time approximation. They are calculated for bulk silicon.

---

Table of Contents

0. Before Starting

1. Theoretical Background

2. Transport coefficients of Bulk Silicon

• Preparation of the calculation

• Calculation of the transport coefficients at room temperature

• Calculation of the transport coefficients at different temperatures

---

0. Before Starting

Read the following paragraphs before starting with the rest of this tutorial!

Before running any Jupyter tutorials, please refer to the 00_before_starting.md document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting's root directory:

source tools/excitingjupyter/venv/excitingvenv/bin/activate

As a first step, you may create a running directory for the notebook.

---

1. Theoretical Background

---

$\Rightarrow$ Transport coefficients by using the Boltzmann equation

The electronic transport coefficients, evaluated in this tutorial, are the electrical conductivity σ, the Seebeck coefficient *S*, and the electronic part of the thermal conductivity $\kappa_e$. They are related to the electronic structure of the material and directly affect the material's thermoelectric efficiency through the dimensionless figure of merit:

\begin{equation} \tag{1} ZT=σS^2\kappa T \end{equation}

Here, T is the temperature and $\kappa$ the total thermal conductivity, which is given by the sum of electronic contribution $\kappa_e$ and the lattice contribution $\kappa_l$, i.e. $\kappa=\kappa_e+\kappa_l$. The product $σS^2$ is also called the power factor.

The transport coefficients are obtained by solving the linearized Boltzmann's equation in the relaxation time approximation (for further information see Ref. [1]. The kernel of the transport coefficients is the transport distribution function,

\begin{equation} \tag{2} \Xi(\epsilon) =\frac{1}{m^2_e}\sum_{n\pmb{k}} \pmb{p}_{n\pmb{k}}\pmb{p}_{n\pmb{k}}\tau_{n\pmb{k}}δ(\epsilon−\epsilon_{n\pmb{k}}) \end{equation}

It is a tensor and depends on the carrier's relaxation time $\tau_{n\pmb{k}}$the dyadic product of the charge carrier's momentum matrix elements $\pmb{p}_{n\pmb{k}}$($m_e$ the electron mass). The sum is over all states with the band index $n$ and the wave vector $\pmb{k}$ in the first Brillouin Zone. In the constant relaxation time approximation, $\tau$ is independent from $n$ and $\pmb{k}$ and the transport distribution function is:

\begin{equation} \tag{3} \Xi(\epsilon) =\frac{\tau}{m^2_e}\sum_{n\pmb{k}} \pmb{p}_{n\pmb{k}}\pmb{p}_{n\pmb{k}}δ(\epsilon−\epsilon_{n\pmb{k}}) \end{equation}

Using Eq. (2), the transport coefficients can be calculated as follows:

\begin{equation} \tag{4} \pmb{\sigma}(T,\mu)= e^2 \int d\epsilon \left(-\frac{\partial{f_0 (T,\mu)}}{\partial{\epsilon}}\right) \Xi(\epsilon) \end{equation} \begin{equation} \tag{5} \pmb{S}(T,\mu)= \frac{ek_B}{\sigma} \int d\epsilon \left(-\frac{\partial{f_0 (T,\mu)}}{\partial{\epsilon}}\right) \Xi(\epsilon) \frac{\epsilon - \mu}{k_BT} \end{equation} \begin{equation} \tag{6} \pmb{\kappa_0}(T,\mu)=k_b^2T \int d\epsilon \left(-\frac{\partial{f_0 (T,\mu)}}{\partial{\epsilon}}\right) \Xi(\epsilon) \left(\frac{\epsilon - \mu}{k_BT}\right)^2 \end{equation}

Here, μ is the chemical potential, $e$ the electron charge, and $−\frac{\partial{f_0}}{\partial{\epsilon}}$ the derivative of the Fermi distribution function

\begin{equation} \tag{7} f_0(T,\mu)=\frac{1}{[exp(\epsilon−\mu)/k_BT]+1} \end{equation}

with respect to the energy $\epsilon$. $\kappa_0$ is the electronic part of the thermal conductivity for zero electrochemical potential gradient inside the sample, while the conventional electronic part of thermal conductivity, $\kappa_e$, is defined for zero electric current and is obtained by the relation $\kappa_e = \kappa_0 - T\sigma S^2$. The lattice part of the thermal conductivity still remains to be calculated by other methods.

Details about the method and its first ab initio implementation can be found in Refs. [2] and [3].

---

2. Transport coefficients of Bulk Silicon

To start the tutorial, create a new, empty directory named BE-Si and move into it.

%%bash
mkdir -p BE-Si

i) Preparation of calculation

Save the following lines as input.xml:

<input>
   <title>Bulk Silicon</title>
   <structure speciespath="$EXCITINGROOT/species/">
      <crystal scale="10.26">
         <basevect>0.0   0.5   0.5</basevect>
         <basevect>0.5   0.0   0.5</basevect>
         <basevect>0.5   0.5   0.0</basevect>
      </crystal>
      <species speciesfile="Si.xml">
         <atom coord="0.00  0.00  0.00"></atom>
         <atom coord="0.25  0.25  0.25"></atom>
      </species>
   </structure>
   <groundstate
      do="fromscratch"
      xctype="GGA_PBE_SOL"
      ngridk="4 4 4"
      rgkmax="8.0"
      gmaxvr="14.0"
      maxscl="100">
   </groundstate>
</input>

Set up $EXCITINGROOT to the correct root directory in speciespath and execute exciting_smp with the commands

%%bash
cd BE-Si
python3 -m excitingscripts.setup.excitingroot
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..

To gain a better understanding of the later calculated transport coefficients of silicon, you may want to look at the band structure. For this, please expand the link below and follow the instructions.

---

$\Rightarrow$ Calculation of the bandstructure of bulk silicon

Add the element properties in your input.xml with the following lines

<properties>
      <bandstructure>
         <plot1d>
            <path steps="300">
               <point coord="1.0     0.0     0.0" label="Gamma"/>
               <point coord="0.625   0.375   0.0" label="K"/>
               <point coord="0.5     0.5     0.0" label="X"/>
               <point coord="0.0     0.0     0.0" label="Gamma"/>
               <point coord="0.5     0.0     0.0" label="L"/>
            </path>
         </plot1d>
      </bandstructure> 
   </properties>

and execute again.

$ time 

After the calculation finished, you can visualize the band structure with the commands

$ python3 -m excitingscripts.plot.bandstructure -e -15 10

For a detailed explanation of calculating the band structure, see the tutorial Electronic band-structure and density of states.

---

ii) Calculation of the transport coefficients at room temperature

In the following, we calculate the transport coefficients with respect to the chemical potential at a constant temperature (here 300K). The equations of the transport coefficients (see Section 1) contain the derivative of the Fermi distribution function $f_0$. That means that only k-points with eigenenergies near the Fermi surface contribute to the coefficients. In order to achieve a sufficient representation of the energy bands in this area, a fine k-grid is needed.

To obtain eigenenergies for a fine k-grid, we perform a single self-consistency loop on top of the previous ground-state calculation by changing in input.xml

• the attribute do = "fromfile";
• the attribute ngridk = "40 40 40";
• the attribute maxscl = "1"

in the element groundstate (for more details see Input Reference.

The element groundstate in input.xml should now look like this:

...
   <groundstate
      do="fromfile"
      xctype="GGA_PBE_SOL"
      ngridk="40 40 40"
      rgkmax="8.0"
      gmaxvr="14.0"
      maxscl="1">
   </groundstate>
...

In order to calculate the transport coefficients, delete in your input.xml the "old" properties block that you used for calculating the electronic band-structure:

<properties>

      <bandstructure>
         ...
      </bandstructure>

   </properties>

and substitute it with

<properties>

      <momentummatrix/>
      <boltzequ
         energyReference="efermi"
         transportDfRange="-0.5 0.5"
         transportDfSpacing="0.001"
         temperatureRange="300 300"
         temperatureSpacing="1"
         chemicalPotentialRange="-0.05 0.05"
         chemicalPotentialSpacing="0.0001"
         evOutputEnergies="true"
         siOutputUnits="true">
         <etCoeffComponents>1 1</etCoeffComponents>
         <etCoeffComponents>2 2</etCoeffComponents>
         <etCoeffComponents>3 3</etCoeffComponents>     
      </boltzequ>

</properties>

With inserting the element momentummatrix, the momentum matrix elements are calculated. Inside the element boltzequ, you define the attributes needed for the calculation of the transport coefficients.

An explanation of the attributes can be found by expanding the link below or in Input Reference.

---

$\Rightarrow$ Table of the attributes used for the transport coefficients

Attribute	Element	Description
energyReference	boltzequ	If "efermi", the energy inputs are given with respect to the Fermi energy. Else "none".
useTransportDf	boltzequ	If "true" (default), the transport coefficients are calculated with the kernel function, the transport distribution function.
transportDfRange	boltzequ	Frequency (or energy) range for which the transport distribution function $\Xi(\epsilon)$ (See Eq. (2) in Section 1) is calculated.
transportDfSpacing	boltzequ	Spacing between the frequency (or energy) points used for the evaluation of the transport distribution function.
temperatureRange	boltzequ	Temperature range in Kelvin for which the transport coefficients are calculated.
temperatureSpacing	boltzequ	Spacing between temperature points in Kelvin.
chemicalPotentialRange	boltzequ	Chemical potential range for which the transport coefficients are calculated.
chemicalPotentialRange	boltzequ	Spacing between the chemical potential points.
useDopingConcentration	boltzequ	If "true", the chemical potential is adjusted to the doping concentration. Input parameters for the chemical potential are discarded.
dopingConcentration	boltzequ	Doping concentration in cm-3.
evOutputEnergies	boltzequ	If "true", output energies are in eV.
siOutputUnits	boltzequ	If "true", output transport coefficients are given in SI units: Seebeck coefficient in V/K, electrical conductivity over relaxation time in S/(m s), thermal conductivity over relaxation time in W/(m K s).

Using the element etCoeffComponents, you specify the tensor component for which the transport coefficients are calculated. The element etCoeffComponents needs, as input, an integer pair, e.g., 1 1 stands for the first diagonal component of the tensor ( 1 = x-direction, 2 = y-direction, and 3 = z-direction).

---

By executing exciting_smp

%%bash
cd BE-Si
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..

the Seebeck coefficient $S$, the electrical conductivity over the relaxation time $\sigma/\tau$, and the thermal conductivity over the relaxation time $\kappa_0/\tau$ are calculated for the temperature range and chemical potential range defined in input.xml. They are found in the files SEEBECK_COMP.OUT, ELECTCOND_COMP.OUT, and THERMALCOND_COMP.OUT, respectively, where, COMP denotes the tensor component specified by the integer pair given in condcomp. Furthermore, the figure of merit estimated by $Z=(\sigma S^2/\kappa_0)$ can be found in file Z_COMP.OUT.

Please note: $\sigma$ and $\kappa_0$ are not directly given, but the output files contain the coefficients $\sigma/\tau$ and $\kappa_0/\tau$. $\tau$ itself is constant (does not depend on electronic structure with $\pmb{k}$ and n), and can be multiplied afterwards.

Now, the transport coefficients are calculated and you can plot them by using excitingscripts.plot.files.

The Seebeck coefficient at the temperature of 300K can be visualized with the following commands:

%%bash
cd BE-Si
cp SEEBECK_11.OUT S11
cp SEEBECK_22.OUT S22
cp SEEBECK_33.OUT S33
python3 -m excitingscripts.plot.files -f S11 S22 S33  -cx 2  -cy 3  -ll '$S_{xx}$' '$S_{yy}$' '$S_{zz}$'  -x -1 1  -lx '$\mu$ [eV]'  -ly 'S [$\mu$V/K]'  -ys 1e6  -rp
cd ..

The output is saved in PLOT.png and looks as follows:

Assuming a relaxation time $\tau$ of 10 fs, the electrical conductivity can be obtained as:

%%bash
cd BE-Si
cp ELECTCOND_11.OUT E11
cp ELECTCOND_22.OUT E22
cp ELECTCOND_33.OUT E33
python3 -m excitingscripts.plot.files -f E11 E22 E33  -cx 2  -cy 3  -ll '$\sigma_{xx}$' '$\sigma_{yy}$' '$\sigma_{zz}$'  -x -0.5 0.5  -y 0 9000  -lx '$\mu$ [eV]'  -ly '$\sigma$ [S/cm]'  -ys 1e-16  -lp 'upper center'  -rp  -mty 5
cd ../

Output:

And the corresponding power factor looks as:

%%bash
cd BE-Si
paste ELECTCOND_11.OUT SEEBECK_11.OUT | awk '{ printf $1 "\t" $2 "\t" $3*$7*$7 "\n" }' > PF11
paste ELECTCOND_22.OUT SEEBECK_22.OUT | awk '{ printf $1 "\t" $2 "\t" $3*$7*$7 "\n" }' > PF22
paste ELECTCOND_33.OUT SEEBECK_33.OUT | awk '{ printf $1 "\t" $2 "\t" $3*$7*$7 "\n" }' > PF33
python3 -m excitingscripts.plot.files -f PF11 PF22 PF33  -cx 2  -cy 3  -ll '${\sigma S^2}_{xx}$' '${\sigma S^2}_{yy}$' '${\sigma S^2}_{zz}$'  -x -0.5 0.5  -y 0 20  -lx '$\mu$ [eV]'  -ly '$\sigma S^2$ [$\mu W\,cm^{-1}\,K^{-2}$]'  -ys 1e-10  -lp 'upper center'  -rp
cd ../

Output:

---

iii) Calculation of the transport coefficients at different temperatures

In the following, we look at the temperature dependence of transport coefficients at a constant chemical potential.

We can use the fine k-grid obtained in part ii) and skip the ground-state calculation by changing in input.xml

• The attribute do to "skip" in the element groundstate
• The momentum matrix elements does not need to be calculated again and, e.g., can be removed or commented.

The element properties of your input.xml can be changed to

<properties>

  <!--momentummatrix/-->
  <boltzequ
     energyReference="efermi"
     transportDfRange="-0.5 0.5"
     transportDfSpacing="0.001"
     temperatureRange="300 1200"
     temperatureSpacing="10"
     chemicalPotentialRange="-0.005 0.005"
     chemicalPotentialSpacing="0.005"
     evOutputEnergies="true"
     siOutputUnits="true">
     <etCoeffComponents>1 1</etCoeffComponents>
  </boltzequ>

</properties>

<properties>

      <momentummatrix/>
      <boltzequ
         energyReference="efermi"
         transportDfRange="-0.5 0.5"
         transportDfSpacing="0.001"
         temperatureRange="300 300"
         temperatureSpacing="1"
         chemicalPotentialRange="-0.05 0.05"
         chemicalPotentialSpacing="0.0001"
         evOutputEnergies="true"
         siOutputUnits="true">
         <etCoeffComponents>1 1</etCoeffComponents>
         <etCoeffComponents>2 2</etCoeffComponents>
         <etCoeffComponents>3 3</etCoeffComponents>     
      </boltzequ>

</properties>

The transport coefficients are calculated from 300 to 1200 K and for the chemical potentials $\mu = -0.136$ eV, 0 eV, and $0.136$ eV with respect to the Fermi energy. Here, we compute only the 1 1-component of the transport coefficient's tensor. With the modified input.xml, we can rerun exciting_smp

%%bash
cd BE-Si
time $EXCITINGROOT/bin/exciting_smp &
cd ../

The Seebeck coefficient for $\mu = -0.136$ eV and $\mu = 0.136$ eV can be visualized as:

%%bash
cd BE-Si
grep ' -0.136' SEEBECK_11.OUT > S11M
grep  ' 0.136' SEEBECK_11.OUT > S11P
python3 -m excitingscripts.plot.files -f S11P S11M  -cy 3  -ll '$S_{xx}$@$\mu$= 0.136eV' '$S_{xx}$@$\mu$=-0.136eV'  -x 300 1200  -lx 'T [K]'  -ly 'S [$\mu$V/K]'  -ys 1e6  -lp 5  -rp
cd ../

Output

Again, using τ = 10 fs, the electrical conductivity for these μ's is obtained as:

%%bash
cd BE-Si
grep ' -0.136' ELECTCOND_11.OUT > E11M
grep  ' 0.136' ELECTCOND_11.OUT > E11P
python3 -m excitingscripts.plot.files -f E11P E11M  -cy 3  -ll '$\sigma_{xx}$@$\mu$= 0.136eV' '$\sigma_{xx}$@$\mu$=-0.136eV'  -x 300 1200  -lx 'T [K]'  -ly '$\sigma$ [S/cm]'  -ys 1e-16  -rp
cd ../

Output:

And the corresponding power factor looks as:

%%bash
cd BE-Si
paste ELECTCOND_11.OUT SEEBECK_11.OUT | grep ' -0.136' | awk '{ printf $1 "\t" $2 "\t" $3*$7*$7 "\n" }' > PF11M
paste ELECTCOND_11.OUT SEEBECK_11.OUT | grep  ' 0.136' | awk '{ printf $1 "\t" $2 "\t" $3*$7*$7 "\n" }' > PF11P
python3 -m excitingscripts.plot.files -f PF11P PF11M  -cy 3  -ll '${\sigma S^2}_{xx}$@$\mu$= 0.136eV' '${\sigma S^2}_{xx}$@$\mu$=-0.136eV'  -x 300 1200  -lx 'T [K]'  -ly '$\sigma S^2$ [$\mu W\,cm^{-1}\,K^{-2}$]'  -ys 1e-10  -rp
cd ../

Output:

---

Bibliography

1. Method described in detail in: B. R. Nag, "Electron Transport in Compound Semiconductors" (Springer Verlag, Berlin, 1996).
2. G. D. Mahan and J. O. Sofo, "The best thermoelectric", Proc. Natl. Aad. Sci. USA 93, 7436 (1996).
3. T. J. Scheidemantel, C. Ambrosch-Draxl, T. Thonhauser, H. V. Badding, and J. O. Sofo, "Transport coefficients from first-rinciples calculations", Phys. Rev. B 68, 125210 (2003).

---