Theory¶

The theory described below and the methods in coolvib are described in detail in [Askerka2016] and [Maurer2016].

Non-adiabatic friction of adsorbate motion on metal surfaces¶

In insulating or most semi-conducting materials the energy spectrum of vibrations and electronic excitations are clearly separated. This is not the case for metals. In metals even the smallest vibrational motion can lead to electronic excitations due to the vanishing gap between occupied and unoccupied states.

Figure 1: Schematic view of adsorbate vibration (here shown for a CO on Cu(100) internal stretch motion) leading to changes in the electronic structure that excite electron hole-pairs from below to above the Fermi level of the metal Density-of-States.

As shown in Figure 1, adsorbate motion can lead to changes in the electronic structure that facilitate low energy electronic excitations in the metal band structure. We can treat this interaction between electrons and vibrations with perturbation theory if the following assumptions are true:

the coupling is weak compared to the individual contributions of electrons and nuclei

electron-hole pair excitations do not lead to a qualitative change in the nuclear dynamics

Non-adiabatic relaxation rates from first order time-dependent perturbation theory¶

The relaxation rate of a vibrational mode $\omega_{\mathbf{q}j}$ due to nonadiabatic coupling is given by Fermi’s Golden Rule as derived from time-dependent perturbation theory and is given by:

(1) $\Gamma_{\mathbf{q}j} = \frac{\pi\hbar}{M}\sum_{s\mathbf{k}}\sum_{\nu,\nu'} (f(\epsilon_{s\mathbf{k+q}\nu'})-f(\epsilon_{s\mathbf{k}\nu}))|g_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{\mathbf{q}j}|^2 \cdot (\epsilon_{s\mathbf{k+q}\nu'}-\epsilon_{s\mathbf{k}\nu}) \cdot\delta(\epsilon_{s\mathbf{k+q}\nu'}-\epsilon_{s\mathbf{k}\nu}-\hbar\omega_{\mathbf{q}j}),$

where the sums run over spin, crystal momentum vectors $\mathbf{k}$ and all pairs of eigenvectors $\nu$ and $\nu'$ with the state occupations $f$ given by the Fermi-Dirac distribution and the coupling elements $g^{\mathbf{q}j}_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}$

(2) $g_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{\mathbf{q}j} = \braket{\psi_{s\mathbf{k+q}\nu'}|\mathbf{e}_{\mathbf{q}j}\nabla_{\mathbf{R}}\psi_{s\mathbf{k}\nu}}.$

For a detailed derivation of these equations see

Derivation

The relaxation rate relates to the lifetime of the mode by $\tau=1/\Gamma$ and to the vibrational linewidth by $\gamma=\Gamma\hbar$ . In eq. (2) $\mathbf{e}_{\mathbf{q}j}$ is the atomic displacement vector corresponding to a given vibrational motion $\omega_{\mathbf{q}j}$ described by a wave-vector $\mathbf{q}$ and a quantum number $j$ . $\nabla_{\mathbf{R}}$ represents the derivative vector with respect to all cartesian atomic positions.

It is often assumed in literature [Head-Gordon1995], [Hellsing1984] that perturbations $\hbar\omega$ due to vibrations are small and the density of state is almost constant in a large energy range around the Fermi level. This so-called quasi-static approximation leads to following simplification:

(3) $\Gamma_{\mathbf{q}j} = \frac{\pi\hbar}{M}\sum_{s\mathbf{k}}\sum_{\nu,\nu'} (f(\epsilon_{s\mathbf{k+q}\nu'})-f(\epsilon_{s\mathbf{k}\nu}))|g_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{\mathbf{q}j}|^2 \cdot (\epsilon_{s\mathbf{k+q}\nu'}-\epsilon_{s\mathbf{k}\nu}) \cdot\delta(\epsilon_{s\mathbf{k+q}\nu'}-\epsilon_{s\mathbf{k}\nu}).$

The Friction Tensor¶

We can reformulate the matrix elements of expression (2) in terms of a cartesian coupling tensor:

$|g_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{\mathbf{q}j}|^2 &= g_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{\mathbf{q}j}\cdot g_{s\mathbf{k}\nu,s\mathbf{k+q}\nu'}^{\mathbf{q}j} = \\ & \mathbf{e}^T_{\mathbf{q}j}\cdot\braket{\psi_{s\mathbf{k}\nu}|\mathbf{\nabla}_{\mathbf{R}}|\psi_{s\mathbf{k+q}\nu'}}\cdot\braket{\psi_{s\mathbf{k+q}\nu'}|\mathbf{\nabla}_{\mathbf{R}}|\psi_{s\mathbf{k}\nu}} \cdot\mathbf{e}_{\mathbf{q}j} = \\ & \mathbf{e}^T_{\mathbf{q}j}\cdot \mathbf{G}_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}\cdot \mathbf{e}_{\mathbf{q}j}.$

The elements of G are defined as:

$\mathbf{G}_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{n'a',na} = \braket{\psi_{s\mathbf{k}\nu}|\frac{\partial}{\partial_{a}\mathbf{R}_{n}}|\psi_{s\mathbf{k+q}\nu'}} \cdot \braket{\psi_{s\mathbf{k+q}\nu'}|\frac{\partial}{\partial_{a'}\mathbf{R}_{n'}}|\psi_{s\mathbf{k}\nu}} ,$

where n and n’ indicate the n-th (n’-th) atom and a and a’ indicate one of the three cartesian directions x, y, and z. Inserting G into eq. (3) we arrive at:

(4) $\Gamma(\omega_{\mathbf{q}j}) &= \mathbf{e}^T_{\mathbf{q}j}\cdot \biggl( \pi\sum_{s,\mathbf{k},\nu,\nu'} \mathbf{G}_{s\mathbf{k+q}\nu',s\mathbf{k}\nu} \cdot(f(\epsilon_{\mathbf{k}\nu})-f(\epsilon_{\mathbf{k+q}\nu'})) \\ \nonumber & \cdot\delta(\epsilon_{\mathbf{k}\nu}-\epsilon_{\mathbf{k+q}\nu'}) \biggr) \cdot \mathbf{e}_{\mathbf{q}j} = \mathbf{e}^T_{\mathbf{q}j}\cdot\Lambda^{\mathbf{q}}\cdot\mathbf{e}_{\mathbf{q}j}.$

In eq. (4) we have defined the friction tensor $\mathbf{\Lambda}^{\mathbf{q}}$ with dimensions $(3N\times3N)$ where N is the number of atoms. The elements of $\mathbf{\Lambda}^{\mathbf{q}}$ are defined as

(5) $\Lambda^{\mathbf{q}}_{n'a',na} = \frac{\pi\hbar}{\sqrt{M_{n'}}\sqrt{M_n}} \sum_{s,\mathbf{k},\nu,\nu'} G_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{n'a',na} \cdot(f(\epsilon_{\mathbf{k}\nu})-f(\epsilon_{\mathbf{k+q}\nu'}))\cdot(\epsilon_{s\mathbf{k+q}\nu'}-\epsilon_{s\mathbf{k}\nu}) \cdot\delta(\epsilon_{\mathbf{k}\nu}-\epsilon_{\mathbf{k+q}\nu'}) .$

Vibrational relaxation rates of adsorbates¶

Eq. (4) shows that using the friction tensor one can calculate relaxation rates and vibrational lifetimes for arbitrary molecular motions.

Figure 2: Collective motion of a periodic overlayer of CO atoms on Cu(100)

Considering periodic boundary conditions, the collective motion of an adsorbate overlayer (as shown in Fig. 2) can be described with a single adsorbate vibrational mode with wave vector $q=0$ . The corresponding relaxatio rate for this vibration is given by:

(6) $\Gamma(\omega_j) = \mathbf{e}^T_j\cdot\Gamma^{\mathbf{q}=0}\cdot\mathbf{e}_j$

Calculation of vibrational cooling using this approach is the most common [Head-Gordon1992] , [Persson1982] and is the default setting when calculating lifetimes in coolvib.

Figure 3: Motion of only one CO atom adsorbed on Cu(100), effectively breaking periodicity.

When considering the lifetime due to motion of a single adsorbate, such as is the case for impinging adsorbates, low coverages or non-collective adlayer motion (see Fig. 3) one effectively needs to integrate over all possible vibrations of the same type with different wave vectors $\mathbf{q}$ , effectively amounting to a Fourier expansion in reciprocal space. Correspondingly, the relaxation rate is given by:

(7) $\Gamma(\omega_{\mathbf{q}j}) = \sum_{\mathbf{q}} \mathbf{e}^T_{\mathbf{q}j}\cdot\Gamma^{\mathbf{q}}\cdot\mathbf{e}_{\mathbf{qj}} \approx \mathbf{e}^T_j\cdot\Gamma^{\sum \mathbf{q}}\cdot\mathbf{e}_j$

In the corresponding relaxation rate, we additionally include excitations that in principle violate momentum-conservation. If we assume that adsorbate vibrations interact only by simple phase modulation we can assume the atomic displacements of individual adsorbate images to be independent of wave vector yielding following expression for the friction tensor in this case:

(8) $\Lambda^{\sum \mathbf{q}}_{n'a',na} = \frac{\pi\hbar}{\sqrt{M_{n'}}\sqrt{M_n}} \sum_{s} \sum_{\mathbf{q}}\sum_{\mathbf{k}} \sum_{\nu,\nu'} G_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{n'a',na} \cdot(f(\epsilon_{\mathbf{k}\nu})-f(\epsilon_{\mathbf{k+q}\nu'}))\cdot(\epsilon_{s\mathbf{k+q}\nu'}-\epsilon_{s\mathbf{k}\nu}) \cdot\delta(\epsilon_{\mathbf{k}\nu}-\epsilon_{\mathbf{k+q}\nu'}) .$

Numerical considerations and working equations¶

Depending on if one wants to calculate vibrational lifetimes of periodic overlayers or single adsorbates one is left with evaluating eqs. (5) or (8). This is by no means a simple task due to the slow convergence of electronic structure close to the Fermi level and intricacies of evaluating a delta function in a non-equidistant space of excitations.

Some of these issues can be solved by transforming the equation into a continuous representation using the the Density-of-State of the system:

$\sum_{\nu} 1= \sum_{\nu} \int_{-\infty}^{\infty} d\epsilon_a\cdot \delta(\epsilon_a-\epsilon_{\nu})= \int_{-\infty}^{\infty} d\epsilon_a\cdot \underbrace{\sum_{\nu} \delta(\epsilon_a-\epsilon_{\nu})}_{\rho(\epsilon)}$

Introducing the Density-of-State for both sums over eigenstates in (5) we arrive at

$\Lambda^{\mathbf{q}}_{n'a',na} = \frac{\pi\hbar}{\sqrt{M_{n'}}\sqrt{M_n}} \int_{-\infty}^{\epsilon_F} d\epsilon_a\int_{\epsilon_F}^{\infty} d\epsilon_b \sum_{s,\mathbf{k}}\sum_{\nu}\sum_{\nu'} G_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{n'a',na} \cdot(f(\epsilon_{\mathbf{k}\nu})-f(\epsilon_{\mathbf{k+q}\nu'}))\cdot\delta(\epsilon_a-\epsilon_{s\mathbf{k}\nu})\cdot\delta(\epsilon_b-\epsilon_{s\mathbf{k+q}\nu'})\cdot\epsilon \cdot\delta(\epsilon_{a}-\epsilon_{b}) .$

The change in variable from the discrete space to the continuous space for the last delta function and the energy difference can be understood as a consequence of sequentially evaluating the first to delta functions (and by equivalence) Now we change the integration variables from occupied and unoccupied states to excitation energies

$\epsilon = \epsilon_b - \epsilon_a \Rightarrow d\epsilon_b=d\epsilon+\epsilon_a$

and arrive at

$\Lambda^{\mathbf{q}}_{n'a',na} = \frac{\pi\hbar}{\sqrt{M_{n'}}\sqrt{M_n}} \int_{-\infty}^{\epsilon_F} d\epsilon_a\int_{\epsilon_F}^{\infty} (d\epsilon+d\epsilon_a) \sum_{s,\mathbf{k}}\sum_{\nu}\sum_{\nu'} G_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{n'a',na} \cdot(f(\epsilon_{\mathbf{k}\nu})-f(\epsilon_{\mathbf{k+q}\nu'}))\cdot\delta(\epsilon_a-\epsilon_{s\mathbf{k}\nu})\cdot\delta(\epsilon+\epsilon_a-\epsilon_{s\mathbf{k+q}\nu'})\cdot\epsilon \cdot\delta(-\epsilon) =$

(9) $&= \frac{\pi\hbar}{\sqrt{M_{n'}}\sqrt{M_n}} \int_{-\infty}^{\epsilon_F} d\epsilon_a\int_{\epsilon_F}^{\infty} d\epsilon_a \sum_{s,\mathbf{k}}\sum_{\nu}\sum_{\nu'} G_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{n'a',na} \cdot(f(\epsilon_{\mathbf{k}\nu})-f(\epsilon_{\mathbf{k+q}\nu'}))\cdot\delta(\epsilon_a-\epsilon_{s\mathbf{k}\nu})\cdot\delta(\epsilon+\epsilon_a-\epsilon_{s\mathbf{k+q}\nu'})\cdot\epsilon \cdot\delta(-\epsilon) = \\ &= \frac{\pi\hbar}{\sqrt{M_{n'}}\sqrt{M_n}} \int_{-\infty}^{\epsilon_F} d\epsilon_a\int_{-\infty}^{\infty} d\epsilon \sum_{s,\mathbf{k}}\sum_{\nu}\sum_{\nu'} G_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{n'a',na} \cdot(f(\epsilon_{\mathbf{k}\nu})-f(\epsilon_{\mathbf{k+q}\nu'}))\cdot\delta(\epsilon_a-\epsilon_{s\mathbf{k}\nu})\cdot\delta(\epsilon+\epsilon_a-\epsilon_{s\mathbf{k+q}\nu'})\cdot\epsilon \cdot\delta(-\epsilon)$

The first term in eq. (9) vanishes because the integration variable for occupied states $d\epsilon_a$ is not defined in the range between $\infty$ and $\epsilon_F$ . In the second term we can use the following properties of delta functions

$\int \delta(\epsilon_a-\epsilon_{s\mathbf{k}\nu})\delta(\epsilon+\epsilon_a-\epsilon_{s\mathbf{k+q}\nu'})=\delta(\epsilon+\epsilon_{s\mathbf{k}\nu}-\epsilon_{s\mathbf{k+q}\nu'})$

and

$\delta(-x) = \delta(x)$

to arrive at

(10) $\Lambda^{\mathbf{q}}_{n'a',na} = \int_{-\infty}^{\infty} d\epsilon \cdot\epsilon \cdot\delta(\epsilon) \cdot\underbrace{\sum_{s,\mathbf{k}}\sum_{\nu}\sum_{\nu'}\frac{\pi\hbar}{\sqrt{M_{n'}}\sqrt{M_n}} G_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{n'a',na} \cdot(f(\epsilon_{\mathbf{k}\nu})-f(\epsilon_{\mathbf{k+q}\nu'}))\cdot\delta(\epsilon_{s\mathbf{k+q}\nu'}-\epsilon_{s\mathbf{k}\nu}-\epsilon) }_{A_{n'a',na}(\epsilon)}$

By rewriting the friction tensor in this way we have gained several things:

The continuous formulation enables simple numerical discretization and accurate evaluation of the delta function $\delta(\epsilon)$ , thereby reducing numerical instabilities.

We have defined a set of electron-phonon spectral functions $A_{n'a',na}(\epsilon)$ for all combinations of cartesian components, which describe the principal spectrum of available coupling channels between electronic and nuclear degrees of freedom.

The second delta function transforms the set of finite transitions into smooth continuous spectral functions and therefore improves convergence with respect to Brillouin zone sampling.

The expression for the spectral functions

(11) $A_{n'a',na}(\epsilon) = \sum_{s,\mathbf{k}}\sum_{\nu}\sum_{\nu'}\frac{\pi\hbar}{\sqrt{M_{n'}}\sqrt{M_n}} G_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{n'a',na} \cdot(f(\epsilon_{\mathbf{k}\nu})-f(\epsilon_{\mathbf{k+q}\nu'}))\cdot\delta(\epsilon_{s\mathbf{k+q}\nu'}-\epsilon_{s\mathbf{k}\nu}-\epsilon)$

and its corresponding equivalent for isolated adsorbates

(12) $A_{n'a',na}^{\sum \mathbf{q}}(\epsilon) = \sum_{s,\mathbf{q},\mathbf{k}}\sum_{\nu}\sum_{\nu'}\frac{\pi\hbar}{\sqrt{M_{n'}}\sqrt{M_n}} G_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{n'a',na} \cdot(f(\epsilon_{\mathbf{k}\nu})-f(\epsilon_{\mathbf{k+q}\nu'}))\cdot\delta(\epsilon_{s\mathbf{k+q}\nu'}-\epsilon_{s\mathbf{k}\nu}-\epsilon)$