Light-Matter Interaction and Dipole Transition Matrix

Quantum theory of light-matter interaction

Let us consider an atom in the presence of an external classical electromagnetic fields, the gauge invariant Schrödinger equation writes ¹

\[\begin{equation} i\hbar\frac{\partial\psi(\mathbf{r}, t)}{\partial t} = \hat{H}\,\psi(\mathbf{r}, t) \label{eq:gauge_inv_schro_eq} \end{equation}\]

where the gauge invariant Hamiltonian

\[\begin{align} \hat{H} &= \frac{1}{2m} \left( {\mathbf{p} - \frac{e}{c}\mathbf{A}} \right)^2 + e\varphi(\mathbf{r}) + \hat{V}(\mathbf{r}) \nonumber \\[6pt] &= \frac{\mathbf{p}^2}{2m} + \hat{V}(\mathbf{r}) - \frac{e}{2mc}\left[ \mathbf{p}\cdot\mathbf{A}+\mathbf{A}\cdot\mathbf{p} \right] + \frac{e^2}{2mc^2}\mathbf{A}^2 + e\varphi(\mathbf{r}) \label{eq:gauge_inv_hamil} \end{align}\]

$\mathbf{A}$ and $\varphi$ are magnetic vector potential and electric scalar potential, recpetively. They are related to the electric and magnetic field by the following equations ²

\[\begin{align} \mathbf{E} &= -\nabla\varphi - \frac{1}{c}\partial_t\mathbf{A} \\[3pt] \mathbf{B} &= \nabla \times \mathbf{A} \end{align}\]

The form of Eq.\eqref{eq:gauge_inv_schro_eq} is invariant under the gauge transformation

\[\begin{align} \mathbf{A}' &= \mathbf{A} + \nabla\lambda(\mathbf{r}, t) \label{eq:gt_vec} \\[9pt] \mathbf{\varphi}' &= \varphi - \frac{1}{c}\frac{\partial\lambda(\mathbf{r}, t)}{\partial t} \label{eq:gt_sca} \\[6pt] \psi(\mathbf{r},t)' &= \exp\left[\frac{ie}{c\hbar}\lambda(\mathbf{r}, t)\right]\psi(\mathbf{r},t) \label{eq:gt_psi} \end{align}\]

i.e.

\[\begin{equation} i\hbar\frac{\partial\psi'(\mathbf{r}, t)}{\partial t} = \left[ \frac{1}{2m} \left( {\mathbf{p} - \frac{e}{c}\mathbf{A}'} \right)^2 + e\varphi' + \hat{V}(\mathbf{r}) \right] \psi'(\mathbf{r}, t) \end{equation}\]

Electric Dipole Approximation

Generally, the wavelength of the electromagnetic field is much larger than the typical size of an atom/molecule, i.e. $2\pi/k = \lambda \gg a$ or $ka \ll 1$ ($a$ is the typical atomic orbital size), we can thus assume that the vector potential is practically position independent for an atom/molecule, i.e.

\[\begin{equation} \mathbf{A}(\mathbf{r}, t) \approx \mathbf{A}(t) \end{equation}\]

The approximation is referred to as the electric dipole approximation (EDA). ³

More specifically, for a monochromatic field of plane wave

\[\begin{equation} \mathbf{A}(\mathbf{r}, t) = 2\mathbf{A}_0 \cos(\mathbf{k}\cdot\mathbf{r} - \omega t) \label{eq:mono_vecpot} \end{equation}\]

Under the EDA, we can thus take $e^{i\mathbf{k}\cdot\mathbf{r}} \approx 1$ with $\mathbf{k}\cdot\mathbf{r} \ll 1$, and as a result the vector potential does not depend on position.

\(\begin{equation} \mathbf{A}(t) = \mathbf{A}_0 [e^{i\omega t} + e^{-i\omega t}] \end{equation}\)

Coulomb Gauge

The Coulomb gauge (also known as the transverse gauge) is defined by the gauge condition

\[\begin{equation} \label{eq:coulomb_gauge} \nabla\cdot\mathbf{A}_C(\mathbf{r}, t) = 0 \end{equation}\]

where “C” in the subscript is used to indicate the gauge. Moreover, if there is no current in the space, then we also have

\[\begin{equation} \varphi_C(\mathbf{r}) = 0 \end{equation}\]

Generally, the operator $\mathbf{p}$ and $\mathbf{A}$ do not commute,

\[\begin{equation} \mathbf{p}\cdot\mathbf{A} - \mathbf{A}\cdot\mathbf{p} = \sum_j \left[ \hat{p}_j A_j - A_j\hat{p}_j \right] = \sum_j [\hat{p}_j, A_j] = -i\hbar \sum_j \partial_j A_j = -i\hbar \nabla\cdot\mathbf{A} \end{equation}\]

Note the difference between $\mathbf{p} \cdot \mathbf{A}$ and $-i\hbar\nabla \cdot \mathbf{A}$ in the above equation, where the former means two operators applying consecutively on the wavefunction and the latter is just the gradient of $\mathbf{A}$, i.e. the nabla operator applies on $\mathbf{A}$.

Apparently, under the Coulomb gauge, the two operators commute, i.e.

\[\begin{equation} \mathbf{p}\cdot\mathbf{A}_C = \mathbf{A}_C\cdot\mathbf{p} \end{equation}\]

Subsituting the Coulomb gauge condition into Eq.\eqref{eq:mono_vecpot}, one has

\[\mathbf{k} \cdot \mathbf{A}_0 = 0\]

i.e. the field propagation direction $\mathbf{k}$ is perpendicular to the field vibration direction $\mathbf{A}_0$, hence the name transverse gauge.

Length Gauge

Starting with the Coulomb gauge and electric dipole approximation (vector potential is position independent and only depends on time), and by choosing $\lambda(\mathbf{r}, t)$ in the gauge transformation Eq.\eqref{eq:gt_vec}, Eq.\eqref{eq:gt_sca} and Eq.\eqref{eq:gt_psi} as

\[\begin{equation} \lambda_L(\mathbf{r}, t) = - \mathbf{r}\cdot \mathbf{A}_C(t) \end{equation}\]

One has the two potentials in length gauge (LG)

\[\begin{align} \mathbf{A}_L &= \mathbf{A}_C(t) + \nabla\lambda_L(\mathbf{r}, t) = 0 \label{eq:LG_vec} \\[9pt] \mathbf{\varphi}_L &= \varphi_C - \frac{1}{c}\frac{\partial\lambda_L(\mathbf{r}, t)}{\partial t} = -\mathbf{r}\cdot\boldsymbol{\mathcal E}(t) \label{eq:LG_sca} \end{align}\]

where $\boldsymbol{\mathcal E}(t) = -c^{-1}\partial_t \mathbf{A}_{C}(t)$ and we have used “C” and “L” in the subscripts to indicate the Coulomb and length gauge, respectively.

Substituting Eq.\eqref{eq:LG_vec} and Eq.\eqref{eq:LG_sca} into Eq.\eqref{eq:gauge_inv_hamil}, one obtains the Hamiltonian in length gauge

\[\begin{align} \hat{H}_L &= \hat{H}_0 + \hat{H}_L' \label{eq:hamil_lg} \\ \hat{H}_0 &= \frac{\mathbf{p}^2}{2m} + \hat{V}(\mathbf{r}) \label{eq:hamil_free_lg} \\ \hat{H}_L' &= - e\mathbf{r} \cdot \boldsymbol{\mathcal E}(t) \label{eq:hamil_I_lg} \end{align}\]

where $\hat{H}_I’$ is the light-matter interaction in the length gauge.

Velocity Gauge

Again, we start with the Coulomb gauge and EDA. However, this time we choose $\lambda$ as

\[\begin{equation} \lambda_V(t) = \frac{e}{2mc} \int_{-\infty}^t A_C^2(\tau) \, \mathrm{d}\tau \end{equation}\]

This leads to the two potentials in velocity gauge (VG)

\[\begin{align} \mathbf{A}_V &= \mathbf{A}_C(t) + \nabla\lambda_V(t) = \mathbf{A}_C \label{eq:VG_vec} \\[9pt] \mathbf{\varphi}_V &= \varphi_C - \frac{1}{c}\frac{\partial\lambda_V(\mathbf{r}, t)}{\partial t} = -\frac{e^2}{2mc^2} \mathbf{A}_C(t) \label{eq:VG_sca} \end{align}\]

Substituting Eq.\eqref{eq:VG_vec} and Eq.\eqref{eq:VG_sca} into Eq.\eqref{eq:gauge_inv_hamil}, one obtains the Hamiltonian in the velocity gauge

\[\begin{align} \hat{H}_L &= \hat{H}_0 + \hat{H}_V' \label{eq:hamil_vg} \\[6pt] \hat{H}_V' &= - \frac{e}{mc} \mathbf{A} \cdot \mathbf{p} \label{eq:hamil_I_vg} \end{align}\]

with $\hat{H}_V’$ being the light-matter interaction in the velocity gauge.

Dipole and Momentum Transition Matrix

According to Fermi’s golden rule ⁴, the transition rate between two states, in the length gauge, is given by

\[\begin{equation} P_{i\to j} = \frac{2\pi}{\hbar} \left|\langle \psi_i |\hat{H}_L' | \psi_j \rangle \right|^2 \delta(E_i - E_j \pm \hbar\omega) \propto \left| \langle \psi_i | e \mathbf{r} | \psi_j \rangle \right|^2 \label{eq:rate_lg} \end{equation}\]

where the matrix element of $e\mathbf{r}$ is usually referred to as the dipole matrix element or transition dipole moment (TDM).

\[\begin{equation} \mathbf{T}_{ij} = e\langle \psi_i | \mathbf{r} | \psi_j \rangle \end{equation}\]

Or in the velocity gauge

\[\begin{equation} P_{i\to j} = \frac{2\pi}{\hbar} \left|\langle \psi_i |\hat{H}_V' | \psi_j \rangle \right|^2 \delta(E_i - E_j \pm \hbar\omega) \propto \left| \langle \psi_i | \mathbf{p} | \psi_j \rangle \right|^2 \label{eq:rate_vg} \end{equation}\]

If the molecular orbitals are eigenfunctions of some Hamiltonian $\hat{H}_0$, i.e.

\[\begin{equation} \hat{H}_0\cdot\psi_j(\mathbf{r}) = \left[ \frac{\mathbf{p}^2}{2m} + \hat{V}(\mathbf{r}) \right] \psi_j(\mathbf{r}) = \varepsilon_j\, \psi_j(\mathbf{r}) \end{equation}\]

Then one can utilize the commutator relation

\[\begin{gather} [\mathbf{r}, \hat{H}_0] = i\hbar\frac{\partial\hat{H}_0}{\partial \mathbf{p}} = i\hbar\frac{\mathbf{p}}{m} \\[6pt] \langle \psi_i | [\mathbf{r}, \hat{H}_0] | \psi_j \rangle = (\varepsilon_j - \varepsilon_i)\cdot \mathbf{T}_{ij} = \frac{i\hbar}{m} \langle \psi_i | \mathbf{p} | \psi_j \rangle \end{gather}\]

which lead to

\[\begin{equation} \mathbf{T}_{ij} = e\langle \psi_i | \mathbf{r} |\psi_j \rangle = \frac{i\hbar e} {m(\varepsilon_j - \varepsilon_i)} \, \langle \psi_i | \mathbf{p} |\psi_j \rangle \label{eq:p-r_relation} \end{equation}\]

We will shorten the above relation between the momentum matrix element and the transition dipole moment as the “p-r” relation.

Transition between Molecular Orbitals

The transition dipole moment between two molecular orbitals

\[\begin{align} \mathbf{T}_{ij} = e\langle \psi_i | \mathbf{r} | \psi_j \rangle &= e\int \psi_i(\mathbf{r})^\star \,\mathbf{r}\, \psi_j(\mathbf{r}) \,\mathrm{d}\mathbf{r} \\[6pt] &= e\begin{pmatrix} \displaystyle \int \psi_i(\mathbf{r})^\star \,x\, \psi_j(\mathbf{r}) \,\mathrm{d}\mathbf{r} \\[3pt] \displaystyle \int \psi_i(\mathbf{r})^\star \,y\, \psi_j(\mathbf{r}) \,\mathrm{d}\mathbf{r} \\[3pt] \displaystyle \int \psi_i(\mathbf{r})^\star \,z\, \psi_j(\mathbf{r}) \,\mathrm{d}\mathbf{r} \end{pmatrix} \label{eq:tdm_def_r_vec} \end{align}\]

where $\psi(\mathbf{r})$ can generally be expanded, within a certain accuracy, with a finite set of spherical harmonics (multipole expansion), i.e. ⁵

\[\begin{equation} \psi(\mathbf{r}) = \sum_{l=0}^L \sum_{m=-l}^l R_{l}(r)\,Y_{l,m}(\hat{\mathbf{r}}) \end{equation}\]

Both real and complex spherical harmonics can be used in the expansion, however, since the molecular orbitals can be real functions, real spherical harmonics are usually used in calculations to save memory and improve the speed. For simplicity, we will just focus on the case where there is only one term in the expansion, i.e.

\[\begin{equation} \psi(\mathbf{r}) = R_{l}(r)\,Y_{l,m}(\hat{\mathbf{r}}) \label{eq:orb_rad_ang} \end{equation}\]

Recall that the real spherical harmonics $Y_{lm}$ of quantum number $l=1$ are defined as ⁶

\[\begin{align*} Y_{1,-1} &= i\sqrt{1\over2}(Y_1^{-1} + Y_1^{1}) = \sqrt{3\over4\pi}\cdot{y\over r} \\[6pt] Y_{1,0} &= Y_1^0 = \sqrt{3\over4\pi}\cdot{z\over r} \\[6pt] Y_{1,1} &= \sqrt{1\over2}(Y_1^{-1} + Y_1^{1}) = \sqrt{3\over4\pi}\cdot{x\over r} \end{align*}\]

or equivalently

\[\begin{equation} x = \sqrt{\frac{4\pi}{3}} r \cdot Y_{1,1};\qquad y = \sqrt{\frac{4\pi}{3}} r \cdot Y_{1,-1};\qquad z = \sqrt{\frac{4\pi}{3}} r \cdot Y_{1,0} \label{eq:real_sph_l1_f2} \end{equation}\]

Substituting Eq.\eqref{eq:real_sph_l1_f2} and Eq.\eqref{eq:orb_rad_ang} into Eq.\eqref{eq:tdm_def_r_vec}, one have

\[\begin{equation} \mathbf{T}_{ij} = e\left[ \int_0^\infty R_{l_1}(r) R_{l_2}(r) r^3 \,\mathrm{d}r \right] \cdot \sqrt{\frac{4\pi}{3}} \begin{pmatrix} {\cal G}(l_1, l_2, 1, m_1, m_2, 1) \\[6pt] {\cal G}(l_1, l_2, 1, m_1, m_2,-1) \\[6pt] {\cal G}(l_1, l_2, 1, m_1, m_2, 0) \\ \end{pmatrix} \label{eq:tdm_real_space} \end{equation}\]

where ${\cal G}(l_1, l_2, l_3, m_1, m_2, m_3)$ is the real Gaunt coefficients defined as

\[\begin{equation} {\cal G}(l_1, l_2, l_3, m_1, m_2, m_3) = \int\limits_0^\pi \sin\theta \mathrm{d}\theta \int\limits_0^{2\pi} \mathrm{d}\phi \, Y_{l_1,m_1}(\hat{\mathbf{r}}) \, Y_{l_2,m_2}(\hat{\mathbf{r}}) \, Y_{l_3,m_3}(\hat{\mathbf{r}}) \end{equation}\]

and can be evaluated analytically. ⁷ Note that $r^3$ rather than $r^2$ is used in the radial integration, where the extra $r$ is introduced by Eq.\eqref{eq:real_sph_l1_f2}. Eq.\eqref{eq:tdm_real_space} shows that the TDM between molecular orbitals can be written as a radial integration term multiplied by an analytical angular term.

TDM from p-r relation

The matrix element of the momentum operator can be easily evaluated in the momentum space, i.e.

\[\begin{align} \langle \psi_i | \mathbf{p} | \psi_j \rangle &= \int \langle \psi_i | \mathbf{k} \rangle \cdot \hbar\mathbf{k} \cdot \langle \mathbf{k} | \psi_j \rangle \,\mathrm{d}\mathbf{k} \\[6pt] &= \hbar \begin{pmatrix} \displaystyle \int \psi_i(\mathbf{k})^\star \,k_x\, \psi_j(\mathbf{k}) \,\mathrm{d}\mathbf{k} \\[3pt] \displaystyle \int \psi_i(\mathbf{k})^\star \,k_y\, \psi_j(\mathbf{k}) \,\mathrm{d}\mathbf{k} \\[3pt] \displaystyle \int \psi_i(\mathbf{k})^\star \,k_z\, \psi_j(\mathbf{k}) \,\mathrm{d}\mathbf{k} \end{pmatrix} \end{align}\]

where the wavefunction in momentum space can be written as

\[\begin{equation} \psi(\mathbf{k}) = \langle \mathbf{k} | \psi \rangle = i^l G_l(k) \, Y_{l,m}(\hat{\mathbf{k}}) \label{eq:orb_momentum_space} \end{equation}\]

with $G_l(k)$ being the spherical Bessel transform (SBT) of the radial function $R_l(r)$

\[\begin{equation} G_l(k) = \sqrt{\frac{2}{\pi}} \int_0^\infty R_l(r)\,j_l(kr) r^2 \,\mathrm{d}r \label{eq:orb_sbt} \end{equation}\]

Combining Eq.\eqref{eq:orb_momentum_space}, Eq.\eqref{eq:orb_sbt} and Eq.\eqref{eq:p-r_relation}, one have the TDM evaluated in momentum space

\[\begin{align} \mathbf{T}_{ij} &= \frac{ i^{l_2 - l_1+1}\hbar^2e }{ m(\varepsilon_j - \varepsilon_i) } \, \left[ \int_0^\infty G_{l_1}(k) G_{l_2}(k) k^3 \,\mathrm{d}k \right] \nonumber\\[9pt] &\quad\times \sqrt{\frac{4\pi}{3}} \begin{pmatrix} {\cal G}(l_1, l_2, 1, m_1, m_2, 1) \\[6pt] {\cal G}(l_1, l_2, 1, m_1, m_2,-1) \\[6pt] {\cal G}(l_1, l_2, 1, m_1, m_2, 0) \\ \end{pmatrix} \label{eq:tdm_momentum_space_final} \end{align}\]

Again, the TDM from momentum space evaluation can be written as a radial integration over $k$, multiplied by the same analytical angular term as in Eq.\eqref{eq:tdm_real_space}.

Selection Rules

• If real molecular orbitals are used, the TDM should also be real vectors. To ensure this, one require that $l_2 - l_1 + 1 = 2n$, where $n\in\mathbb{N}$ in Eq.\eqref{eq:tdm_momentum_space_final}.

• The properties of the real Gaunt coefficients require that $l_1 + l_2 + 1=2n$, where $n\in\mathbb{N}$. ⁸

• The spherical harmonics have definite parity, i.e. they are either even or odd with respect to inversion about the origin ⁹

  \[\begin{equation*} Y_{l,m}(\hat{\mathbf{r}}) = (-1)^l \, Y_{l,m}(-\hat{\mathbf{r}}) \end{equation*}\]

  Moreover, $\mathbf{r}$ is an odd function. To make sure that TDM is non-zero, one require that the $l_1 + l_2$ is ODD so that the integration is not zero.

  If complex spherical harmonics are used, then in this case

  \[\begin{equation} x = \sqrt{\frac{2\pi}{3}} r \cdot [Y_1^1 + Y_1^{-1}];\qquad y = \sqrt{\frac{2\pi}{3}} r \cdot i[Y_1^{-1} + Y_1^1];\qquad z = \sqrt{\frac{4\pi}{3}} r \cdot Y_1^0 \label{eq:complex_sph_l1_f2} \end{equation}\]

  The Gaunt coefficients defined by complex spherical harmonics are only non-zero for $m_1 + m_2 + m_3 = 0$. ¹⁰ Therefore, the TDM is only non-zero for

  \[\begin{equation} \Delta m = m_1 - m_2 = 0\ \mathrm{or}\ \pm 1 \end{equation}\]

TDM between hydrogen orbitals

Below, we use hydrogen orbital, which can be obtained analytically, as a simple example calculate the TDM and verify the equivalence of the two method. We choose the initial orbital as hydrogen 1s and the final state being hydrogen 2p_x orbitals.

The script follows, where atomics units are used

From the outputs, one can see that transition dipole moment from the two methods are indeed identical.

TDM between Bloch states

For crystals, the wavefunctions take the form of a Bloch wave

\[\begin{equation} | \psi_{n\mathbf{k}} \rangle = e^{i\mathbf{k}\cdot\mathbf{r}} | u_{n\mathbf{k}} \rangle \end{equation}\]

For finite crystals with periodic boundary condition (PBC), the matrix elements of the momentum operator is related to the gradient-k operator by ¹¹

\[\begin{equation} \langle \psi_{n\mathbf{k}} | \mathbf{p} | \psi_{m\mathbf{k}} \rangle = \frac{ im (\varepsilon_{n\mathbf{k}} - \varepsilon_{m\mathbf{k}}) }{\hbar} \langle u_{n\mathbf{k}} | i\nabla_\mathbf{k} | u_{m\mathbf{k}} \rangle \end{equation}\]

Moreover, there’s no simple p-r relation between momentum and position operator matrix elements.

\[\begin{equation} \langle \psi_{n\mathbf{k}} | \mathbf{r} | \psi_{m\mathbf{k}} \rangle \ne \langle u_{n\mathbf{k}} | i\nabla_\mathbf{k} | u_{m\mathbf{k}} \rangle \end{equation}\]

We will focus on the matrix elements of the momentum opeartor in the Bloch waves. Within the PAW formalism, the momentum opeartor matrix element is given by

\[\begin{align} \langle \psi_{c\mathbf{k}} | \mathbf{p} | \psi_{v\mathbf{k}} \rangle &= \hbar \langle u_{c,\mathbf{k}} | \mathbf{k} - i\nabla | u_{v,\mathbf{k}} \rangle \\[9pt] \langle u_{c,\mathbf{k}} | \mathbf{k} - i\nabla | u_{v,\mathbf{k}} \rangle &= \langle \tilde{u}_{c,\mathbf{k}} | \mathbf{k} - i\nabla | \tilde{u}_{v,\mathbf{k}} \rangle - i \sum_{nm} \langle \tilde{u}_{c,\mathbf{k}} | \tilde{p}_{n,\mathbf{k}} \rangle \langle \tilde{p}_{m,\mathbf{k}} | \tilde{u}_{v,\mathbf{k}} \rangle \nonumber \\[6pt] &\quad \times \left[ \langle \phi_n | \nabla_\mathbf{r} | \phi_m \rangle - \langle \tilde\phi_n | \nabla_\mathbf{r} | \tilde\phi_m \rangle \right] \end{align}\]

where $\tilde{p}_{n\mathbf{k}}$ is the k-dependent projector function as defined in a previous post, $\phi$ and $\tilde\phi$ are the PAW AE and PS partial waves.

Note that the imaginary part of the dielectric function in the longitudinal expression is given by ¹²

\[\begin{align} \epsilon''_{\alpha\beta}(\omega) &= \frac{4\pi^2 e^2}{\Omega} \lim_{\mathbf{q}\to0} \frac{1}{q^2} \sum_{c,v,\mathbf{k}} 2\mathrm{w}_\mathbf{k} \delta(\varepsilon_{c,\mathbf{k}} - \varepsilon_{v,\mathbf{k}} - \omega) \nonumber\\[6pt] &\qquad \times \langle u_{c,\mathbf{k}+q\mathbf{e}_\alpha} | u_{v,\mathbf{k}} \rangle \langle u_{c,\mathbf{k}} | u_{v,\mathbf{k}+q\mathbf{e}_\beta} \rangle^* \end{align}\]

and in the transversal expression by ¹²

\[\begin{align} \epsilon''_{\alpha\beta}(\omega) &= \frac{4\pi^2 e^2 \hbar^4}{\Omega \omega^2 m_e^2} \lim_{\mathbf{q}\to0} \sum_{c,v,\mathbf{k}} 2\mathrm{w}_\mathbf{k} \delta(\varepsilon_{c,\mathbf{k}+\mathbf{q}} - \varepsilon_{v,\mathbf{k}} - \omega) \nonumber\\[6pt] &\qquad \times \langle u_{c,\mathbf{k}} | i\nabla_\alpha - \mathbf{k}_\alpha | u_{v,\mathbf{k}} \rangle \langle u_{c,\mathbf{k}} | i\nabla_\beta - \mathbf{k}_\beta | u_{v,\mathbf{k}} \rangle^* \end{align}\]

In principle, both longitudinal and transversal expression yield almost exactly indetical results in the limit of a complete PAW one-center basis. However, this limit is in practice not always reached, and hence the longitudinal expression should be preferred over the simpler transversal expression.

Matrix element of nabla operator — pseudo-wavefunction

The matrix elements of the nabla operator can be readily evaluated in reciprocal space

\[\begin{equation} \langle u_{c,\mathbf{k}} | \mathbf{k} - i\nabla | u_{v,\mathbf{k}} \rangle = \sum_{\mathbf{G}} c_{c,\mathbf{k}}^*(\mathbf{G}) \cdot \left[ \mathbf{k} + \mathbf{G} \right] \cdot c_{v,\mathbf{k}} (\mathbf{G}) \end{equation}\]

where $c(\mathbf{G})$ is the plane-wave coefficients for $u_{n\mathbf{k}}(\mathbf{r})$.

Matrix element of nabla operator — one-center contribution

\[\begin{equation} \boldsymbol\tau_{nm} = \langle \phi_n | \nabla_\mathbf{r} | \phi_m \rangle - \langle \tilde\phi_n | \nabla_\mathbf{r} | \tilde\phi_m \rangle \end{equation}\]

Recall that AE/PS partial waves consist of a radial part and an angular part, ¹³

\[\begin{align} \phi_n(\mathbf{r}) &= U_{n,l}(r) \cdot Y_{l,m}(\hat{\mathbf{r}}) \\ \tilde\phi_n(\mathbf{r}) &= \tilde{U}_{n,l}(r) \cdot Y_{l,m}(\hat{\mathbf{r}}) \end{align}\]

and their Fourier transform

\[\begin{align} \phi_n(\mathbf{k}) &= i^l G_{n,l}(k) \cdot Y_{l,m}(\hat{\mathbf{k}}) \\ \tilde\phi_n(\mathbf{k}) &= i^l \tilde{G}_{n,l}(k) \cdot Y_{l,m}(\hat{\mathbf{k}}) \end{align}\]

where $G_{n,l}(k)$ is the spherical Bessel transform of $U_{n,l}(r)$. ¹⁴

The matrix element of the nabla operator can be evaluated by either real-space or momentum-space integration.

Real-space integration

\[\begin{align} \langle \phi_n | \nabla_\mathbf{r} | \phi_m \rangle &= \int \phi_n(\mathbf{r}) \nabla_\mathbf{r} \phi_m(\mathbf{r}) \,\mathrm{d}\mathbf{r} \nonumber\\ &= \int U_{n,l'}(r) Y_{l',m}(\hat{\mathbf{r}}) \, \nabla_\mathbf{r} \left[ \frac{U_{n,l}(r)}{r^{l}} \cdot r^{l} Y_{l,m}(\hat{\mathbf{r}}) \right] \,r^2 \mathrm{d}r \mathrm{d}\Omega \nonumber\\ &= \int U_{n,l'}(r) r^{l+2} \, \frac{\partial}{\partial r} \, \left[ \frac{U_{n,l}(r)}{r^{l}} \right] Y_{l',m}(\hat{\mathbf{r}}) Y_{l,m}(\hat{\mathbf{r}}) \, \vec{e}_\mathbf{r} \, \mathrm{d}r \mathrm{d}\Omega \nonumber\\ &\quad + \int U_{n,l'}(r) U_{n,l}(r) r^{2-l} \, Y_{l',m}(\hat{\mathbf{r}}) \nabla_\mathbf{r}[r^{l} Y_{l,m}(\hat{\mathbf{r}})] \, \mathrm{d}r \mathrm{d}\Omega \label{eq:nabla_real_0} \end{align}\]

• For the first term in the last line of Eq.\eqref{eq:nabla_real_0}, remember that

  \[\begin{equation} \vec{e}_\mathbf{r} = \begin{pmatrix} \displaystyle \frac{x}{r}, & \displaystyle \frac{y}{r}, & \displaystyle \frac{z}{r} \end{pmatrix} \end{equation}\]

  which can be written as real spherical harmonics with $l=1$, i.e.

  \[\begin{equation} \vec{e}_\mathbf{r} = \sqrt{\frac{4\pi}{3}} \begin{pmatrix} Y_{1, 1}(\hat{\mathbf{r}}), & Y_{1,-1}(\hat{\mathbf{r}}),& Y_{1, 0}(\hat{\mathbf{r}}) \end{pmatrix} \end{equation}\]

  As a result, the first term reduces to

  \[\begin{equation} \sqrt{\frac{4\pi}{3}} \int_0^\infty U_{n,l'}(r) r^{l+2} \, \frac{\partial}{\partial r} \, \left[ \frac{U_{n,l}(r)}{r^{l}} \right] \mathrm{d}r \times \begin{pmatrix} {\cal G}(l', l, 1, m', m, 1) \\[6pt] {\cal G}(l', l, 1, m', m,-1) \\[6pt] {\cal G}(l', l, 1, m', m, 0) \\ \end{pmatrix} \label{eq:nabla_real_1} \end{equation}\]

• For the second term in Eq.\eqref{eq:nabla_real_0}, one can show that ⁶

  \[\begin{equation} \nabla_\mathbf{r} [r^l Y_{l,m}(\hat{\mathbf{r}})] = r^{l-1} \sum_{m=1-l}^{l-1} c_m Y_{l-1,m}(\hat{\mathbf{r}}) \end{equation}\]

  Therefore, the second term reduces to

  \[\begin{equation} \int_0^\infty U_{n,l'}(r) U_{n,l}(r) r \, \mathrm{d}r \times \int Y_{l',m}(\hat{\mathbf{r}}) r^{1-l} \nabla_\mathbf{r}[r^{l} Y_{l,m}(\hat{\mathbf{r}})] \, \mathrm{d}\Omega \label{eq:nabla_real_2} \end{equation}\]

  The integration with respect to the angle can be performed using GPAW. ¹⁵

Note that in VASP, $r \, U(r) = R(r)$ instead of $U(r)$ for partial waves are stored in POTCAR. By taking this into account and combine Eq.\eqref{eq:nabla_real_1} and Eq.\eqref{eq:nabla_real_2}, one have

\[\begin{align} \langle \phi_n | \nabla_\mathbf{r} | \phi_m \rangle &= \sqrt{\frac{4\pi}{3}} \int_0^\infty \left[ R_{n,l'} \frac{\partial R_{n,l}}{\partial r} - (l+1) \frac{R_{n,l'} R_{n,l}}{r} \right] \mathrm{d}r \times \begin{pmatrix} {\cal G}(l', l, 1, m', m, 1) \\[6pt] {\cal G}(l', l, 1, m', m,-1) \\[6pt] {\cal G}(l', l, 1, m', m, 0) \\ \end{pmatrix} \nonumber \\[9pt] &\quad + \int_0^\infty \frac{R_{n,l'} R_{n,l}}{r} \, \mathrm{d}r \times \int Y_{l',m}(\hat{\mathbf{r}}) r^{1-l} \nabla_\mathbf{r}[r^{l} Y_{l,m}(\hat{\mathbf{r}})] \, \mathrm{d}\Omega \label{eq:nabla_real_final} \end{align}\]

All the terms with spherical harmonics in Eq.\eqref{eq:nabla_real_final} can be obtained analytically.

The subroutine NABLA_RADIAL in VASP source file optics.F calculates $\boldsymbol\tau_{nm}$ by real space integration. The results are stored in NABLA of a potcar type. VASP version 5.x caluclate the matrix element up to lmax = 2 while version 6.x increases lmax to 3.

Since the AE/PS partial waves are not necessarily zero at $r=r_\mathrm{max}$, the matrix element of nabla operator is asymmetric with respect to AE and PS partial waves. However, the difference of the two is symmetric.

One can add the following code to main.F of VASP source and perform a recompile. Afterwards, VASP can outputs $\boldsymbol\tau_{nm}$ to file “nabla_X”, where X is the name of the elements.

Momentum-space integration

Similar to the evaluaton of Eq.\eqref{eq:tdm_momentum_space_final}

\[\begin{align} \langle \phi_n | \nabla_\mathbf{r} | \phi_m \rangle &= i^{l_2 - l_1+1} \, \int_0^\infty G_{n,l_1}(k) G_{m,l_2}(k) k^3 \,\mathrm{d}k \nonumber \\[9pt] &\quad\times \sqrt{\frac{4\pi}{3}} \begin{pmatrix} {\cal G}(l_1, l_2, 1, m_1, m_2, 1) \\[6pt] {\cal G}(l_1, l_2, 1, m_1, m_2,-1) \\[6pt] {\cal G}(l_1, l_2, 1, m_1, m_2, 0) \\ \end{pmatrix} \end{align}\]

Therefore

\[\begin{align} \boldsymbol\tau_{nm} &= \langle \phi_n | \nabla_\mathbf{r} | \phi_m \rangle - \langle \tilde\phi_n | \nabla_\mathbf{r} | \tilde\phi_m \rangle \nonumber\\[9pt] &= i^{l_2 - l_1+1} \, \int_0^\infty \left[ G_{n,l_1}(k) G_{m,l_2}(k) - \tilde{G}_{n,l_1}(k) \tilde{G}_{m,l_2}(k) \right] k^3 \,\mathrm{d}k \nonumber \\[9pt] &\quad\times \sqrt{\frac{4\pi}{3}} \begin{pmatrix} {\cal G}(l_1, l_2, 1, m_1, m_2, 1) \\[6pt] {\cal G}(l_1, l_2, 1, m_1, m_2,-1) \\[6pt] {\cal G}(l_1, l_2, 1, m_1, m_2, 0) \\ \end{pmatrix} \end{align}\]

The matrix element of the nabla operator from real/momentum space integration is implemented in the VaspBandUnfolding package.

The outputs of the script are shown below, where the numbers in the patentheses of the first column and first line of each section correspond to (n, l, m) for AE/PS partial waves.

Reference