Physical and mathematical justification of the numerical Brillouin zone integration of the Boltzmann rate equation by Gaussian smearing

Scatterings of electrons at quasiparticles or photons are very important for many topics in solid state physics, e.g., spintronics, magnonics or photonics, and therefore a correct numerical treatment of these scatterings is very important. For a quantum-mechanical description of these scatterings Fermi's golden rule is used in order to calculate the transition rate from an initial state to a final state in a first-order time-dependent perturbation theory. One can calculate the total transition rate from all initial states to all final states with Boltzmann rate equations involving Brillouin zone integrations. The numerical treatment of these integrations on a finite grid is often done via a replacement of the Dirac delta distribution by a Gaussian. The Dirac delta distribution appears in Fermi's golden rule where it describes the energy conservation among the interacting particles. Since the Dirac delta distribution is a not a function it is not clear from a mathematical point of view that this procedure is justified. We show with physical and mathematical arguments that this numerical procedure is in general correct, and we comment on critical points.


Introduction
In solid-state physics, scatterings of electrons at periodic perturbations (quasiparticles or photons) are very important for many research fields and we give three examples in the following: 1. In all-optical switching experiments [1] a thin ferrimagnetic film, e.g., GdFeCo, is irradiated by a femtosecond laser pulse which can be linearly or circularly polarized and thereafter a demagnetization with subsequent switching of the magnetization can be observed under certain preconditions. The fundamental mechanisms are strongly debated at the moment, however, electron-photon scatterings, electron-phonon scatterings and electron-magnon scatterings certainly play a big role for the demagnetization of the ferrimagnetic film. 2. In ultrafast demagnetization experiments [2] a thin ferromagnetic film, e.g., Ni or Fe, is irradiated by a femtosecond laser pulse which is normally linearly polarized and thereafter an ultrafast demagnetization (on the timescale of about 100 fs) without switching of the magnetization can be observed. The magnetization recovers on a timescale of several picoseconds. Despite many years of research the fundamental

A r c h i v e o f S I D
first-order time-dependent perturbation theory. One can calculate the total transition rate from all initial states to all

A r c h i v e o f S I D
that this procedure is justified. We show with physical and mathematical arguments that this numerical procedure is in

A r c h i v e o f S I D
mathematical arguments that this numerical procedure is in general correct, and we comment on critical points.

A r c h i v e o f S I D Introduction
In solid-state physics, scatterings of electrons at periodic

A r c h i v e o f S I D
In solid-state physics, scatterings of electrons at periodic

A r c h i v e o f S I D
faehnle@is.mpg.de mechanisms are still unclear but scatterings of electrons at phonons [3,4] or at magnons [5] or at electrons [6] have been discussed intensively. 3. Spin-polarized currents are important for devices in spintronics [7], e.g., spin-transistors or spin-diodes. The lifetime of the spin-polarized electrons is crucial for the spintronics devices. The lifetimes are determined by scatterings of electrons at quasiparticles and at interfaces or defects.
A correct numerical calculation of the various scattering processes is important for the understanding of these effects in solid-state physics. In quantum mechanics, Fermi's golden rule gives the transition rate W k jk;j 0 k 0 from an initial electronic state W jk in a solid with energy e jk to a final electronic state W j 0 k 0 with energy e j 0 k 0 (j,j 0 : band indices; k, k 0 : wave vectors) due to a periodic perturbation arising from a (quasi)particle [8] AE hx qk is the energy of the involved (quasi)particle (q: wave vector, k: polarization) which may be, e.g., photons, phonons, magnons, plasmons etc., with frequency x qk for absorption (plus sign) or emission (minus sign), and M k jk;j 0 k 0 is the scattering matrix element where F j i and F 0 j i are the initial and final (quasi)particle states and W qk is the scattering operator. Thereby, momentum conservation k AE q ¼ k 0 þ G is demanded (G: reciprocal lattice vector). Fermi's golden rule is the firstorder approximation of the time-dependent quantum-mechanical perturbation theory. It implies that the scattering processes are Markovian which means that a scattering process does not depend on preceding scattering processes. Fermi's golden rule is only valid in a time window where the perturbation time on the one hand must be short enough because of the first-order approximation and on the other hand must be long enough to replace the sinðxÞ=x-function appearing in the derivation of Fermi's golden rule by the Dirac delta distribution. The validity of Fermi's golden rule for a magnetization dynamics on the 100 fs timescale is critically discussed in Ref. [4].
Normally, one is not interested in a specific transition rate W k jk;j 0 k 0 from an initial state W jk to a final state W j 0 k 0 but in the total transition rate W total from all initial states to all final states. Thereby k and k 0 are related via k AE q ¼ k 0 þ G if the scattering is at a quasiparticle with wave vector q. This is calculated with Boltzmann rate equations [4,9] where X BZ is the Brillouin zone (BZ) volume and n is the distribution function for the electrons. Often one is also interested in the rate of change of the distribution function n jk due to scattering which is also calculated with Boltzmann rate equations [10] Because the quantities e jk , e j 0 k 0 , W k j 0 k 0 ;jk , W k jk;j 0 k 0 can be calculated numerically only for a finite number of kpoints, finite k-point grids have to be used for the numerical calculation of the total transition rate W total or of the rate of change of the distribution function dn jk =dt. Thereby, energy conservation e j 0 k 0 ¼ e jk AE hx qk and momentum conservation k AE q ¼ k 0 þ G have to be fulfilled; however, energy conservation in combination with momentum conservation is in general never fulfilled for a finite k-point grid. Therefore, the Dirac delta distribution has to be replaced by a ''smeared'' delta function to obtain a result which approximates the integral (which is done, e.g., in Refs. [3,4,[11][12][13]

A r c h i v e o f S I D :
So we have to calculate Brillouin zone integrals of the form

A r c h i v e o f S I D
So we have to calculate Brillouin zone integrals of the form k g

A r c h i v e o f S I D
: Because the quantities

A r c h i v e o f S I D
Because the quantities calculated numerically only for a finite number of

A r c h i v e o f S I D
numerical calculation of the total transition rate of the rate of change of the distribution function d

A r c h i v e o f S I D
filled; however, energy conservation in combination with respect to the energy e but the integration is with respect to the wave vector k. The problem is explained in more detail in Sect. 2.
Mathematical proofs of Eq. (8) under certain preconditions can be found in Ref. [14], Theorem 7.2.1, and in Ref. [15], Theorem 6.1.5; however, the proofs are for general distributions and are very abstract. We want to show in this article that Eq. (8) is correct using also physical arguments.
In Sect. 2 we explain in detail the problem which arises when in a Brillouin zone integration Dirac's delta distribution is approximated by a Gaussian. This is done in many papers without giving any justification. We, therefore, think that the outline of this problem is a novelty per se. Then we give a justification of the Gaussian smearing method by mathematical and physical arguments. Each of these arguments has been used in other contexts in previous papers. The novelty of our paper is that we use these arguments to justify the Gaussian smearing method for the Brillouin zone integration. First, we consider a coordinate transformation from the wave vector variables k ¼ k x ; k y ; k z À Á to the variables ; #; u where is the energy and #; u are variables for the surface of constant energy. This transformation involves the Jacobian J k x ; k y ; k z À Á . The inverse function theorem [14] says that if this Jacobian is non-zero at a k-point, this transformation is invertible. Then the integration in k-space including d ð Þ can be represented by an integration over ; #; u of a function which involves the Jacobian e J ; #; u ð Þ ¼ J k x ; k y ; k z À Á Â Ã À1 and, which now can without any problem be replaced by a Gaussian. The problem is that there are special k-points where r k k ð Þ ¼ 0. For these special points the Jacobian J k x ; k y ; k z À Á is zero, and the reverse transformation involving e J ; #; u ð Þ is not defined in a rigorous mathematical interpretation. According to a general theorem of M. Morse the dispersion relation k ð Þ exhibits such special points because it is periodic in all components. There are special points which can be identified easily, e.g., the C-point and points on the Brillouin zone boundary. These points can be avoided by shifting the grid of k-points considered in the Brillouin zone integration accordingly [16]. Other special points cannot be easily found, and they might be in the shifted kpoint grid. Van Hove has shown [17] that for three dimensions the appearance of these special points does not appreciably modify the result of a numerical integration. We motivate these steps by physical reasoning.
In Sect. 3 we give practical hints for the appropriate choice of the smearing parameter r. Finally, our results are summarized in Sect. 4.

Numerical integration of the Dirac delta distribution
It is very well known that in integrals involving the Dirac delta distribution, the distribution can be replaced by a Gaussian for the limes r ! 0. It reads Z þ1 À1 de gðeÞ dðeÞ ¼ lim where gðeÞ is a continuously differentiable function which depends on the energy e. The Dirac delta distribution is approximated by a Gaussian Z þ1 for a numerical calculation of the integral and r has to be chosen appropriately, see Sect. 3. However, the integrals in Eqs. (3)-(6) are not over the energy e but over the wave vector k. For the sake of simplicity we discuss the following integral Z where gðkÞ is a continuously differentiable function of the wave vector k, and the generalization to the expression d e j 0 k 0 À ðe jk AE hx qk Þ À Á used in Eqs. (3)-(6) is straightforward. In explicit numerical calculations it is always assumed that also the relation holds without giving any justification, reference or comment and that this may be approximated by Z However, the Dirac delta distribution is defined by Eq. (9) and not by Eq. (12). We show in the following how the use of Eq. (13) can be justified. We consider a coordinate transformation from the wave vector variables k ¼ ðk x ; k y ; k z Þ to the variables e; #; u (e: energy; #; u: variables for the surface of constant energy)

A r c h i v e o f S I D
] says that if this Jaco-

A r c h i v e o f S I D 3. A r c h i v e o f S I D .
However, the integrals in Eqs. (

A r c h i v e o f S I D
where g

A r c h i v e o f S I D
ward. In explicit numerical calculations it is always assumed that also the relation

A r c h i v e o f S I D
assumed that also the relation www.SID.ir www.SID.ir 123 www.SID.ir 123 where the energy dispersion relation eðk x ; k y ; k z Þ is known and the surfaces of constant energy can be parametrized with two variables # and u. The inverse function theorem says [18] that every continuously differentiable, vectorvalued function which maps values from an open set of R n to other values of an open set of R n [so-called coordinate transformation, e.g., Eq. (14)] and whose Jacobian determinant If Eq. (14) is invertible in the neighborhood of every point k ¼ ðk x ; k y ; k z Þ-whereby only points k with eðkÞ ¼ 0 are relevant because of the Dirac delta distribution in Eq. (11)-it is possible to make for this neighborhood a local coordinate transformation (using Eq. (16)) for the function gðkÞ ¼ e gðe; #; uÞ appearing in Eq. (11). Then, the integral over the wave vector k can be replaced by the integral over the variables e; #; u Z BZ d 3 k gðkÞ dðeðkÞÞ Note that e J ðe; #; uÞ ¼ J À1 ðk x ; k y ; k z Þ with ðe; #; uÞ expressed by Eq. (14). It is now definitely allowed to approximate the Dirac delta distribution by a Gaussian in analogy to Eqs. (9) and (10) where in the last step the integration variables are changed back again to an integration over the wave vector using Eq. (14). This is exactly what we wanted to show in Eq. (13). One must keep in mind that the Jacobian determinant Jðk x ; k y ; k z Þ given by Eq. (15) is zero for special points k ¼ ðk x ; k y ; k z Þ, where r k eðkÞ ¼ 0. This is the case for the C-point and usually for points on the Brillouin zone boundary [19], and even the transformation (14) could be not continuously differentiable for special points. Then, the reverse transformation (16) used in Eq. (19) is not defined anymore in a rigorous mathematical interpretation. However, these problems arise because of two idealizations, the long-time idealization and the infinite-solid idealization, and the following remarks have to be considered:

A r c h i v e o f S I D
). This is exactly what we wanted to show in One must keep in mind that the Jacobian determinant

A r c h i v e o f S I D
One must keep in mind that the Jacobian determinant given by Eq. (

A r c h i v e o f S I D
) is zero for special points z

A r c h i v e o f S I D
, where r 19], and even the transformation (

A r c h i v e o f S I D
], and even the transformation ( not continuously differentiable for special points. Then, the

A r c h i v e o f S I D
)-it is possible to make for this neighborhood a local coordinate transformation (using Eq. (

A r c h i v e o f S I D
)) for the appearing in Eq. (

A r c h i v e o f S I D
). Then, the can be replaced by the

A r c h i v e o f S I D
can be replaced by the Z

A r c h i v e o f S I D
; #; u

A r c h i v e o f S I D
x k x k   16 anymore in a rigorous mathematical interpretation. How-

A r c h i v e o f S I D
anymore in a rigorous mathematical interpretation. However, these problems arise because of two idealizations, the

A r c h i v e o f S I D
ever, these problems arise because of two idealizations, the long-time idealization and the infinite-solid idealization,

A r c h i v e o f S I D
long-time idealization and the infinite-solid idealization, and the following remarks have to be considered:

A r c h i v e o f S I D
and the following remarks have to be considered: 1. In a physical interpretation the Dirac delta distribution appearing in Fermi's golden rule, Eq. (

A r c h i v e o f S I D
appearing in Fermi's golden rule, Eq. ( defined critical points (see also Ref. [16]). One could argue that a k-point very close to a critical k-point could yield an extremely large contribution because the Jacobian determinant e Jðe; #; uÞ (18) appearing in Eq. (17) may be very large for this point. This means that the second line of Eq. (19) would not be a good approximation for a choice of the k-point grid which contains points very close to critical points, but the near equality of the first line with the third line of Eq. (19) still holds because for the transition between the first and the third line the integration variables have been changed back again and this corresponds to the multiplication with e J À1 ðe; #; uÞ ¼ Jðk x ; k y ; k z Þ, so altogether, a possibly large value of e J ðe; #; uÞ does not matter. However, it is extremely complicated to identify all critical points and, therefore, it is not clear whether a chosen k-point grid contains critical points or not. To avoid these cumbersome investigations one can also do the following: one performs calculations for denser and denser grids and/or for shifted and rotated grids and compares the results. If the results are very similar, this means that the grids either do not contain critical points or that the critical points do not make a big contribution so that they do not falsify the results, in agreement with Ref. [17].
Practical hints for the appropriate choice of the smearing parameter The appropriate choice of the smearing parameter r appearing in Eq. 10 is crucial for the correct numerical calculation of the Boltzmann rate equation. The appropriate choice of r depends on two quantities: 1. First, the smearing parameter r depends on the energy scale of the involved (quasi)particle which may be, e.g., a photon, phonon, magnon, plasmon (see Sect. 1). The energy scale for a phonon is in the order of some mRy (about 40 meV) and for a magnon the energy scale is a factor 10 larger than for a phonon. The energy scale for a plasmon is much larger, in the order of 700 mRy (about 10 eV). An appropriate choice of the smearing parameter is in the same order as the energy scale of the involved (quasi)particle. 2. Second, the smearing parameter r depends on the grid spacing. A typical ansatz is r ¼ p=N 1 where N 1 is the number of k-points in one direction and the total number of k-points is N 3 1 . For a fixed proportionality constant p it is guaranteed that the smearing parameter is smaller, the larger the N 1 .
In the following we discuss the choice of r for the case of electron-phonon scatterings. In Ref. [3] the smearing parameter is fixed to 15 meV (about 1 mRy) for the numerical calculation of the electron-phonon Boltzmann rate equation. In Ref. [4] we tested our numerical results of the electron-phonon Boltzmann rate equation for many different grids and smearing parameters. In this publication we considered the case of ultrafast demagnetization, see Sect. 1. Among other quantities we calculated the rate of the magnetic moment change per atom dM=dtðt s Þ for a time t s (see Eq. (14) of Ref. [4]). t s is the time after the laser pulse irradiation where the electron system has thermalized, i.e., the electron distribution can be described by a Fermi-Dirac distribution with the electron temperature T e : Figure 1 shows the rate of the magnetic moment change per atom dM=dtðt s Þ for Fe and an electron temperature of T e ¼ 2000 K. We calculated dM=dtðt s Þ for different grids (number of k-points in one direction N 1 ) and for different smearing parameters r. One can see that the results for dM=dtðt s Þ depend hardly on the chosen grid and on the chosen smearing parameter except for N 1 ¼ 10. Therefore, the above-discussed critical points do not falsify our results and the smearing parameter is in the right order of magnitude (about several mRy). Of course it is trivial that increasing the number of k-points increases the convergence. By Fig. 1 we just want to show that the results depend only very slightly on the specific choice of r if the number of k-points is above a certain value.

Conclusions
Scattering processes of electrons at periodic perturbations are very important in solid-state physics and also a correct numerical treatment is crucial for the quantitative analysis of scattering processes in many research activities, e.g.,

Practical hints for the appropriate choice
The appropriate choice of the smearing parameter

A r c h i v e o f S I D
The appropriate choice of the smearing parameter r

A r c h i v e o f S I D
depends on the energy scale of the involved (quasi)particle which may be,

A r c h i v e o f S I D
scale of the involved (quasi)particle which may be, e.g., a photon, phonon, magnon, plasmon (see Sect.

A r c h i v e o f S I D
e.g., a photon, phonon, magnon, plasmon (see Sect. The energy scale for a phonon is in the order of some

A r c h i v e o f S I D
The energy scale for a phonon is in the order of some mRy (about 40 meV) and for a magnon the energy

A r c h i v e o f S I D
mRy (about 40 meV) and for a magnon the energy scale is a factor 10 larger than for a phonon. The

A r c h i v e o f S I D
scale is a factor 10 larger than for a phonon. The energy scale for a plasmon is much larger, in the order spintronics. In quantum mechanics, Fermi's golden rule is used which contains the Dirac delta distribution. The Dirac delta distribution is usually replaced by a Gaussian to integrate numerically the Boltzmann rate equations on a finite grid of k-points. It is not obvious from the very beginning that this numerical treatment is correct since the Dirac delta distribution is not a function and the smearing variable differs from the integration variable. We have shown in the present article that this procedure is in general correct. There are special k-points for which it is in principle not justified to replace the Dirac delta distribution by a Gaussian; however, we have given mathematical and physical arguments why this procedure is nevertheless a good approximation for the integration of the Boltzmann rate equation and should not falsify the results, at least for three dimensions. It is not clear whether the same holds also for d ¼ 2 or even for d ¼ 1. In conclusion, the naive replacement of the Dirac delta distribution by a Gaussian gives in general correct results for the Boltzmann rate equation but this has to be checked for denser and/or for shifted and rotated grids to avoid wrong contributions from the above-described special k-points where the replacement is critical.

A r c h i v e o f S I D
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://crea tivecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.