Basic particle properties

Particle statistics

For electrons: The Fermi–Dirac distribution for a Bloch state (\(\boldsymbol{k}\), \(b\)) (\(\boldsymbol{k}\), \(b\) are Bloch indices) represents the occupation of electron states, and is:

\[f_{\boldsymbol{k},b} = \frac{1}{e^{(\frac{\epsilon_{\boldsymbol{k},b}-\mu}{k_BT})}+1}\]

For phonons: The Bose-Einstein distribution, which represents phonon occupation, is:

\[n_{\boldsymbol{k},b} = \frac{1}{e^{(\frac{\epsilon_{\boldsymbol{k},b}}{k_BT})}-1}\]

The derivatives of these functions show up in many transport equations. A computationally stable way to evaluate derivatives of these occupation functions is:

\[\frac{\partial n_{\boldsymbol{k},b}}{\partial T} = \frac{\epsilon_{\boldsymbol{k},b}}{4k_BT^2} \sinh( \frac{\epsilon_{\boldsymbol{k},b}}{2T} )\]

\[\frac{\partial f_{\boldsymbol{k},b}}{\partial T} = \frac{\epsilon_{\boldsymbol{k},b}-\mu}{4k_BT^2} \cosh( \frac{\epsilon_{\boldsymbol{k},b}-\mu}{2T} )\]

\[\frac{\partial n_{\boldsymbol{k},b}}{\partial \epsilon} = - \frac{1}{4k_BT} \sinh( \frac{\epsilon_{\boldsymbol{k},b}}{2T} )\]

\[\frac{\partial f_{\boldsymbol{k},b}}{\partial \epsilon} = - \frac{1}{4k_BT} \cosh( \frac{\epsilon_{\boldsymbol{k},b}-\mu}{2T} )\]

Specific heat

The specific heat at constant volume for phonons (or electrons, substituting in the appropriate occupation factors) is evaluated as:

\[C_v = \frac{1}{V N_k} \sum_{\boldsymbol{k},b} \epsilon_{\boldsymbol{k},b} \frac{\partial n_{\boldsymbol{k},b}}{\partial T}\]

where \(V\) is the volume of the primitive crystal unit cell, \(N_k\) is the number of wavevectors used, \(\boldsymbol{k}\) is a wavevector index, \(b\) is a branch/band index, \(\epsilon\) is the phonon energy (for electrons, this is replaced by \(\epsilon-\mu\)), \(T\) the temperature and \(n\) is the Bose–Einstein (Fermi–Dirac) distribution function.

Density of States

The density of states is defined as the number of states that are available to a particle at a certain energy, \(E\):

\[DOS(E) = \frac{1}{(2\pi)^d V} \sum_b \int_{BZ} \delta(E-\epsilon_{\boldsymbol{k}b}) d\boldsymbol{k}\]

where \(d\) is the dimensionality, \(V\) is the crystal unit cell volume, \(b\) is a Bloch index over bands, \(\boldsymbol{k}\) is a Bloch index over wavevectors, \(\epsilon\) is the particle energy and the integration is carried over the Brillouin zone.

The definition holds for both phonon and electrons.

The integral is sampled with a uniform grid of points in the Brillouin zone. In the DoS apps, the integral of the Dirac delta is implemented using the tetrahedron method.

Dirac delta approximations

We offer two possible methods for the approximation of the Dirac delta functions used in Phoebe’s transport calculations.

Gaussian Approximation

The delta function for the energy conservation can be replaced by a Gaussian:

\[\delta(\hbar \omega)=\frac{1} {\sqrt{\pi} \sigma} \exp{\left[-(\hbar \omega/ \sigma )^2 \right]} \;,\]

where the width, \(\sigma\) is a constant decided by user input. It is important to note that when the delta function is substituted with a Gaussian, the detailed balance condition is only valid under approximation. The definition used above guarantees that the scattering matrix is symmetric and non-negative.

Adaptive Gaussian

Another method is the adaptive Gaussian smearing scheme (see https://link.aps.org/doi/10.1103/PhysRevB.75.195121). We use this method to approximate a Dirac delta function of the form:

\[\delta(\hbar (\omega_1+\omega_2-\omega_3)=\frac{1} {\sqrt{\pi} \sigma} \exp{(-(\hbar (\omega_1+\omega_2-\omega_3)/ \sigma )^2)}.\]

In this method, \(\sigma\) is now dependent on the energies, and is not a user input. Specifically, we build it as:

\[\sigma = \frac{1}{\sqrt{12}} \sqrt{ \sum_{\beta} \left(\sum_{\alpha} (v_2-v_3) \frac{M_{\alpha \beta}}{N_{\beta}} \right)^2 }\]

where \(M\) is a matrix comprised of the primitive cell lattice vectors (each column is a lattice vector), \(v_2\) and \(v_3\) are the electron or phonon group velocities, and \(N_{\beta}\) is the number of wavevectors sampled along direction \(\beta\).

Note that the adaptive scheme may be critical in the case where the velocity sum to zero: in that case, we skip the scattering event, unless we have an exact energy conservation taking place.

Dynamical matrix

A density functional theory code can compute the following force constant matrix:

\[M(ls\alpha | l's'\alpha') = \frac{\partial^2 \mathcal{E}}{\partial u_{ls\alpha} \partial u_{l's'\alpha'}}\]

where \(M\) is a matrix of second order derivatives of the total crystal energy \(\mathcal{E}\) with respect to an ionic displacement \(u_{ls\alpha}\), where \(l\) labels a unit cell in a supercell, \(s\) is an index over the ionic basis, and \(\alpha\) denotes the direction in which the displacement is made. This matrix can either be computed with density functional perturbation theory (DFPT) or with a frozen-phonon approach. Due to the periodicity of the crystal, one has the freedom to set \(l=0\).

The dynamical matrix is the Fourier transform of this matrix. Excluding polar corrections, the dynamical matrix is:

\[D(s\alpha | s'\alpha')(\boldsymbol{q}) = \sum_{l'} M(0s\alpha | l's'\alpha') e^{i \boldsymbol{q} \cdot \boldsymbol{R}_{l'}}\]

Note that the Bravais lattice vectors are defined as the Bravais lattice vectors belonging to the Wigner-Seitz zone (not the Brillouin zone!) of a supercell, whose size is \(N_{qx}\times N_{qy}\times N_{qz}\) that of the primitive unit cell and this is the size of the q-point mesh used to compute the phonons in the DFT code.

The phonon energy and phonon eigenvectors are defined from by the eigenvalue problem,

\[D(s\alpha | s'\alpha')(\boldsymbol{q}) z_{s'\alpha'j}(\boldsymbol{q}) = \omega_{j}^2(\boldsymbol{q}) z_{s\alpha j}(\boldsymbol{q}).\]

where \(\omega_j(q)\) are the phonon energies, and \(z_{saj}(q)\) are the phonon eigenmodes.

Polar Correction

If ions carry a charge, one must not forget to add an additional term to D, representing a polar correction:

\[D(s\alpha | s'\alpha')(\boldsymbol{q}) += \frac{4\pi}{\Omega} e^2 \frac{ (\boldsymbol{q} \cdot Z^*_s)_{\alpha} (\boldsymbol{q} \cdot Z^*_{s'})_{\alpha'} } { (\boldsymbol{q} \cdot \epsilon_{\infty} \cdot \boldsymbol{q}) }\]

where \(Z_{s,\alpha,\beta}\) is the Born charge tensor of atom \(s\), and \(\epsilon_{\infty}\) is the static dielectric constant.