Responses to external fields

Compared to atomic fields, the externally applied electromagnetic fields are several orders of magnitude weaker and perturbation theory is called for and a main merit of such an approach is that the response functions defined in perturbation theory provide the natural meeting point between experiment and theory with a distinct separation of one-, two-, three-photon, etc., optical processes

Response theory Norman et al. (2018) may be thought of as a reformulation of standard time-dependent perturbation theory into a form suitable for approximate state theory. Virtually all spectroscopic properties are encompassed by the theory as possible perturbations include

time-independent or time-dependent
electric or magnetic
internal or external
geometric distortions

The derivation of expressions for response functions can appear very different from one source to another and the vast number of technicalities can at first appear overwhelming. The principle, however, is quite straightforward and clear: Form a well-defined quantity of interest, e.g., the electric dipole moment, $\mu(t)$ , and identify response functions in the order expansions of these quantities.

Defining the system¶

In the following, we assume that the system is described by a Hamiltonian of the form

\hat{H} = \hat{H}_0 + \hat{V}(t)

(1)

where $\hat{H}_0$ is the Hamiltonian of the isolated system and $\hat{V}(t)$ described the quantum mechanical field-coupling. In the electric-dipole approximation, the coupling between the molecular system and external electric field is provided by minus the electric dipole moment operator, $-\hat{\mu}$ , such that

\hat{V}(t) = - \hat{\mu} F(t)

(2)

where $F(t)$ is the electric field amplitude. For simplicity, we assume the external electric field to be monochromatic and adiabatically switched on at time zero according to

F(t) = F^\omega \sin \omega t \times \mathrm{erf}(a t)

(3)

We further assume the solutions to the eigenvalue problem of $\hat{H}_0$ to be known

\hat{H}_0 | n \rangle = E_n | n \rangle

(4)

such that it can be used as a basis in the calculation. Before being exposed to the perturbation, we assume the molecule to be in the molecular ground state, $| 0 \rangle$ .

Time propagation¶

With a time step $\Delta t$ that is small compared to the time variations in the Hamiltonian, we can approximate the Hamiltonian to be constant from time $t_0$ to $t_0 + \Delta t$ . It is then possible to perform a time integration of the Schröinger equation using

\Psi(t_0 + \Delta t) = \hat{U}(t_0+\Delta t, t_0) \Psi(t_0)

where the time-evolution propagator is defined as

\hat{U}(t_0+\Delta t, t_0) = e^{-i \hat{H}(t_0) \Delta t / \hbar}

Repeated applications of $\hat{U}$ allow us to determine the wave function in the region $t>0$ and, given the time-dependent wave function, the dipole moment is obtained as the expectation value of the electric dipole operator

\mu(t) = \langle \Psi(t) | \hat{\mu} | \Psi(t) \rangle

Response functions¶

Response functions are defined as the expansion coefficients of Fourier expansions of time-dependent expectation values of an observable associated with an operator $\hat{\Omega}$

\begin{align*} \langle \Psi(t) | \hat{\Omega} | \Psi(t) \rangle & = & \langle 0 | \hat{\Omega} | 0 \rangle + \sum_{\omega} \langle \langle \hat{\Omega}; \hat{V} \rangle \rangle_{\omega} F^{\omega} e^{-i\omega t} + \cdots \end{align*}

(5)

In the case of the observable being the electric dipole moment, the linear coupling coefficient to the external electric field is known as the electric-dipole polarizability

\mu(t) = \mu^0 + \sum_\omega \alpha(\omega) F^\omega e^{-i\omega t} + \cdots

(6)

such that

\alpha(\omega) = - \langle \langle \hat{\mu}; \hat{\mu} \rangle \rangle_{\omega}

(7)

Two-level system¶

import matplotlib.pyplot as plt
import numpy as np
from scipy.linalg import expm
from scipy.special import erf

The representation of an atomic system with two levels takes the form

| \Psi(t) \rangle = a(t) | \Psi_a \rangle + b(t) | \Psi_b \rangle = \begin{pmatrix} a(t) \\ b(t) \end{pmatrix}

(8)

\hat{H}_0 = \begin{pmatrix} E_a & 0 \\ 0 & E_b \\ \end{pmatrix}; \qquad \hat{\mu} = \begin{pmatrix} 0 & \mu_{ab} \\ \mu_{ba} & 0 \\ \end{pmatrix}

(9)

# Time
delta = 0.1
t = np.arange(0, 500, delta)
n = len(t)

# External electric field
Fw = 0.00005
w = 0.1
F = Fw * np.sin(w * t) * erf(0.01 * t)

# Time-dependent wave function
psi = np.zeros((2, n), dtype=complex)
psi[:, 0] = [1, 0]

# Energy of levels
Ea = 1
Eb = 1.5

# Transition dipole moment
muab = 1
muba = muab

# Hamiltonian and dipole moment operators
H0 = np.array([[Ea, 0], [0, Eb]])
mu = np.array([[0, muab], [muba, 0]])
V = -mu

# Time-dependent dipole moment (polarization)
P = np.zeros(n, dtype=float)
P[0] = np.dot(psi[:, 0], np.matmul(mu, psi[:, 0])).real

Time propagation¶

for k in range(1, n):

    psi[:, k] = np.matmul(expm(-1j * (H0 + V * F[k - 1]) * delta), psi[:, k - 1])

    P[k] = np.dot(np.conjugate(psi[:, k]), np.matmul(mu, psi[:, k])).real

Plotting observables¶

fig = plt.figure(figsize=(10, 6))

ax = plt.subplot(2, 1, 1, xlim=(0, 500))
plt.plot(t, P, "k-", linewidth=2.0)
plt.grid(True)
plt.ylabel(r"$\mu(t)$")

ax = plt.subplot(2, 1, 2, xlim=(0, 500))
plt.plot(t, F, "k-", linewidth=2.0)
plt.grid(True)
plt.ylabel(r"$F(t)$")

plt.xlabel("Time (a.u.)")

plt.show()

Determining response functions¶

As nonlinear responses are negligible, the value of the linear response function (or polarizability) at the given optical frequency is simply the ratio of the induced dipole moment and the amplitude of the electric field

\alpha(\omega) = \frac{\max[ \mu(t) ]}{F^\omega}

alpha_numerical = np.max(P) / Fw

print(f"Polarizability: {alpha_numerical : 4.2f} a.u.")

Analytic response functions¶

In practice, we derive analytic expressions for the response functions with use of perturbation theory. This results in sum-over-state formulas such as that for the linear response function

\langle \langle \hat{\Omega}; \hat{V} \rangle \rangle_\omega = - \frac{1}{\hbar} \sum_{n>0} \left[ \frac{\langle 0 | \hat{\Omega} | n \rangle \langle n | \hat{V} | 0 \rangle } {\omega_{n0}-\omega} + \frac{\langle 0 | \hat{V} | n \rangle \langle n | \hat{\Omega} | 0 \rangle } {\omega_{n0}+\omega} \right]

Eba = Eb - Ea

alpha_analytic = muab * muba * (1.0 / (Eba - w) + 1.0 / (Eba + w))

print(f"Polarizability: {alpha_analytic : 4.2f} a.u.")

Poles of response functions¶

Response functions have poles at transition energies and the associated residues are products of transition moments.

Most excited state methods in computational chemistry are based on the principle of identifying excitation energies and spectrum intensities from a residue analysis of response functions.

Dispersion¶

The variations of response functions with optical frequencies are known as dispersions and values for non-oscillating fields ( $\omega = 0$ ) are known as static values.

epsilon = 0.01
wa = np.arange(0, Eba - epsilon, 0.001)
wb = np.arange(Eba + epsilon, 1, 0.001)

alpha_wa = muab * muba * (1.0 / (Eba - wa) + 1.0 / (Eba + wa))
alpha_wb = muab * muba * (1.0 / (Eba - wb) + 1.0 / (Eba + wb))

ax = plt.subplot(1, 1, 1, xlim=(0, 1), ylim=(-40, 40))

plt.plot(wa, alpha_wa, lw=2)
plt.plot(wb, alpha_wb, lw=2)

plt.vlines(Eba, -40, 40, color="grey", linestyles="--")

plt.ylabel(r"$\alpha(\omega)$ (a.u.)")
plt.xlabel(r"$\omega$ (a.u.)")

plt.annotate(
    "Static value",
    (0, alpha_wa[0]),
    xytext=(0.05, 15),
    arrowprops={"width": 1, "headwidth": 6, "headlength": 6, "color": 'black'},
)

plt.annotate(
    "Pole",
    (Eba, 0),
    xytext=(0.55, 10),
    arrowprops={"width": 1, "headwidth": 6, "headlength": 6, 'color': 'black'},
)

plt.grid(True)

plt.show()