Phononlimited mobility in graphene using the Boltzmann transport equation¶
Version: 2017.0
In this tutorial you will learn how to calculate the phononlimited mobility in graphene. The mobility will be calculated using the Boltzmann transport equation (BTE) with the electronic structure, phonons and electronphonon 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:
 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.
 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 momentumspace. In the following, we will refer to this as the Edependent 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, add a graphene configuration by clicking and searching for ‘Graphene‘ in the structure database.
Next, you have to increase the vacuum gap above and below the graphene. Click on
and set the lattice parameter along the Cdirection to C = 20 Å.Center the configuration by clicking on button to send the structure to the Script Generator.
and then click on theTo 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:

Set the Occupation method to MethfesselPaxton
Set the Density mesh cutoff to 90 Ha
Set the Exchange correlation to LDA
Set the kpoint sampling to:
 \(\mathrm{k_A}\) = 33
 \(\mathrm{k_B}\) = 33
 \(\mathrm{k_C}\) = 1

 Set the Force tolerance to 0.001 eV/Å
 Untick Constrain lattice vectors in the \(x\) and \(y\) directions.

 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, set the Default output file to ‘Graphene_relax.hdf5‘, send the script to the Job manager, save it as ‘Graphene_relax.py‘, and click on the button to run the calculation.
When the calculation is done, click on the Bandstructure object contained in the file ‘Graphene_relax.hdf5‘ in the LabFloor and use the Bandstructure Analyzer to visualize the band structure.
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, and modify the script as follows:
In the block, set the kpoint sampling to:
 \(\mathrm{k_A}\) = 3
 \(\mathrm{k_B}\) = 3
 \(\mathrm{k_C}\) = 1
Add an 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
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 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 set the Default output file to ‘Graphene_dynmat.hdf5‘, send the script to the Job manager, save it as ‘Graphene_dynmat.py‘ and click on the button to run the calculation.
When the calculation is done, go back to the in the LabFloor, and inspect the 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.
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:
 Calculation of the Hamiltonian derivatives
 Calcualtion of the Electronphonon couplings
 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 energydependent method is outlines in in sections 1, 2B and 3B
1. Hamiltonian derivatives¶
In order to calculate the electronphonon 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, and modify the script as follows:
In the block, set the kpoint sampling to:
 \(\mathrm{k_A}\) = 3
 \(\mathrm{k_B}\) = 3
 \(\mathrm{k_C}\) = 1
Add an block and set the Number of repetitions to:
 \(\mathrm{n_A}\) = 11
 \(\mathrm{n_B}\) = 11
 \(\mathrm{n_C}\) = 1
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, set the Default output file to ‘Graphene_dHdR.hdf5‘, send the script to the Job manager, save it as ‘Graphene_dHdR.py‘ and click on the 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. ElectronPhonon 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 electronphonon 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. Add the following block, and modify the parameters as follows:

 In the kpoints 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 qpoints 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
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, set the Default output file to ‘Graphene_eph.hdf5‘, transfer the script to the 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 electronphonon 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, save the input file as ‘Graphene_eph.py‘ and click on the 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. Add the following block, and modify the parameters as follows:

 Set the Method to Full angular (k,q)dependence
 Set the Fermi shift to 0.13 eV
 Untick Calculate Hall coefficients
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, set the Default output file to ‘Graphene_mufull.hdf5‘, send the script to the Job manager, save it as ‘Graphene_mufull.py‘ and click on the button to run the calculation.
Once the calculation is done, open the file ‘Graphene_mufull.hdf5‘ in the LabFloor, select the 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 
 qgrid 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.11e05 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.81e01 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. ElectronPhonon couplings: energydependent method¶
We will now calculate the electronphonon couplings to be used for the calculation of the energydependent 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 electronphonon 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. Add the following block, and modify the parameters as follows:

In the kpoints 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 qpoints 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
In the main panel of the Script Generator, set the Default output file to ‘Graphene_ephline.hdf5‘, send the script to the 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, save it as ‘Graphene_ephline.py‘ and click on the 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: energydependent method¶
We will now use a twostep procedure to evaluate the roomtemperature mobility \(\mu\) based on the energydependent 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 Kpoint.
 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. Add the following block, and modify the parameters as follows:

 Set the Method to Full angular (k,q)dependence
 Set the Fermi shift to 0.13 eV
 Untick Calculate Hall coefficients
In the main panel of the Script Generator, set the Default output file to ‘Graphene_mulinefull.hdf5‘, send the script to the 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_ephline.hdf5',ElectronPhononCoupling)[1]
Then, send the script to the Job manager,save it as ‘Graphene_mulinefull.py‘ and click on the button to run the calculation.
Once the calculation is done, drag and drop again the relaxed configuration to the Script generator. Add the following block, and modify the parameters as follows:

 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 kpoint Sampling to \(99 \times 99 \times 1\)
In the main panel of the Script Generator, set the Default output file to ‘Graphene_mulineiso.hdf5‘, send the script to the Editor, remove the Dynamical Matrix, Hamiltonian Derivatives and Electron Phonon Couplings blocks, and modifiy the Mobility block as follows:
# 
# Mobility
# 
mobility_full = nlread('Graphene_mulinefull.hdf5', Mobility)[1]
# Generate energydependent 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_mulineiso.hdf5', mobility)
Send the script to the Job manager, save it as ‘Graphene_mulineiso.py‘ and click on the button to run the calculation.
Once the calculation is done, open the file ‘Graphene_mulineiso.hdf5‘, select the 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 
 qgrid 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.91e05 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.25e03 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. energydependent method¶
A more stringent test for the reliability of the energydependent 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 Mobility analysis object.
 In the energydependent method, the twostep 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 energydependent 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\).
Theory section¶
The mobility \(\mu\) is related to the conductivity \(\sigma\) via
where \(n\) is the carrier density, \(e\) is the electronic charge.
In this tutorial, you will calculate the conductivity using the semiclassical Boltzmann transport equation (BTE). In the case of zero magnetic field and a timeindependent electric fields in the steady state limit the BTE simplifies to:
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 handside is approximated to linear order in the electric field by changing to the equilibrium distribution.
The electronphonon scattering is described by the collision integral. Most commonly this is approximated by a relaxation time:
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 electronphonon coupling matrix elements \(g_{\mathbf{k}\mathbf{k'}}^{\lambda n n'}\) are in ATK calculated using the ElectronPhononCoupling analysis module.
The general electronphonon collision integral is given by
Relaxation time approximation¶
In the relaxationtime approximation (RTA) one assumes a special simplified form of the RHS/collision integral:
where \(\delta f_{\mathbf{k}n} = f_{\mathbf{k}n}f^0_{\mathbf{k}n}\). The linearized BTE becomes:
Hereby we find the solution:
The relaxationtime 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:
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 relaxationtime:
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
where \(\mathbf{v}_{\mathbf{k}n}\) is the electron velocity. Following [GMSB16] we obtain
Mobility¶
Once knowing the relaxation times one obtain the mobility as:
where the factor 2 accounts for spin degeneracy.
References¶
[GMSB16]  (1, 2) T. Gunst, T. Markussen, K. Stokbro, and M. Brandbyge. Firstprinciples method for electronphonon coupling and electron mobility: Applications to twodimensional 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 electronphonon interaction in graphene. Physical Review B, 85(16):165440, 2012. 