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{matrix} (1) & i ℏ \frac{\partial ψ (r, t)}{\partial t} = \hat{H} ψ (r, t) \end{matrix}

where the gauge invariant Hamiltonian

\begin{aligned} \hat{H} & = \frac{1}{2 m} {(p - \frac{e}{c} A)}^{2} + e φ (r) + \hat{V} (r) \\ (2) & = \frac{p^{2}}{2 m} + \hat{V} (r) - \frac{e}{2 m c} [p \cdot A + A \cdot p] + \frac{e^{2}}{2 m c^{2}} A^{2} + e φ (r) \end{aligned}

$A$ and $φ$ are magnetic vector potential and electric scalar potential, recpetively. They are related to the electric and magnetic field by the following equations ²

\begin{aligned} (3) & E & = - \nabla φ - \frac{1}{c} \partial_{t} A \\ (4) & B & = \nabla \times A \end{aligned}

The form of Eq. $(1)$ is invariant under the gauge transformation

\begin{aligned} (5) & A^{'} & = A + \nabla λ (r, t) \\ (6) & φ^{'} & = φ - \frac{1}{c} \frac{\partial λ (r, t)}{\partial t} \\ (7) & ψ (r, t)^{'} & = \exp [\frac{i e}{c ℏ} λ (r, t)] ψ (r, t) \end{aligned}

i.e.

\begin{matrix} (8) & i ℏ \frac{\partial ψ^{'} (r, t)}{\partial t} = [\frac{1}{2 m} {(p - \frac{e}{c} A^{'})}^{2} + e φ^{'} + \hat{V} (r)] ψ^{'} (r, t) \end{matrix}

Electric Dipole Approximation

Generally, the wavelength of the electromagnetic field is much larger than the typical size of an atom/molecule, i.e. $2 π / k = λ ≫ a$ or $k a ≪ 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{matrix} (9) & A (r, t) \approx A (t) \end{matrix}

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

More specifically, for a monochromatic field of plane wave

\begin{matrix} (10) & A (r, t) = 2 A_{0} \cos (k \cdot r - ω t) \end{matrix}

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

$\begin{matrix} (11) & A (t) = A_{0} [e^{i ω t} + e^{- i ω t}] \end{matrix}$

Coulomb Gauge

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

\begin{matrix} (12) & \nabla \cdot A_{C} (r, t) = 0 \end{matrix}

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{matrix} (13) & φ_{C} (r) = 0 \end{matrix}

Generally, the operator $p$ and $A$ do not commute,

\begin{matrix} (14) & p \cdot A - A \cdot p = \sum_{j} [{\hat{p}}_{j} A_{j} - A_{j} {\hat{p}}_{j}] = \sum_{j} [{\hat{p}}_{j}, A_{j}] = - i ℏ \sum_{j} \partial_{j} A_{j} = - i ℏ \nabla \cdot A \end{matrix}

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

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

\begin{matrix} (15) & p \cdot A_{C} = A_{C} \cdot p \end{matrix}

Subsituting the Coulomb gauge condition into Eq. $(10)$ , one has

k \cdot A_{0} = 0

i.e. the field propagation direction $k$ is perpendicular to the field vibration direction $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 $λ (r, t)$ in the gauge transformation Eq. $(5)$ , Eq. $(6)$ and Eq. $(7)$ as

\begin{matrix} (16) & λ_{L} (r, t) = - r \cdot A_{C} (t) \end{matrix}

One has the two potentials in length gauge (LG)

\begin{aligned} (17) & A_{L} & = A_{C} (t) + \nabla λ_{L} (r, t) = 0 \\ (18) & φ_{L} & = φ_{C} - \frac{1}{c} \frac{\partial λ_{L} (r, t)}{\partial t} = - r \cdot E (t) \end{aligned}

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

Substituting Eq. $(17)$ and Eq. $(18)$ into Eq. $(2)$ , one obtains the Hamiltonian in length gauge

\begin{aligned} (19) & {\hat{H}}_{L} & = {\hat{H}}_{0} + {\hat{H}}_{L}^{'} \\ (20) & {\hat{H}}_{0} & = \frac{p^{2}}{2 m} + \hat{V} (r) \\ (21) & {\hat{H}}_{L}^{'} & = - e r \cdot E (t) \end{aligned}

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 $λ$ as

\begin{matrix} (22) & λ_{V} (t) = \frac{e}{2 m c} \int_{- \infty}^{t} A_{C}^{2} (τ) d τ \end{matrix}

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

\begin{aligned} (23) & A_{V} & = A_{C} (t) + \nabla λ_{V} (t) = A_{C} \\ (24) & φ_{V} & = φ_{C} - \frac{1}{c} \frac{\partial λ_{V} (r, t)}{\partial t} = - \frac{e^{2}}{2 m c^{2}} A_{C} (t) \end{aligned}

Substituting Eq. $(23)$ and Eq. $(24)$ into Eq. $(2)$ , one obtains the Hamiltonian in the velocity gauge

\begin{aligned} (25) & {\hat{H}}_{L} & = {\hat{H}}_{0} + {\hat{H}}_{V}^{'} \\ (26) & {\hat{H}}_{V}^{'} & = - \frac{e}{m c} A \cdot p \end{aligned}

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{matrix} (27) & P_{i \to j} = \frac{2 π}{ℏ} {| ⟨ ψ_{i} | {\hat{H}}_{L}^{'} | ψ_{j} ⟩ |}^{2} δ (E_{i} - E_{j} \pm ℏ ω) \propto {| ⟨ ψ_{i} | e r | ψ_{j} ⟩ |}^{2} \end{matrix}

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

\begin{matrix} (28) & T_{i j} = e ⟨ ψ_{i} | r | ψ_{j} ⟩ \end{matrix}

Or in the velocity gauge

\begin{matrix} (29) & P_{i \to j} = \frac{2 π}{ℏ} {| ⟨ ψ_{i} | {\hat{H}}_{V}^{'} | ψ_{j} ⟩ |}^{2} δ (E_{i} - E_{j} \pm ℏ ω) \propto {| ⟨ ψ_{i} | p | ψ_{j} ⟩ |}^{2} \end{matrix}

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

\begin{matrix} (30) & {\hat{H}}_{0} \cdot ψ_{j} (r) = [\frac{p^{2}}{2 m} + \hat{V} (r)] ψ_{j} (r) = ε_{j} ψ_{j} (r) \end{matrix}

Then one can utilize the commutator relation

\begin{matrix} (31) & [r, {\hat{H}}_{0}] = i ℏ \frac{\partial {\hat{H}}_{0}}{\partial p} = i ℏ \frac{p}{m} \\ (32) & ⟨ ψ_{i} | [r, {\hat{H}}_{0}] | ψ_{j} ⟩ = (ε_{j} - ε_{i}) \cdot T_{i j} = \frac{i ℏ}{m} ⟨ ψ_{i} | p | ψ_{j} ⟩ \end{matrix}

which lead to

\begin{matrix} (33) & T_{i j} = e ⟨ ψ_{i} | r | ψ_{j} ⟩ = \frac{i ℏ e}{m (ε_{j} - ε_{i})} ⟨ ψ_{i} | p | ψ_{j} ⟩ \end{matrix}

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{aligned} (34) & T_{i j} = e ⟨ ψ_{i} | r | ψ_{j} ⟩ & = e \int ψ_{i} (r)^{⋆} r ψ_{j} (r) d r \\ (35) & = e (\begin{array}{c} \int ψ_{i} (r)^{⋆} x ψ_{j} (r) d r \\ \int ψ_{i} (r)^{⋆} y ψ_{j} (r) d r \\ \int ψ_{i} (r)^{⋆} z ψ_{j} (r) d r \end{array}) \end{aligned}

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

\begin{matrix} (36) & ψ (r) = \sum_{l = 0}^{L} \sum_{m = - l}^{l} R_{l} (r) Y_{l, m} (\hat{r}) \end{matrix}

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{matrix} (37) & ψ (r) = R_{l} (r) Y_{l, m} (\hat{r}) \end{matrix}

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

\begin{aligned} Y_{1, - 1} & = i \sqrt{\frac{1}{2}} (Y_{1}^{- 1} + Y_{1}^{1}) = \sqrt{\frac{3}{4 π}} \cdot \frac{y}{r} \\ Y_{1, 0} & = Y_{1}^{0} = \sqrt{\frac{3}{4 π}} \cdot \frac{z}{r} \\ Y_{1, 1} & = \sqrt{\frac{1}{2}} (Y_{1}^{- 1} + Y_{1}^{1}) = \sqrt{\frac{3}{4 π}} \cdot \frac{x}{r} \end{aligned}

or equivalently

\begin{matrix} (38) & x = \sqrt{\frac{4 π}{3}} r \cdot Y_{1, 1}; y = \sqrt{\frac{4 π}{3}} r \cdot Y_{1, - 1}; z = \sqrt{\frac{4 π}{3}} r \cdot Y_{1, 0} \end{matrix}

Substituting Eq. $(38)$ and Eq. $(37)$ into Eq. $(35)$ , one have

\begin{matrix} (39) & T_{i j} = e [\int_{0}^{\infty} R_{l_{1}} (r) R_{l_{2}} (r) r^{3} d r] \cdot \sqrt{\frac{4 π}{3}} (\begin{matrix} G (l_{1}, l_{2}, 1, m_{1}, m_{2}, 1) \\ G (l_{1}, l_{2}, 1, m_{1}, m_{2}, - 1) \\ G (l_{1}, l_{2}, 1, m_{1}, m_{2}, 0) \end{matrix}) \end{matrix}

where $G (l_{1}, l_{2}, l_{3}, m_{1}, m_{2}, m_{3})$ is the real Gaunt coefficients defined as

\begin{matrix} (40) & G (l_{1}, l_{2}, l_{3}, m_{1}, m_{2}, m_{3}) = \int_{0}^{π} \sin θ d θ \int_{0}^{2 π} d ϕ Y_{l_{1}, m_{1}} (\hat{r}) Y_{l_{2}, m_{2}} (\hat{r}) Y_{l_{3}, m_{3}} (\hat{r}) \end{matrix}

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. $(38)$ . Eq. $(39)$ 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{aligned} (41) & ⟨ ψ_{i} | p | ψ_{j} ⟩ & = \int ⟨ ψ_{i} | k ⟩ \cdot ℏ k \cdot ⟨ k | ψ_{j} ⟩ d k \\ (42) & = ℏ (\begin{array}{c} \int ψ_{i} (k)^{⋆} k_{x} ψ_{j} (k) d k \\ \int ψ_{i} (k)^{⋆} k_{y} ψ_{j} (k) d k \\ \int ψ_{i} (k)^{⋆} k_{z} ψ_{j} (k) d k \end{array}) \end{aligned}

where the wavefunction in momentum space can be written as

\begin{matrix} (43) & ψ (k) = ⟨ k | ψ ⟩ = i^{l} G_{l} (k) Y_{l, m} (\hat{k}) \end{matrix}

with $G_{l} (k)$ being the spherical Bessel transform (SBT) of the radial function $R_{l} (r)$

\begin{matrix} (44) & G_{l} (k) = \sqrt{\frac{2}{π}} \int_{0}^{\infty} R_{l} (r) j_{l} (k r) r^{2} d r \end{matrix}

Combining Eq. $(43)$ , Eq. $(44)$ and Eq. $(33)$ , one have the TDM evaluated in momentum space

\begin{aligned} T_{i j} & = \frac{i^{l_{2} - l_{1} + 1} ℏ^{2} e}{m (ε_{j} - ε_{i})} [\int_{0}^{\infty} G_{l_{1}} (k) G_{l_{2}} (k) k^{3} d k] \\ (45) & \times \sqrt{\frac{4 π}{3}} (\begin{array}{c} G (l_{1}, l_{2}, 1, m_{1}, m_{2}, 1) \\ G (l_{1}, l_{2}, 1, m_{1}, m_{2}, - 1) \\ G (l_{1}, l_{2}, 1, m_{1}, m_{2}, 0) \end{array}) \end{aligned}

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

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 = 2 n$ , where $n \in N$ in Eq. $(45)$ .
The properties of the real Gaunt coefficients require that $l_{1} + l_{2} + 1 = 2 n$ , where $n \in N$ . ⁸
The spherical harmonics have definite parity, i.e. they are either even or odd with respect to inversion about the origin ⁹
$Y_{l, m} (\hat{r}) = (- 1)^{l} Y_{l, m} (- \hat{r})$
Moreover, $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{matrix} (46) & x = \sqrt{\frac{2 π}{3}} r \cdot [Y_{1}^{1} + Y_{1}^{- 1}]; y = \sqrt{\frac{2 π}{3}} r \cdot i [Y_{1}^{- 1} + Y_{1}^{1}]; z = \sqrt{\frac{4 π}{3}} r \cdot Y_{1}^{0} \end{matrix}$
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{matrix} (47) & Δ m = m_{1} - m_{2} = 0 or \pm 1 \end{matrix}$

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

#!/usr/bin/env python

import numpy as np
from sympy import lambdify
from sympy.abc import r
from sympy.physics.hydrogen import R_nl, E_nl

from pysbt import sbt, GauntTable

# All the quantities are in atomic units, i.e.
# Hartree for energy and Bohr for length.
# hbar = 1; m = 1
N    = 512
rmin = 2.0 / 1024 / 32
rmax = 30
# Logarithmic radial grid
rr   = np.logspace(np.log10(rmin), np.log10(rmax), N, endpoint=True)
# init spherical Bessel transform
ss   = sbt(rr)

n1 = 1; l1 = 0; m1 = 0 # 1s
n2 = 2; l2 = 1; m2 = 1 # 2px

R1 = lambdify((r), R_nl(n1, l1, r), 'numpy')(rr)
R2 = lambdify((r), R_nl(n2, l2, r), 'numpy')(rr)

G1 = ss.run(R1, l=l1, norm=True)
G2 = ss.run(R2, l=l2, norm=True)

# energy difference between the two states
dE = E_nl(n2) - E_nl(n1)

# Radial integration over r and k
Tr = np.sum(R1 * R2 * ss.rr**3 * ss.simp_wht_rr)
Tp = np.sum(G1 * G2 * ss.kk**3 * ss.simp_wht_kk)

print(f'R-space radial integration: Tr = {Tr:8.6f}')
print(f'K-space radial integration: Tp = {Tp:8.6f}')
print(f'Tp / dE = {Tp / dE:8.6f}')

Ta  = np.sqrt(4*np.pi / 3) * np.array([
    GauntTable(l1, l2, 1, m1, m2, m) for m in [1, -1, 0]
])
Tij =  Tr * Ta
print(f'Angular term: Ta = ({Ta[0]:8.6f}, {Ta[1]:8.6f}, {Ta[2]:8.6f})')
print(f'   Tij = Ta * Tr = ({Tij[0]:8.6f}, {Tij[1]:8.6f}, {Tij[2]:8.6f})')

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

R-space radial integration: Tr = 1.290266
K-space radial integration: Tp = 0.483850
Tp / dE = 1.290266
Angular term: Ta = (0.577350, 0.000000, 0.000000)
   Tij = Ta * Tr = (0.744936, 0.000000, 0.000000)

TDM between Bloch states

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

\begin{matrix} (48) & | ψ_{n k} ⟩ = e^{i k \cdot r} | u_{n k} ⟩ \end{matrix}

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

\begin{matrix} (49) & ⟨ ψ_{n k} | p | ψ_{m k} ⟩ = \frac{i m (ε_{n k} - ε_{m k})}{ℏ} ⟨ u_{n k} | i \nabla_{k} | u_{m k} ⟩ \end{matrix}

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

\begin{matrix} (50) & ⟨ ψ_{n k} | r | ψ_{m k} ⟩ \neq ⟨ u_{n k} | i \nabla_{k} | u_{m k} ⟩ \end{matrix}

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{aligned} (51) & ⟨ ψ_{c k} | p | ψ_{v k} ⟩ & = ℏ ⟨ u_{c, k} | k - i \nabla | u_{v, k} ⟩ \\ ⟨ u_{c, k} | k - i \nabla | u_{v, k} ⟩ & = ⟨ {\tilde{u}}_{c, k} | k - i \nabla | {\tilde{u}}_{v, k} ⟩ - i \sum_{n m} ⟨ {\tilde{u}}_{c, k} | {\tilde{p}}_{n, k} ⟩ ⟨ {\tilde{p}}_{m, k} | {\tilde{u}}_{v, k} ⟩ \\ (52) & \times [⟨ ϕ_{n} | \nabla_{r} | ϕ_{m} ⟩ - ⟨ {\tilde{ϕ}}_{n} | \nabla_{r} | {\tilde{ϕ}}_{m} ⟩] \end{aligned}

where ${\tilde{p}}_{n k}$ is the k-dependent projector function as defined in a previous post, $ϕ$ and $\tilde{ϕ}$ 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{aligned} ϵ_{α β}^{″} (ω) & = \frac{4 π^{2} e^{2}}{Ω} lim_{q \to 0} \frac{1}{q^{2}} \sum_{c, v, k} 2 w_{k} δ (ε_{c, k} - ε_{v, k} - ω) \\ (53) & \times ⟨ u_{c, k + q e_{α}} | u_{v, k} ⟩ ⟨ u_{c, k} | u_{v, k + q e_{β}} ⟩^{*} \end{aligned}

and in the transversal expression by ¹²

\begin{aligned} ϵ_{α β}^{″} (ω) & = \frac{4 π^{2} e^{2} ℏ^{4}}{Ω ω^{2} m_{e}^{2}} lim_{q \to 0} \sum_{c, v, k} 2 w_{k} δ (ε_{c, k + q} - ε_{v, k} - ω) \\ (54) & \times ⟨ u_{c, k} | i \nabla_{α} - k_{α} | u_{v, k} ⟩ ⟨ u_{c, k} | i \nabla_{β} - k_{β} | u_{v, k} ⟩^{*} \end{aligned}

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{matrix} (55) & ⟨ u_{c, k} | k - i \nabla | u_{v, k} ⟩ = \sum_{G} c_{c, k}^{*} (G) \cdot [k + G] \cdot c_{v, k} (G) \end{matrix}

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

Matrix element of nabla operator — one-center contribution

\begin{matrix} (56) & τ_{n m} = ⟨ ϕ_{n} | \nabla_{r} | ϕ_{m} ⟩ - ⟨ {\tilde{ϕ}}_{n} | \nabla_{r} | {\tilde{ϕ}}_{m} ⟩ \end{matrix}

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

\begin{aligned} (57) & ϕ_{n} (r) & = U_{n, l} (r) \cdot Y_{l, m} (\hat{r}) \\ (58) & {\tilde{ϕ}}_{n} (r) & = {\tilde{U}}_{n, l} (r) \cdot Y_{l, m} (\hat{r}) \end{aligned}

and their Fourier transform

\begin{aligned} (59) & ϕ_{n} (k) & = i^{l} G_{n, l} (k) \cdot Y_{l, m} (\hat{k}) \\ (60) & {\tilde{ϕ}}_{n} (k) & = i^{l} {\tilde{G}}_{n, l} (k) \cdot Y_{l, m} (\hat{k}) \end{aligned}

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{aligned} ⟨ ϕ_{n} | \nabla_{r} | ϕ_{m} ⟩ & = \int ϕ_{n} (r) \nabla_{r} ϕ_{m} (r) d r \\ = \int U_{n, l^{'}} (r) Y_{l^{'}, m} (\hat{r}) \nabla_{r} [\frac{U_{n, l} (r)}{r^{l}} \cdot r^{l} Y_{l, m} (\hat{r})] r^{2} d r d Ω \\ = \int U_{n, l^{'}} (r) r^{l + 2} \frac{\partial}{\partial r} [\frac{U_{n, l} (r)}{r^{l}}] Y_{l^{'}, m} (\hat{r}) Y_{l, m} (\hat{r}) {\vec{e}}_{r} d r d Ω \\ (61) & + \int U_{n, l^{'}} (r) U_{n, l} (r) r^{2 - l} Y_{l^{'}, m} (\hat{r}) \nabla_{r} [r^{l} Y_{l, m} (\hat{r})] d r d Ω \end{aligned}

For the first term in the last line of Eq. $(61)$ , remember that
$\begin{matrix} (62) & {\vec{e}}_{r} = (\begin{matrix} \frac{x}{r}, & \frac{y}{r}, & \frac{z}{r} \end{matrix}) \end{matrix}$
which can be written as real spherical harmonics with $l = 1$ , i.e.
$\begin{matrix} (63) & {\vec{e}}_{r} = \sqrt{\frac{4 π}{3}} (\begin{matrix} Y_{1, 1} (\hat{r}), & Y_{1, - 1} (\hat{r}), & Y_{1, 0} (\hat{r}) \end{matrix}) \end{matrix}$
As a result, the first term reduces to
$\begin{matrix} (64) & \sqrt{\frac{4 π}{3}} \int_{0}^{\infty} U_{n, l^{'}} (r) r^{l + 2} \frac{\partial}{\partial r} [\frac{U_{n, l} (r)}{r^{l}}] d r \times (\begin{matrix} G (l^{'}, l, 1, m^{'}, m, 1) \\ G (l^{'}, l, 1, m^{'}, m, - 1) \\ G (l^{'}, l, 1, m^{'}, m, 0) \end{matrix}) \end{matrix}$
For the second term in Eq. $(61)$ , one can show that ⁶
$\begin{matrix} (65) & \nabla_{r} [r^{l} Y_{l, m} (\hat{r})] = r^{l - 1} \sum_{m = 1 - l}^{l - 1} c_{m} Y_{l - 1, m} (\hat{r}) \end{matrix}$
Therefore, the second term reduces to
$\begin{matrix} (66) & \int_{0}^{\infty} U_{n, l^{'}} (r) U_{n, l} (r) r d r \times \int Y_{l^{'}, m} (\hat{r}) r^{1 - l} \nabla_{r} [r^{l} Y_{l, m} (\hat{r})] d Ω \end{matrix}$
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. $(64)$ and Eq. $(66)$ , one have

\begin{aligned} ⟨ ϕ_{n} | \nabla_{r} | ϕ_{m} ⟩ & = \sqrt{\frac{4 π}{3}} \int_{0}^{\infty} [R_{n, l^{'}} \frac{\partial R_{n, l}}{\partial r} - (l + 1) \frac{R_{n, l^{'}} R_{n, l}}{r}] d r \times (\begin{array}{c} G (l^{'}, l, 1, m^{'}, m, 1) \\ G (l^{'}, l, 1, m^{'}, m, - 1) \\ G (l^{'}, l, 1, m^{'}, m, 0) \end{array}) \\ (67) & + \int_{0}^{\infty} \frac{R_{n, l^{'}} R_{n, l}}{r} d r \times \int Y_{l^{'}, m} (\hat{r}) r^{1 - l} \nabla_{r} [r^{l} Y_{l, m} (\hat{r})] d Ω \end{aligned}

All the terms with spherical harmonics in Eq. $(67)$ can be obtained analytically.

The subroutine NABLA_RADIAL in VASP source file optics.F calculates $τ_{n m}$ 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_{\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 $τ_{n m}$ to file “nabla_X”, where X is the name of the elements.

 IF (IO%LOPTICS) CALL SET_NABIJ_AUG(P,T_INFO%NTYP)
 ! Add the code after the previous line
 io_begin
 do nt = 1, T_INFO%NTYP
   open(unit=124, file='nabla_' // trim(P(NT)%ELEMENT), status='unknown', action='write')
   do i = 1, 3
     do j = 1, P(NT)%LMMAX
       write(124, fmt=*) (P(NT)%NABLA(i,j,n), n=1, P(NT)%LMMAX)
     end do
     write(124,*) 
   end do
   close(unit=124)
 end do
 io_end

Momentum-space integration

Similar to the evaluaton of Eq. $(45)$

\begin{aligned} ⟨ ϕ_{n} | \nabla_{r} | ϕ_{m} ⟩ & = i^{l_{2} - l_{1} + 1} \int_{0}^{\infty} G_{n, l_{1}} (k) G_{m, l_{2}} (k) k^{3} d k \\ (68) & \times \sqrt{\frac{4 π}{3}} (\begin{array}{c} G (l_{1}, l_{2}, 1, m_{1}, m_{2}, 1) \\ G (l_{1}, l_{2}, 1, m_{1}, m_{2}, - 1) \\ G (l_{1}, l_{2}, 1, m_{1}, m_{2}, 0) \end{array}) \end{aligned}

Therefore

\begin{aligned} τ_{n m} & = ⟨ ϕ_{n} | \nabla_{r} | ϕ_{m} ⟩ - ⟨ {\tilde{ϕ}}_{n} | \nabla_{r} | {\tilde{ϕ}}_{m} ⟩ \\ = i^{l_{2} - l_{1} + 1} \int_{0}^{\infty} [G_{n, l_{1}} (k) G_{m, l_{2}} (k) - {\tilde{G}}_{n, l_{1}} (k) {\tilde{G}}_{m, l_{2}} (k)] k^{3} d k \\ (69) & \times \sqrt{\frac{4 π}{3}} (\begin{array}{c} G (l_{1}, l_{2}, 1, m_{1}, m_{2}, 1) \\ G (l_{1}, l_{2}, 1, m_{1}, m_{2}, - 1) \\ G (l_{1}, l_{2}, 1, m_{1}, m_{2}, 0) \end{array}) \end{aligned}

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

#!/usr/bin/env python

from paw import pawpotcar

pp    = pawpotcar(potfile='pot.C')
# evaluate tau_nm in real-space, False for momentum-space evaluation
nij   = pp.get_nablaij(lreal=True)
xyz   = 'xyz'

for ii in range(3):
    head = ' '.join(
        [f'({n}, {l},{m:2d})'.rjust(10) for n, l, m in pp.ilm]
    )
    print(f'tau_nm({xyz[ii]})'.center(10) + head)
    for jj in range(pp.lmmax):
        n, l, m = pp.ilm[jj]
        c1 = f'({n}, {l},{m:2d})'.center(10)
        line = ' '.join([
            f'{x:10.6f}' for x in nij[ii,jj,:]
        ])
        print(c1 + line)
    print()

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.

tau_nm(x)  (0, 0, 0)  (1, 0, 0)  (2, 1,-1)  (2, 1, 0)  (2, 1, 1)  (3, 1,-1)  (3, 1, 0)  (3, 1, 1)
(0, 0, 0)   0.000000   0.000000  -0.000000  -0.000000  -0.417279  -0.000000  -0.000000   0.941495
(1, 0, 0)   0.000000   0.000000  -0.000000  -0.000000  -0.557957  -0.000000  -0.000000   1.261615
(2, 1,-1)   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(2, 1, 0)   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(2, 1, 1)   0.417279   0.557957   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(3, 1,-1)   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(3, 1, 0)   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(3, 1, 1)  -0.941495  -1.261615   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000

tau_nm(y)  (0, 0, 0)  (1, 0, 0)  (2, 1,-1)  (2, 1, 0)  (2, 1, 1)  (3, 1,-1)  (3, 1, 0)  (3, 1, 1)
(0, 0, 0)   0.000000   0.000000  -0.417279  -0.000000  -0.000000   0.941495  -0.000000  -0.000000
(1, 0, 0)   0.000000   0.000000  -0.557957  -0.000000  -0.000000   1.261615  -0.000000  -0.000000
(2, 1,-1)   0.417279   0.557957   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(2, 1, 0)   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(2, 1, 1)   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(3, 1,-1)  -0.941495  -1.261615   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(3, 1, 0)   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(3, 1, 1)   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000

tau_nm(z)  (0, 0, 0)  (1, 0, 0)  (2, 1,-1)  (2, 1, 0)  (2, 1, 1)  (3, 1,-1)  (3, 1, 0)  (3, 1, 1)
(0, 0, 0)   0.000000   0.000000  -0.000000  -0.417279  -0.000000  -0.000000   0.941495  -0.000000
(1, 0, 0)   0.000000   0.000000  -0.000000  -0.557957  -0.000000  -0.000000   1.261615  -0.000000
(2, 1,-1)   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(2, 1, 0)   0.417279   0.557957   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(2, 1, 1)   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(3, 1,-1)   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(3, 1, 0)  -0.941495  -1.261615   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000
(3, 1, 1)   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000   0.000000

Light-Matter Interaction and Dipole Transition Matrix

Quantum theory of light-matter interaction

Electric Dipole Approximation

Coulomb Gauge

Length Gauge

Velocity Gauge

Dipole and Momentum Transition Matrix

Transition between Molecular Orbitals

TDM from p-r relation

Selection Rules

TDM between hydrogen orbitals

TDM between Bloch states

Matrix element of nabla operator — pseudo-wavefunction

Matrix element of nabla operator — one-center contribution

Real-space integration

Momentum-space integration

Reference

Trending Tags

Light-Matter Interaction and Dipole Transition Matrix

Quantum theory of light-matter interaction

Electric Dipole Approximation

Coulomb Gauge

Length Gauge

Velocity Gauge

Dipole and Momentum Transition Matrix

Transition between Molecular Orbitals

TDM from p-r relation

Selection Rules

TDM between hydrogen orbitals

TDM between Bloch states

Matrix element of nabla operator — pseudo-wavefunction

Matrix element of nabla operator — one-center contribution

Real-space integration

Momentum-space integration

Reference

Further Reading

PAW All-Electron Wavefunction in VASP

Band Unfolding Tutorial

VASP: Stacking-dependent Potential Energy Surface of Bilayer MoS2

Trending Tags

VASP: Stacking-dependent Potential Energy Surface of Bilayer MoS₂