Phonon-limited mobility in graphene using the Boltzmann transport equation

Version: 2017.0

In this tutorial you will learn how to calculate the phonon-limited mobility in graphene. The mobility will be calculated using the Boltzmann transport equation (BTE) with the electronic structure, phonons and electron-phonon coupling calculated using density functional theory (DFT).

The mobility \(\mu\) will be calculated using two different methods to calculate the relaxation times \(\tau\) entering the BTE:

  1. Full angular (k,q)-dependence: in this method, the full dependency of \(\tau\) on the electron and phonon wave vectors \(\mathbf{k}\) and \(\mathbf{q}\) is taken into account, so that \(\tau = \tau(\mathbf{k},\mathbf{q})\). In the following, we will refer to this as the (k,q)-dependent method.
  2. Isotropic scattering rate: in this method, only the energy dependence of \(\tau\) is considered, so that \(\tau = \tau(E)\), and \(\tau(E)\) is assumed to vary isotropically in momentum-space. In the following, we will refer to this as the E-dependent method.

As it will be shown, the two methods give essentially the same results, but the second method is considerably faster than the first one.

In the theory section, you can read about the theoretical background. A more comprehensive description can also be found in the paper [GMSB16].

Geometry and electronic structure of graphene

In the builder_icon Builder, add a graphene configuration by clicking Add ‣ From Database and searching for ‘Graphene‘ in the structure database.

Next, you have to increase the vacuum gap above and below the graphene. Click on Bulk Tools ‣ Lattice parameters and set the lattice parameter along the C-direction to C = 20 Å.

Center the configuration by clicking on Coordinate Tools ‣ Center and then click on the sendto_icon button to send the structure to the script_generator_icon Script Generator.

To perform an accurate calculation of the relaxation times \(\tau\) and the mobility \(\mu\), the first important step is to optimize the lattice parameters A and B and calculate the electronic band structure. In order to do this, add the following blocks and set the associated parameters as follows:

  • calculator_icon New Calculator

    • Set the Occupation method to Methfessel-Paxton

    • Set the Density mesh cutoff to 90 Ha

    • Set the Exchange correlation to LDA

    • Set the k-point sampling to:

      • \(\mathrm{k_A}\) = 33
      • \(\mathrm{k_B}\) = 33
      • \(\mathrm{k_C}\) = 1
    ../../_images/unitcell_calculator1.png
  • optimization_icon Optimization ‣ OptimizeGeometry

    • Set the Force tolerance to 0.001 eV/Å
    • Untick Constrain lattice vectors in the \(x\) and \(y\) directions.
    ../../_images/optimizegeometry.png
  • analysis_icon Analysis ‣ Bandstructure

    • Set the number of Points per segment to 100
    • Set the Brillouin zone route to G, K, M, G

In the main panel of the script_generator_icon Script Generator, set the Default output file to ‘Graphene_relax.hdf5‘, send the script to the job_manager_icon Job manager, save it as ‘Graphene_relax.py‘, and click on the start_icon button to run the calculation.

When the calculation is done, click on the labfloor_bandstructure_icon Bandstructure object contained in the file ‘Graphene_relax.hdf5‘ in the LabFloor and use the Bandstructure Analyzer to visualize the band structure.

../../_images/bandstructure_analyzer.png

By placing the mouse cursor on top of a band, information about the band is shown. You see that the valence band (highlighted in orange in the figure above) is number 3 and the conduction band is number 4. These are the two electronic bands relevant for the calculation of the mobility, and in the following we will concentrate on these two bands.

Phonons in Graphene

The next step is to calculate the dynamical matrix of graphene. In order to test the quality of the result you will also calculate the phonon band structure, which is based on the calculated dynamical matrix.

In the LabFloor, drag and drop the relaxed configuration of graphene (gID001) contained in the file ‘Graphene_relax.hdf5‘ to the script_generator_icon Script generator, and modify the script as follows:

  • In the calculator_icon New Calculator block, set the k-point sampling to:

    • \(\mathrm{k_A}\) = 3
    • \(\mathrm{k_B}\) = 3
    • \(\mathrm{k_C}\) = 1
  • Add an analysis_icon Analysis ‣ DynamicalMatrix block and modify the following settings:

    • Set Repetitions to Custom
    • Set the Number of repetitions to:
      • \(\mathrm{n_A}\) = 11
      • \(\mathrm{n_B}\) = 11
      • \(\mathrm{n_C}\) = 1
../../_images/dynamicalmatrix.png

Note

To ensure an equivalent \(\mathbf{k}\)-points sampling between the calculation of the graphene unit cell, used for the electronic structure, and the graphene supercell, used for the phonons calculation, it is important to set \(\mathrm{k_i^{unit cell}} = \mathrm{k_i^{supercell}} \times \mathrm{N_i^{supercell}}\), for \(\mathrm{i} = \mathrm{A},\mathrm{B},\mathrm{C}\).

Tip

Due to the large dimensions of the \(11 \times 11\) graphene super cell (242 atoms), the calculation of the DynamicalMatrix object is rather time consuming. However, the calculation can be parallelized over the atomic displacements. Since there are two atoms in the graphene unit cell, and each is displaced in the \(x\), \(y\) and \(z\) directions, there are in total 6 calculations to be performed. Maximum efficiency is obtained when the number of calculations times the value of the parameter Processes per displacement matches the total number of cores used for the calculations. In the present case, the calculation takes only 30 minutes if Processes per displacement = 4 and the calculation is run on 24 cores.

  • Add an analysis_icon Analysis ‣ PhononBandstructure block and set the following parameters:

    • Set the number of Points per segment to 100
    • Set the Brillouin zone route to G, M, K, G

Finally, in the main panel of the script_generator_icon Script Generator set the Default output file to ‘Graphene_dynmat.hdf5‘, send the script to the job_manager_icon Job manager, save it as ‘Graphene_dynmat.py‘ and click on the start_icon button to run the calculation.

When the calculation is done, go back to the in the LabFloor, and inspect the labfloor_phononbandstructure_icon PhononBandstructure object contained in the file ‘Graphene_dynmat.hdf5‘ using Compare data. The calculated phonon band structure should match with that shown in the figure below.

../../_images/phonon_band.png

Fig. 31 Phonon bandstructure. There are three acoustic modes: The lowest, out-of-plane mode (ZA) has a \(\propto \mathbf{q}^2\) for small \(\mathbf{q}\), while the two next have linear \(\mathbf{q}\)-dependence with a constant velocity for small \(\mathbf{q}\). The transverse acoustic (TA) mode is lower in energy than the longitudinal acoustic (LA) mode.

Mobility of graphene

In the following, the procedure to calculate the electron mobility in graphene will be described. Provided that one has already calculated the electronic structure and dynamical matrix of the system, the following three steps are necessary to evaluate the mobility:

  1. Calculation of the Hamiltonian derivatives
  2. Calcualtion of the Electron-phonon couplings
  3. Calculation of the Mobility

The two methods to calculate \(\mu\) described above differ in the way in which steps 2 and 3 are carried out. These two steps will therefore be described separately for each method. The procedure for the (k, q)-dependent method is outlined in sections 1, 2A and 3A, whereas the procedure for the energy-dependent method is outlines in in sections 1, 2B and 3B

1. Hamiltonian derivatives

In order to calculate the electron-phonon coupling matrix, it is necessary to calculate the derivative of the Hamiltonian \(\partial \hat{H}/\partial R_{i,\alpha}\) with respect to the coordinate of the \(i\)-th atom along the Cartesian direction \(\alpha\).

In the LabFloor, drag and drop the relaxed configuration of graphene (gID001) contained in the file ‘Graphene_relax.hdf5‘ to the script_generator_icon Script generator, and modify the script as follows:

  • In the calculator_icon New Calculator block, set the k-point sampling to:

    • \(\mathrm{k_A}\) = 3
    • \(\mathrm{k_B}\) = 3
    • \(\mathrm{k_C}\) = 1
  • Add an analysis_icon Analysis ‣ HamiltonianDerivatives block and set the Number of repetitions to:

    • \(\mathrm{n_A}\) = 11
    • \(\mathrm{n_B}\) = 11
    • \(\mathrm{n_C}\) = 1
../../_images/hamiltonianderivatives.png

Warning

It is important that the number of repetitions matches that used for the calculation of the dynamical matrix.

In the main panel of the script_generator_icon Script Generator, set the Default output file to ‘Graphene_dHdR.hdf5‘, send the script to the job_manager_icon Job manager, save it as ‘Graphene_dHdR.py‘ and click on the start_icon button to run the calculation.

Note

As for the DynamicalMatrix, the calculation of the HamiltonianDerivatives is rather time consuming, but can be parallelized over the atomic displacements. In the present case, the calculation takes around 1 hour if it is run on 24 cores.

2A. Electron-Phonon couplings: (k,q)-dependent method

In order to calculate the lifetimes \(\tau(\mathbf{k},\mathbf{q})\) and the mobility \(\mu\) we need to calculate the electron-phonon coupling matrix on a fine grid of \(\mathbf{k}\)- and \(\mathbf{q}\)-points.

In the LabFloor, select the relaxed configuration (gID001) object in the file ‘Graphene_relax.hdf5‘, and drag and drop it to the script_generator_icon Script generator. Add the following block, and modify the parameters as follows:

  • analysis_icon Analysis ‣ ElectronPhononCoupling
    • In the k-points sampling section, set the k grid to Custom (Cartesian), and the sampling parameters min, max and number as shown in the figure below
    • In the q-points sampling section, set the q grid to Custom (Cartesian), and the sampling parameters min, max and number as shown in the figure below
    • Set Electron bands to 3, 4
../../_images/electron_phonon_fine.png

Note

The sampling ranges are chosen so that only the relevant regions of the Brillouin zone in \(\mathbf{k}\)-space and \(\mathbf{q}\)-space are sampled. For the former, this region is that in the vicinity of the \(\mathrm{K}\)-point, which is located at [1.7063, 0.0, 0.0]. For the latter, this region is that in the vicinity of the \(\mathrm{\Gamma}\)-point at [0.0, 0.0, 0.0].

Warning

DynamicalMatrix and HamiltonianDerivatives are automatically also added, when you add an ElectronPhononCoupling object. However, in the present case, we do not need to calculate them again, as they have been calculated previously.

In the main panel of the script_generator_icon Script Generator, set the Default output file to ‘Graphene_eph.hdf5‘, transfer the script to the editor_icon Editor and replace the Dynamical Matrix and Hamiltonian Derivatives blocks with the following lines:

# -------------------------------------------------------------
# Dynamical matrix
# -------------------------------------------------------------
dynamical_matrix = nlread('Graphene_relax.hdf5', DynamicalMatrix)[-1]

# -------------------------------------------------------------
# Hamiltonian derivatives
# -------------------------------------------------------------
hamiltonian_derivatives = nlread('Graphene_dHdR.hdf5', HamiltonianDerivatives)[-1]

Furthermore, in the Electron Phonon Coupling block modify the parameters energy_tolerance and initial_state_energy_range as follows:

 electron_phonon_coupling = ElectronPhononCoupling(
     configuration=bulk_configuration,
     dynamical_matrix=dynamical_matrix,
     hamiltonian_derivatives=hamiltonian_derivatives,
     kpoints_cartesian=kpoints,
     qpoints_cartesian=qpoints,
     electron_bands=[3,4],
     phonon_modes=All,
     energy_tolerance=0.01*eV,
     initial_state_energy_range=[-0.5,0.5]*eV,
     store_dense_coupling_matrices=False,
     )

Note

Sampling the \(\mathbf{k}\)- and \(\mathbf{q}\)-space using dense meshes means that the calculation might become very time consuming. The two parameters energy_tolerance and initial_state_energy_range allow to consider only the relevant electron-phonon coupling element and therefore decrease considerably the computational time and the memory requirements. In the present case, the calculation takes 1.5 hours on 24 cores. More information about these two parameters can be found in the section “Usage examples” in the manual page ElectronPhononCoupling.

Send the script to the job_manager_icon Job manager, save the input file as ‘Graphene_eph.py‘ and click on the start_icon button to run the calculation.

3A. Electron mobility: (k,q)-dependent method

In the LabFloor, select again the relaxed configuration (gID001) object in the file Graphene_relax.hdf5, and drag and drop it to the script_generator_icon Script generator. Add the following block, and modify the parameters as follows:

  • analysis_icon Analysis ‣ Mobility
    • Set the Method to Full angular (k,q)-dependence
    • Set the Fermi shift to 0.13 eV
    • Untick Calculate Hall coefficients
../../_images/mobility_1.png

Note

The Fermi shift of 0.13 eV corresponds to a carrier concentration of \(n = 10^{12} \mathrm{cm^{-2}}\)

In the main panel of the script_generator_icon Script Generator, set the Default output file to ‘Graphene_mu-full.hdf5‘, send the script to the job_manager_icon Job manager, save it as ‘Graphene_mu-full.py‘ and click on the start_icon button to run the calculation.

Once the calculation is done, open the file ‘Graphene_mu-full.hdf5‘ in the LabFloor, select the labfloor_mobility_icon Mobility object, and click on Text Representation.

 +------------------------------------------------------------------------------+
 | Mobility Report                                                              |
 | ---------------------------------------------------------------------------- |
 | Input parameters:                                                            |
 |     Temperature       = 300.00 K                                             |
 |     Fermi level shift = 0.13 eV                                              |
 |     Energy broadening = 0.0030 eV                                            |
 |     q-grid refinement = 1                                                    |
 |                                                                              |
 +------------------------------------------------------------------------------+
 | Trace of linear responce tensors:                                            |
 +------------------------------------------------------------------------------+
 |                                                                              |
 | Electrons:                                                                   |
 |                                                                              |
 |     Mobility                 =  2.56e+05 cm^2/(V*s)                          |
 |     Conductivity             =  7.71e+06 S/m                                 |
 |     Seebeck coefficient      = -6.11e-05 V/K                                 |
 |     Thermal conductivity     =  4.46e+01 W/(m*K)                             |
 |     Carrier density (2D, xy) =  3.76e+11 cm^-2                               |
 |                                                                              |
 |                                                                              |
 | Holes:                                                                       |
 |                                                                              |
 |     Mobility                 =  2.56e+06 cm^2/(V*s)                          |
 |     Conductivity             =  4.27e+04 S/m                                 |
 |     Seebeck coefficient      =  1.81e-01 V/K                                 |
 |     Thermal conductivity     =  4.02e+00 W/(m*K)                             |
 |     Carrier density (2D, xy) =  2.09e+08 cm^-2                               |
 |                                                                              |
 +------------------------------------------------------------------------------+

The calculated electron mobility, highlighted in yellow, is \(2.56 \cdot 10^{5}\ \mathrm{cm^{2} V^{-1} s^{-1}}\), and agrees well with previously reported data at room temperature and \(n = 10^{12} \mathrm{cm^{-2}}\) [KTJ12b] .

2B. Electron-Phonon couplings: energy-dependent method

We will now calculate the electron-phonon couplings to be used for the calculation of the energy-dependent relaxation times \(\tau(E)\). We will assume that the relaxation times vary isotropically in \(\mathbf{k}\)-space. This means that it will be sufficient to evaluate the electron-phonon coupling matrix for a number of \(\mathbf{k}\)-points along a line outward from one of the Dirac (K) points, thereby reducing significantly the computational workload of the calculation.

In the LabFloor, drag and drop the relaxed graphene configuration (gID001) in the file ‘Graphene_relax.hdf5‘ to the script_generator_icon Script generator. Add the following block, and modify the parameters as follows:

  • analysis_icon Analysis ‣ ElectronPhononCoupling

    • In the k-points sampling, set the k grid to Custom (Cartesian), and the sampling parameters min, max and number as shown in the figure below

      Note

      As you can see from the figure, the \(\mathbf{k}\)-space is sampled only along a line!

    • In the q-points sampling, set the k grid to Custom (Cartesian), and the sampling parameters min, max and number as shown in the figure below

    • Set Electron bands to 3, 4

../../_images/epc_k_line.png

In the main panel of the script_generator_icon Script Generator, set the Default output file to ‘Graphene_eph-line.hdf5‘, send the script to the editor_icon Editor, and replace the Dynamical Matrix and Hamiltonian Derivatives blocks with the following lines:

# -------------------------------------------------------------
# Dynamical matrix
# -------------------------------------------------------------
dynamical_matrix = nlread('Graphene_dynmat.hdf5', DynamicalMatrix)[-1]

# -------------------------------------------------------------
# Hamiltonian derivatives
# -------------------------------------------------------------
hamiltonian_derivatives = nlread('Graphene_dHdR.hdf5', HamiltonianDerivatives)[-1]

Furthermore, in the Electron Phonon Coupling block modify the parameters energy_tolerance and initial_state_energy_range as follows:

 electron_phonon_coupling = ElectronPhononCoupling(
     configuration=bulk_configuration,
     dynamical_matrix=dynamical_matrix,
     hamiltonian_derivatives=hamiltonian_derivatives,
     kpoints_cartesian=kpoints,
     qpoints_cartesian=qpoints,
     electron_bands=[3,4],
     phonon_modes=All,
     energy_tolerance=0.01*eV,
     initial_state_energy_range=[-0.5,0.5]*eV,
     store_dense_coupling_matrices=False,
     )

Send the script to the job_manager_icon Job manager, save it as ‘Graphene_eph-line.py‘ and click on the start_icon button to run the calculation.

Note

Since here we sample the \(\mathbf{k}\)-space only along a line, the calcualtion is much faster than the one described in section 2A, and takes only 5 minutes to complete using 24 cores!

3B. Electron mobility: energy-dependent method

We will now use a two-step procedure to evaluate the room-temperature mobility \(\mu\) based on the energy-dependent relaxation times \(\tau(E)\):

  • In the first step, the \(\mathbf{k}\)- and \(\mathbf{q}\)-dependent relaxation times \(\tau(\mathbf{k},\mathbf{q})\) are evaluated along the line outward from the K-point.
  • In the second step, the values of \(\tau(\mathbf{k},\mathbf{q})\) are averaged in \(\mathbf{k}\)-space to obtain the values of \(\tau(E)\), which are then be used to calculate the mobility \(\mu\).

In the LabFloor, drag and drop the relaxed configuration (gID001) object in the file ‘Graphene_relax.hdf5‘ to the script_generator_icon Script generator. Add the following block, and modify the parameters as follows:

  • analysis_icon Analysis ‣ Mobility
    • Set the Method to Full angular (k,q)-dependence
    • Set the Fermi shift to 0.13 eV
    • Untick Calculate Hall coefficients
../../_images/mobility_2a.png

In the main panel of the script_generator_icon Script Generator, set the Default output file to ‘Graphene_mu-line-full.hdf5‘, send the script to the editor_icon Editor, and replace the Dynamical Matrix, the Hamiltonian Derivatives and the Electron Phonon Coupling blocks with the following lines:

# -------------------------------------------------------------
# Electron Phonon Coupling
# -------------------------------------------------------------
electron_phonon_coupling = nlread('Graphene_eph-line.hdf5',ElectronPhononCoupling)[-1]

Then, send the script to the job_manager_icon Job manager,save it as ‘Graphene_mu-line-full.py‘ and click on the start_icon button to run the calculation.

Once the calculation is done, drag and drop again the relaxed configuration to the script_generator_icon Script generator. Add the following block, and modify the parameters as follows:

  • analysis_icon Analysis ‣ Mobility
    • Set the Method to Isotropic scattering rate
    • Set the Fermi shift to 0.13 eV
    • Untick Calculate Hall coefficients
    • Set the following values of the Energy Range:
      • \(\mathrm{E_0}\) = -0.24 eV
      • \(\mathrm{E_0}\) = 0.24 eV
      • Points = 100
    • Set the k-point Sampling to \(99 \times 99 \times 1\)
../../_images/mobility_2b.png

In the main panel of the script_generator_icon Script Generator, set the Default output file to ‘Graphene_mu-line-iso.hdf5‘, send the script to the editor_icon Editor, remove the Dynamical Matrix, Hamiltonian Derivatives and Electron Phonon Couplings blocks, and modifiy the Mobility block as follows:

# -------------------------------------------------------------
# Mobility
# -------------------------------------------------------------
mobility_full = nlread('Graphene_mu-line-full.hdf5', Mobility)[-1]

# Generate energy-dependent inverse relaxation times
inverse_relaxation_time=mobility_full.generateEnergyDependentInverseRelaxationTime(
    energies=numpy.linspace(-0.24, 0.24, 100)*eV
    )

kpoint_grid = MonkhorstPackGrid(
    na=99,
    nb=99,
    )

mobility = Mobility(
    configuration=bulk_configuration,
    method=Isotropic,
    energies=numpy.linspace(-0.24, 0.24, 100)*eV,
    inverse_relaxation_time=inverse_relaxation_time,
    electron_bands=[3,4],
    kpoints=kpoint_grid,
    temperature=300*Kelvin,
    fermi_shift=0.13*eV,
    energy_broadening=0.003*eV,
    calculate_hall_coefficients=False,
    )
nlsave('Graphene_mu-line-iso.hdf5', mobility)

Send the script to the job_manager_icon Job manager, save it as ‘Graphene_mu-line-iso.py‘ and click on the start_icon button to run the calculation.

Once the calculation is done, open the file ‘Graphene_mu-line-iso.hdf5‘, select the labfloor_mobility_icon Mobility object in the LabFloor and click on Text Representation.

 +------------------------------------------------------------------------------+
 | Mobility Report                                                              |
 | ---------------------------------------------------------------------------- |
 | Input parameters:                                                            |
 |     Temperature       = 300.00 K                                             |
 |     Fermi level shift = 0.13 eV                                              |
 |     Energy broadening = 0.0030 eV                                            |
 |     q-grid refinement = 1                                                    |
 |                                                                              |
 +------------------------------------------------------------------------------+
 | Trace of linear responce tensors:                                            |
 +------------------------------------------------------------------------------+
 |                                                                              |
 | Electrons:                                                                   |
 |                                                                              |
 |     Mobility                 =  2.58e+05 cm^2/(V*s)                          |
 |     Conductivity             =  1.79e+07 S/m                                 |
 |     Seebeck coefficient      =  7.91e-05 V/K                                 |
 |     Thermal conductivity     =  1.16e+02 W/(m*K)                             |
 |     Carrier density (2D, xy) =  8.63e+11 cm^-2                               |
 |                                                                              |
 |                                                                              |
 | Holes:                                                                       |
 |                                                                              |
 |     Mobility                 =  1.68e+06 cm^2/(V*s)                          |
 |     Conductivity             =  4.94e+04 S/m                                 |
 |     Seebeck coefficient      = -1.25e-03 V/K                                 |
 |     Thermal conductivity     =  5.98e+00 W/(m*K)                             |
 |     Carrier density (2D, xy) =  3.66e+08 cm^-2                               |
 |                                                                              |
 +------------------------------------------------------------------------------+

The calculated value for the electron mobility, highlighted in yellow above, is \(2.58 \cdot 10^{5}\ \mathrm{cm^{2} V^{-1} s^{-1}}\), in excellent agreement with the value obtained by using the (k,q)-dependent method.

Temperature dependence of the mobility: (k,q)-dependent method vs. energy-dependent method

A more stringent test for the reliability of the energy-dependent method is to calculate the temperature dependence of the mobility in an energy range around room temperature.

  • In the (k,q)-dependent method, this can be done by simply modifying the value of the target temperature in the analysis_icon Mobility analysis object.
  • In the energy-dependent method, the two-step procedure must be repeated for each temperature, and the target temperature has to be set both when calculating the values of the \(\mathbf{k}\)- and \(\mathbf{q}\)-dependent relaxation times \(\tau(\mathbf{k},\mathbf{q})\) along the line, and when calculating the energy-dependent relaxation times \(\tau(E)\).

As it can be seen from the figure below, where the \(T\)-dependency of \(\mu\) is calculated in the temperature range \(100 \mathrm{K} \leq T \leq 300 \mathrm{K}\), both methods reproduce well the expected \(1/T\) behaviour of \(\mu\).

../../_images/mu_vs_T.png

Theory section

The mobility \(\mu\) is related to the conductivity \(\sigma\) via

(1)\[\mu = \frac{\sigma}{ne},\]

where \(n\) is the carrier density, \(e\) is the electronic charge.

In this tutorial, you will calculate the conductivity using the semi-classical Boltzmann transport equation (BTE). In the case of zero magnetic field and a time-independent electric fields in the steady state limit the BTE simplifies to:

(2)\[\frac{q\mathbf{E}}{\hbar} \cdot \nabla_{\mathbf{k}} f_{\mathbf{k} n} \approx \frac{q\mathbf{E}}{\hbar} \cdot \nabla_{\mathbf{k}} f^0_{\mathbf{k} n} = \left.\frac{\partial f_{\mathbf{k}n}}{\partial t}\right|_{coll}\]

The right hand side (RHS), often called the collision integral, describes different sources of scattering and dissipation that drives the system towards steady state. The left hand-side is approximated to linear order in the electric field by changing to the equilibrium distribution.

The electron-phonon scattering is described by the collision integral. Most commonly this is approximated by a relaxation time:

\[\left.\frac{\partial f_{\mathbf{k}n}}{\partial t}\right|_{coll} = \frac{\delta f_{\mathbf{k}n}}{\tau_{\mathbf{k}n}}\]

describing quasielastic scattering on acoustic phonons (relaxation time approximation). We let \(\lambda\) label the phonon modes, \(n\) the electronic bands, \(\mathbf{k}\) the electron momentum. The transition rate from a state \(|\mathbf{k}n\rangle\) to \(|\mathbf{k}'n'\rangle\) is obtained from Fermi’s golden rule (FGR):

\[\begin{split}P_{\mathbf{k}\mathbf{k'}}^{\lambda nn'} &= \frac{2\pi}{\hbar} |g_{\mathbf{k}\mathbf{k'}}^{\lambda n n'}|^2 [ n_{\mathbf{q}}^{\lambda} \delta \left(\epsilon_{\mathbf{k}'n'}-\epsilon_{\mathbf{k}n}-\hbar \omega_{q \lambda} \right) \delta_{\mathbf{k}',\mathbf{k}+\mathbf{q}} \\[.5cm] &+ (n_{\mathbf{-q}}^{\lambda}+1) \delta \left( \epsilon_{\mathbf{k}'n'}-\epsilon_{\mathbf{k} n}+\hbar \omega_{-q \lambda} \right) \delta_{\mathbf{k}',\mathbf{k}-\mathbf{q}} ],\end{split}\]

where the first and last terms in the square brackets describes phonon absorption and emission, respectively. The electron-phonon coupling matrix elements \(|g_{\mathbf{k}\mathbf{k'}}^{\lambda n n'}|\) are in ATK calculated using the ElectronPhononCoupling analysis module.

The general electron-phonon collision integral is given by

(3)\[\begin{split}\left.\frac{\partial f_{\mathbf{k}n}}{\partial t}\right|_{coll} &= - \sum_{\mathbf{k}'n'} \left[f_{\mathbf{k}n}\left(1-f_{\mathbf{k}'n'}\right)P_{\mathbf{k}\mathbf{k}'}^{nn'}-f_{\mathbf{k}'n'}\left(1-f_{\mathbf{k}n}\right)P_{\mathbf{k}'\mathbf{k}}^{n'n}\right]\end{split}\]

Relaxation time approximation

In the relaxation-time approximation (RTA) one assumes a special simplified form of the RHS/collision integral:

(4)\[\begin{split}\left.\frac{\partial f_{\mathbf{k}n}}{\partial t}\right|_{coll} &=& -\frac{\delta f_{\mathbf{k}n}}{\tau_{\mathbf{k} n}}\end{split}\]

where \(\delta f_{\mathbf{k}n} = f_{\mathbf{k}n}-f^0_{\mathbf{k}n}\). The linearized BTE becomes:

\[q \mathbf{E} \cdot \mathbf{v}_{\mathbf{k} n} \frac{\partial f^0_{\mathbf{k}n}}{\partial \epsilon_{\mathbf{k}n}} = -\frac{\delta f_{\mathbf{k}n}}{\tau_{\mathbf{k} n}}.\]

Hereby we find the solution:

\[\delta f_{\mathbf{k}n} = -q \mathbf{E} \cdot \mathbf{v}_{\mathbf{k} n} \frac{\partial f^0_{\mathbf{k}n}}{\partial \epsilon_{\mathbf{k}n}} \tau_{\mathbf{k} n}.\]

The relaxation-time can be evaluated from the general collision integral, (3). However, the specific form in eqn. (4) only applies in the limit of quasielastic scattering (\(\omega_{q \lambda}\rightarrow 0\)), in which case the collision integral, in eqn. (3), simplifies significantly:

\[\left.\frac{\partial f_{\mathbf{k}n}}{\partial t}\right|_{coll}^{RTA} = \sum_{\mathbf{k}'n'} \left(f_{\mathbf{k}'n'}-f_{\mathbf{k}n}\right) P_{\mathbf{k}\mathbf{k'}}^{nn'} \equiv -\frac{\delta f_{\mathbf{k}n}}{\tau_{\mathbf{k} n}},\]

since \(P_{\mathbf{k}\mathbf{k'}}^{nn'} = P_{\mathbf{k}'\mathbf{k}}^{n'n}\) in this limit.

We then obtain the expression for the scattering rate or inverse relaxation-time:

(5)\[\frac{1}{\tau_{\mathbf{k} n}} = \sum_{\mathbf{k}'n'} \left(1-\frac{\delta f_{\mathbf{k}'n'}}{\delta f_{\mathbf{k}n}}\right) P_{\mathbf{k}\mathbf{k'}}^{nn'}\]

where we applied the detailed balance equation in this limit, \(P_{\mathbf{k}\mathbf{k'}}^{nn'}(f^0_{\mathbf{k}'n'}-f^0_{\mathbf{k}n})=0\), and assumed isotropic scattering.

To evaluate the relaxation time we introduce the scattering angle

\[cos(\theta_{\mathbf{k}\mathbf{k'}}) = \frac{\mathbf{v}_{\mathbf{k}'n'} \cdot \mathbf{v}_{\mathbf{k}n}}{|\mathbf{v}_{\mathbf{k}'n'}| |\mathbf{v}_{\mathbf{k}n}|}\]

where \(\mathbf{v}_{\mathbf{k}n}\) is the electron velocity. Following [GMSB16] we obtain

\[\frac{1}{\tau_{\mathbf{k} n}} = =\sum_{\mathbf{k}'n'} \frac{(1-f^0_{\mathbf{k}'n'})}{(1-f^0_{\mathbf{k}n})}\left(1-cos(\theta_{\mathbf{k}\mathbf{k'}})\right) P_{\mathbf{k}\mathbf{k'}}^{nn'}\]

Mobility

Once knowing the relaxation times one obtain the mobility as:

(6)\[\mu = -2q \frac{\sum_{\mathbf{k}n} v_{\mathbf{k} n}^2 \frac{\partial f^0_{\mathbf{k}n}}{\partial \epsilon_{\mathbf{k}n}} \tau_{\mathbf{k} n}}{\sum_{\mathbf{k}n}f^0_{\mathbf{k}n}}\]

where the factor 2 accounts for spin degeneracy.

References

[GMSB16](1, 2) T. Gunst, T. Markussen, K. Stokbro, and M. Brandbyge. First-principles method for electron-phonon coupling and electron mobility: Applications to two-dimensional materials. Phys. Rev. B, 93:035414, Jan 2016. doi:10.1103/PhysRevB.93.035414.
[KTJ12b]Kristen Kaasbjerg, Kristian S Thygesen, and Karsten W Jacobsen. Unraveling the acoustic electron-phonon interaction in graphene. Physical Review B, 85(16):165440, 2012.