Derivation¶

The following derivations are directly taken from [Maurer2015].

The electron-phonon interaction¶

Our starting point is the molecular many-body problem in the Born-Oppenheimer approximation. Herein we define an electronic Schrödinger equation (SGE):

$\hat{H}_{\mathrm{e}}\Phi(t,\mathbf{r};\mathbf{R}) = i\hbar\frac{d}{dt} \Phi(t,\mathbf{r};\mathbf{R})$

with

(1) $\hat{H}_{\mathrm{e}} &= -\sum_{i=1}^N\frac{\hbar^2}{2}\nabla_i^2 +\sum_{i=1}^N\sum_{j>i}^N\frac{1}{r_{ij}}-\sum_{i=1}^N\sum_{A=1}^M \frac{Z_A}{r_{iA}}+ \\ \nonumber &+ \sum_{A=1}^M\sum_{B>A}^M \frac{Z_AZ_B}{R_{AB}}$

and a nuclear Schrödinger equation as

$\hat{H}_{\mathrm{n}}\chi(t,\mathbf{R}) = i\hbar\frac{d}{dt} \chi(t,\mathbf{R})$

with

$\hat{H}_{\mathrm{n}} = -\sum_{A=1}^M\frac{\hbar^2}{2M_A}\nabla_A^2 + E_{\mathrm{e}}(\mathbf{R}) .$

All of the above is given in atomic units, with the exception of explicitly stating $\hbar(=1)$ . The electronic SGE and its eigenfunctions depend parametrically on the nuclear coordinates. Due to separation of time-scales at which electrons and nuclei move, these equations are solved independently.

We can introduce coupling between the two systems by constructing the full molecular wavefunction as linear combination of products of nuclear and adiabatic electronic wavefunctions in a so-called Born expansion

$\Psi(t,\mathbf{r},\mathbf{R}) = \sum_j \chi_j(t,\mathbf{R}) \Phi_j(\mathbf{r};\mathbf{R})$

The coupling is introduced by the third term in eq. (1). Insertion into the full time-dependent Schrödinger equation, multiplication from the left with an electronic wavefunction $\Phi_i$ and integration over electronic degrees of freedom we find an additional term in the nuclear Schrödinger equation:

$(\underbrace{\hat{T}_{{\mathrm{n}}} + E^{\mathrm{e}}_j}_{\hat{H}_{\mathrm{n}}}+ \hat{H}_{\mathrm{ep}})\chi_j = i\hbar\frac{d}{dt}\chi_j$

with the action of the electron-phonon coupling term on the nuclear wavefunction defined as

(2) $\hat{H}_{\mathrm{ep}}\chi_j =-\frac{\hbar^2}{2M}\sum_{i\neq j}\left[2\langle\Phi_i|\nabla_{\mathbf{R}}\Phi_j\rangle\cdot\nabla_{\mathbf{R}}+\langle\Phi_i|\nabla^2_{\mathbf{R}}\Phi_j\rangle \right]\chi_i$

One way to understand this interaction is by viewing it as hybridization between electronic and nuclear wavefunctions. Each nuclear wavefunction is split into as many states as there are electronic excitations (and vice versa). Throughout the remainder of this documentation we will assume $\hat{H}_{\mathrm{ep}}$ to be small compared to $\hat{H}_n$ , i.e. it does not change the nuclear dynamics qualitatively. This enables application of perturbation theory. In the following we will find an expression for the interaction matrix elements and subsequently use Time Dependent Perturbation Theory to derive our working expressions.

Derivation of coupling matrix elements¶

Next we will find an expression for the matrix elements

$g_{fj} = \braket{\chi_f|\hat{H}_{\mathrm{ep}}|\chi_j} = -\frac{\hbar^2}{2M}\sum_{i}\bra{\chi_f}\biggl[2\braket{\Phi_i|\nabla_{\mathbf{R}}\Phi_j}\cdot\nabla_{\mathbf{R}} +\braket{\Phi_i|\nabla^2_{\mathbf{R}}\Phi_j} \biggr]\ket{\chi_i}.$

Our starting point is a periodic condensed phase material in which the electronic structure underlying the nuclear motion is given by effectively independent quasi-particles such as calculated with semi-local approximations to Density Functional Theory or the GW method. The unperturbed nuclear motion is described in terms of harmonic vibrations as given by the following Hamiltonian:

(3) $H_0 = \sum_{\mathbf{q}j} \hbar \omega_{\mathbf{q}j}\left(b^\dagger_{\mathbf{q}j}b_{\mathbf{qj}} +\frac{1}{2} \right) .$

Herein $b^\dagger_{\mathbf{q}j}$ and $b_{\mathbf{qj}}$ are the usual particle annihilation and creation operators for a phonon $j$ at wavevector $\mathbf{q}$ . Phonon normal mode displacement vectors are given as $\mathbf{e}_{\mathbf{q}j}$ . With this choice the scalar coupling term in eq. (2) vanishes due to orthonormality of the nuclear wavefunctions. The matrix element between two many-body nuclear wavefunction states, which only differ by a single one-particle state reduces to

$\braket{\chi_f|\hat{H}_{\mathrm{ep}}|\chi_j} \Rightarrow \braket{\chi_{\mathbf{q}f}|\hat{H}_{\mathrm{ep}}|\chi_{\mathbf{q}j}}.$

Furthermore, by taking the many-body electronic wavefunction as a product of quasi-particle states with quantum states specified by the quantum numbers $s$ , $\mathbf{k}$ , and $\nu$ representing spin, wave vector, and band index, all possible single-particle excited states are defined with the above indices for electrons and phonons. Correspondingly the matrix elements are

$\braket{\chi_{\mathbf{q}f}|\hat{H}_{\mathrm{ep}}|\chi_{\mathbf{q}j}} = -\frac{\hbar^2}{M} \sum_{i}\cdot \braket{\chi_{\mathbf{q}f}|\sum_{s}\sum_{\mathbf{k}\nu\mathbf{k}'\nu'} \braket{\psi_{s\mathbf{k}\nu}|\mathbf{e}_{\mathbf{q}j}\nabla_{\mathbf{R}}\psi_{s\mathbf{k}'\nu'}}\mathbf{e}_{\mathbf{q}j}\nabla_{\mathbf{R}}|\chi_{\mathbf{q}i}} ,$

where the different sums run over occupied and unoccupied states in the ground state, respectively. The occupations of the one-particle electronic states and thereby the Fermi level are defined by Fermi-Dirac statistics. With harmonic displacements as starting point we have defined the nuclear derivative in $H_{ep}$ in the basis of these displacements.

The matrix element describing the coupling between two nuclear wavefunctions which differ by replacing a single phonon $\mathbf{q}j$ with a phonon $\mathbf{q}j'$ therefore reads

$g^{\mathbf{q}f,\mathbf{q}j} = -\frac{\hbar^2}{M}\sum_i \braket{\chi_{\mathbf{q}f}|\mathbf{e}_{\mathbf{q}j} \nabla_{\mathbf{R}}\chi_{\mathbf{q}i}} \cdot \sum_{s\mathbf{k}}w_{\mathbf{k}}\sum_{\nu<\epsilon_F}\sum_{\nu'>\epsilon_F} \braket{\psi_{s\mathbf{k+q}\nu'}|\mathbf{e}_{\mathbf{q}j}\nabla_{\mathbf{R}}\psi_{s\mathbf{k}\nu}} ,$

where Brillouin sampling weights $w_{\mathbf{k}}$ have been accounted for and final state indices follow from momentum conservation. Utilizing the ladder operator technique for harmonic oscillator eigenstates [Mahan2013], [Head-Gordon1992] we find

(4) $\braket{\chi_{\mathbf{q}f}| \mathbf{e}_{\mathbf{q}i}\nabla_{\mathbf{R}}\chi_{\mathbf{q}i}} &= \sqrt{\frac{\omega_{\mathbf{q}i} M}{2\hbar}} (\sqrt{n_{\mathbf{q},f}}\hat{\delta}_{\mathbf{q}f,\mathbf{q},i+1}- \sqrt{n_{\mathbf{q},f}+1}\hat{\delta}_{\mathbf{q}f,\mathbf{q},i-1}) .$

As a result transitions can only occur between vibrational states with $f=i\pm1$ . If we assume that we start out in the first vibrationally excited state and are interested in decay to the ground state we arrive at following expression for the matrixelement

(5) $g^{\mathbf{q}j} &= \sqrt{\frac{\hbar^3\omega_{\mathbf{q}j}}{2M}}\sum_{s\mathbf{k}}w_{\mathbf{k}} \sum_{\nu<\epsilon_F}\sum_{\nu'>\epsilon_F} g_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{\mathbf{q}j} .$

with

(6) $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}}.$

Therefore, the excitation or deexcitation probability of a phonon $\mathbf{q}j$ is given by the coupling of all possible electronic excitations that do not violate momentum conservation.

Often it is more desirable to express the matrix elements $g_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{\mathbf{q}j}$ in a different form, namely in terms of nuclear derivatives of the unperturbed Hamiltonian and overlap matrices. In the case of a local atomic orbital representation of the wavefunctions, such a representation was derived by Head-Gordon and Tully [Head-Gordon1992] as well as Savrasov and Savrasov [Savrasov1996] in the current context. To review this derivation we express the wavefunctions in terms of a non-orthogonal, incomplete basis set

(7) $\ket{\psi_{s\mathbf{k}i}} = \sum_a c_{s\mathbf{k}i}^a \ket{\phi_{\mathbf{k}}^a}$

such that following generalized eigenvalue matrix equation is satisfied (In the following intermediate results we will drop the spin index):

(8) $\mathbf{H}_{\mathbf{k}} \mathbf{c}_{\mathbf{k}i} = \mathbf{S}_{\mathbf{k}} \mathbf{c}_{\mathbf{k}i}\epsilon_{\mathbf{k}i}$

with

$H_{\mathbf{k}}^{ab} = \braket{\phi_{\mathbf{k}}^b|\hat{H}_e|\phi_{\mathbf{k}}^a}$

and

$S_{\mathbf{k}}^{ab} = \braket{\phi_{\mathbf{k}}^b|\phi_{\mathbf{k}}^a} .$

Substituting (7) into the Dirac-braket of eq. (5) gives

(9) $\braket{\psi_{s\mathbf{k+q}\nu'}|\mathbf{e}_{\mathbf{q}j}\nabla_{\mathbf{R}}\psi_{s\mathbf{k}\nu}} = \sum_{a,b} c^b_{s\mathbf{k+q}\nu'}\braket{\phi_{\mathbf{k+q}}^b|\mathbf{e}_{\mathbf{q}j}\nabla_{\mathbf{R}}\phi_{\mathbf{k}}^a}c^a_{s\mathbf{k}\nu} +c^b_{s\mathbf{k+q}\nu'}\braket{\phi_{\mathbf{k+q}}^b|\phi_{\mathbf{k}}^a}\mathbf{e}_{\mathbf{q}j}\nabla_{\mathbf{R}}c^a_{s\mathbf{k}\nu}$

The second term of this equality can be reexpressed as follows using the nuclear derivative of eq. (8)

(10) $c^b_{s\mathbf{k+q}\nu'}S_{\mathbf{k+q}\mathbf{k}}^{ab}\mathbf{e}_{\mathbf{q}j}\nabla_{\mathbf{R}} c^a_{s\mathbf{k}\nu} = \frac{1}{\epsilon_{\mathbf{k+q}\nu'}-\epsilon_{\mathbf{k}\nu}} c^b_{s\mathbf{k+q}\nu'}\left(H_{\mathbf{k+q},\mathbf{k}}^{\mathbf{q}j,ab}-\epsilon_{\mathbf{k}\nu} S_{\mathbf{k+q},\mathbf{k}}^{\mathbf{q}j,ab} \right) c^a_{s\mathbf{k}\nu}$

with

$S_{\mathbf{k+q},\mathbf{k}}^{ab} = \braket{\phi_{\mathbf{k+q}}^b|\phi_{\mathbf{k}}^a}$

and the nuclear derivative of an operator in matrix representation defined as

$O_{\mathbf{k+q},\mathbf{k}}^{\mathbf{q}j,ab} = \mathbf{e}_{\mathbf{q}j}\nabla_{\mathbf{R}}\braket{\phi_{\mathbf{k+q}}^b|O|\phi_{\mathbf{k}}^a} .$

Inserting this result into eq. (5) and changing to matrix representation we find

(11) $g_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{\mathbf{q}j} \approx\frac{ \mathbf{c}_{s\mathbf{k+q}\nu'}\left(\mathbf{H}_{\mathbf{k+q},\mathbf{k}}^{\mathbf{q}j}-\epsilon_{F} \mathbf{S}_{\mathbf{k+q},\mathbf{k}}^{\mathbf{q}j} \right) \mathbf{c}_{s\mathbf{k}\nu}} {\epsilon_{s\mathbf{k+q}\nu'}-\epsilon_{s\mathbf{k}\nu}}$

In eq. (11) we assume that the energy difference between states $\epsilon_{s\mathbf{k}\nu}$ and $\epsilon_{s\mathbf{k+q}\nu'}$ is very small and that both are close to the Fermi level and can be replaced therewith [Head-Gordon1992].

Unfortunately the above terms can not straightforwardly be extracted from an electronic structure calculation in reciprocal space representation due to the dependence on different wavevectors $\mathbf{k}$ and $\mathbf{k+q}$ and $\braket{\mathbf{k}|\mathbf{H}|\mathbf{k}'}\neq \delta_{\mathbf{k}\mathbf{k}'}\cdot\braket{\mathbf{k}|\mathbf{H}|\mathbf{k}}$ . In the case of a plane-wave representation a solution to this problem has been given by Trail et al. with coupling elements given in terms of eq. (?) [Trail2001]. In the current case, we can approach this by transforming the Hamiltonian and overlap matrices themselves to their real-space representations using a Fourier sum over wavevectors within the first Brillouin zone:

(12) $\mathbf{H}(|\mathbf{N}-\mathbf{N}'|)=\sum_{\mathbf{k}} w_{\mathbf{k}} e^{i\mathbf{k}\cdot(\mathbf{R}_{\mathbf{N}'}-\mathbf{R}_{\mathbf{N}})}\mathbf{H}_{\mathbf{k}}$

and

(13) $\mathbf{S}(|\mathbf{N}-\mathbf{N}'|)= \sum_{\mathbf{k}} w_{\mathbf{k}} e^{i\mathbf{k}\cdot(\mathbf{R}_{\mathbf{N}'}-\mathbf{R}_{\mathbf{N}})} \mathbf{S}_{\mathbf{k}} .$

In eq. (12) and (13), $\mathbf{N}$ and $\mathbf{N}'$ are integer index-vectors that define the absolute position of a unit cell image in the crystal volume $\Omega$ and $w_{\mathbf{k}}$ is the normalization weight of wavevector $\mathbf{k}$ . We have furthermore defined crystal translation vectors as

$\mathbf{R}_\mathbf{N} = \mathbf{R}_0 + \mathbf{T}_\mathbf{N},$

where $\mathbf{R}_0$ is the center of the principal unit cell and $\mathbf{T}_\mathbf{N}$ is a lattice translation vector to the unit cell image with index vector $\mathbf{N}$ [Soler2002], [Blum2009].

From the above results $\mathbf{H}^{\mathbf{q}j}(|\mathbf{N}-\mathbf{N}'|)$ and $\mathbf{S}^{\mathbf{q}j}(|\mathbf{N}-\mathbf{N}'|)$ can be directly calculated using the finite-difference method or Density Functional Perturbation Theory (Coupled Perturbed Kohn-Sham) [Savrasov1996], [Baroni2001]. With the definition of the inverse transformation

$\mathbf{H}^{\mathbf{q}j}_{\mathbf{k},\mathbf{k}'}=\sum_{\mathbf{N},\mathbf{N}'\in\Omega} e^{i(\mathbf{k}'\cdot\mathbf{R}_{\mathbf{N}'}-\mathbf{k}\cdot\mathbf{R}_{\mathbf{N}})} \mathbf{H}^{\mathbf{q}j}(|\mathbf{N}-\mathbf{N}'|)$

and

$\mathbf{S}^{\mathbf{q}j}_{\mathbf{k},\mathbf{k}'}= \sum_{\mathbf{N},\mathbf{N}'\in\Omega} e^{i(\mathbf{k}'\cdot\mathbf{R}_{\mathbf{N}'}-\mathbf{k}\cdot\mathbf{R}_{\mathbf{N}})} \mathbf{S}^{\mathbf{q}j}(|\mathbf{N}-\mathbf{N}'|)$

we may finally express the coupling element of eq. (11) as

$& g_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{\mathbf{q}j} \approx\sqrt{\frac{\hbar}{M}}\frac{1}{\epsilon_{\mathbf{k+q}\nu'}-\epsilon_{\mathbf{k}\nu}}\sum_{\mathbf{N,N}'}e^{i((\mathbf{k+q})\cdot\mathbf{R}_{\mathbf{N}'}-\mathbf{k}\cdot\mathbf{R}_{\mathbf{N}})} \cdot \\ \nonumber & \cdot \mathbf{c}_{s\mathbf{k+q}\nu'}\left(\mathbf{H}^{\mathbf{q}j}(|\mathbf{N}-\mathbf{N}'|)-\epsilon_{F} \mathbf{S}^{\mathbf{q}j}(|\mathbf{N}-\mathbf{N}'|) \right) \mathbf{c}_{s\mathbf{k}\nu}$

Derivation using first order time dependent perturbation theory¶

Our starting point is given by eqs. (3) and (5). The harmonic nuclear wavefunctions are perturbed by a potential $H_{\mathrm{ep}}$ , which we assume to be harmonic in time, namely with the frequency $\omega_{\mathbf{q}j}$ of the vibration itself. This can be understood by considering that at the equilibrium position of the nuclei the friction effects due to the electronic structure are zero, while they are strongest at the maximum displacement. The effect, as already discussed, is that electronic and vibrational states hybridize. Therefore all ‘vibronic’ states that can be excited by the phonon frequency $\omega_{\mathbf{q}j}$ will contribute to the decay of the phonon. As a result, in the following we will not distinguish between energies of electronic or vibrational states. The potential introduces mixing between the nuclear eigenfunctions $\chi$ characterized with coefficients $c_n(t)$ :

$\ket{\chi(t)} = \sum_n c_n(t) \operatorname{e}^{-iE_n t/\hbar} \ket{\chi_n}$

Expanding the coefficients and wavefunctions in a perturbation series and truncating after the first perturbation term, we can define the probability of excitation between two states (both differing by occupation change between phonons $\mathbf{q}j$ to $\mathbf{q}f$ ) by the absolute square of the first order coefficients:

$&P_{\mathbf{q}j\rightarrow \mathbf{q}f}(t)=|c^{(1)}_{jf}(t)|^2 = \\ \nonumber & =\frac{1}{\hbar^2}|\braket{\chi_{\mathbf{q}f}|H_{\mathrm{ep}}|\chi_{\mathbf{q}j}}|^2 \left( \frac{\sin((\Delta\omega-\omega_{\mathbf{q}j})t/2)}{(\Delta\omega-\omega_{\mathbf{q}j})} \right)^2 ,$

with $\Delta\omega$ as energy difference between all vibronic states generated by the perturbation. Due to our previous assumption of harmonic oscillator states every phonon state can only decay into the next lowest state with the result that the overall probability of finding state $\mathbf{q}j$ being the excitation probability into this state from the next lowest state

(14) $P_{\mathbf{q}j}(t) =|\braket{\chi_{\mathbf{q}f}|H_{\mathrm{ep}}|\chi_{\mathbf{q}i}}|^2 \left( \frac{\sin((\Delta\epsilon-\hbar\omega_{\mathbf{q}j})t/2\hbar)}{(\Delta\epsilon-\hbar\omega_{\mathbf{q}j})} \right)^2$

At this point we can introduce our previous result for the coupling matrix element by inserting eq. (5) into eq. (14)

(15) $P_{\mathbf{q}j}(t) = \frac{\hbar^3\omega_{\mathbf{q}j}}{2M}\sum_{s\mathbf{k}}w_{\mathbf{k}} \sum_{\nu<\epsilon_F}\sum_{\nu'>\epsilon_F}|g_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{\mathbf{q}j}|^2 \left( \frac{\sin((\Delta\epsilon-\hbar\omega_{\mathbf{q}j})t/2\hbar)}{(\Delta\epsilon-\hbar\omega_{\mathbf{q}j})} \right)^2$

We are interested in the decay rate in the limit of small perturbing frequencies or, equally, the limit of long time ( $t>>1/(\Delta\omega-\omega_{\mathbf{q}j})$ ). Taking this limit we can assume exponential decay and we arrive at the well known Fermi’s Golden Rule expression:

$\Gamma_{\mathbf{q}j} = \lim_{t\rightarrow\infty} \frac{P_{\mathbf{q}j}(t)}{t} = \frac{\pi\hbar^2\omega_{\mathbf{q}j}}{M}\sum_{s\mathbf{k}}w_{\mathbf{k}} \sum_{\nu<\epsilon_F}\sum_{\nu'>\epsilon_F}|g_{s\mathbf{k+q}\nu',s\mathbf{k}\nu}^{\mathbf{q}j}|^2 \delta(\Delta\epsilon-\hbar\omega_{\mathbf{q}j})$

However, even without taking the limit of $t\rightarrow\infty$ , applying the $\sin$ -function in eq. (15) with a reasonably chosen time would yield a similarly narrow integration kernel for the typical range of vibrational frequencies found for atomic or molecular adsorbates. Finally, we can generalize to non-integer finite-temperature occupations and arrive at:

(16) $\Gamma_{\mathbf{q}j} = \frac{\pi\hbar^2\omega_{\mathbf{q}j}}{M}\sum_{s\mathbf{k}}\sum_{\nu}\sum_{\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 \delta(\epsilon_{s\mathbf{k+q}\nu'}-\epsilon_{s\mathbf{k}\nu}-\hbar\omega_{\mathbf{q}j}) .$

With $\hbar\omega_{\mathbf{q}j}$ being the perturbing energy for which we apply the delta function in (16) we can re-express the above equation as

$\Gamma_{\mathbf{q}j} = \frac{\pi\hbar}{M}\sum_{s\mathbf{k}}\sum_{\nu}\sum_{\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}) .$