# Relativistic Real-Time Methods

Marius Kadek  Lukas Konecny  Michal Repisky  [<EMAIL_ADDRESS>
###### Abstract

Recent advances in laser technology enable to follow electronic motion at its
natural time-scale with ultrafast pulses, leading the way towards atto- and
femtosecond spectroscopic experiments of unprecedented resolution.
Understanding of these laser-driven processes, which almost inevitably involve
non-linear light–matter interactions and non-equilibrium electron dynamics, is
challenging and requires a common effort of theory and experiment. Real-time
electronic structure methods provide the most straightforward way to simulate
experiments and to gain insights into non-equilibrium electronic processes. In
this Chapter, we summarize the fundamental theory underlying the relativistic
particle–field interaction Hamiltonian as well as equation-of-motion for
exact-state wave function in terms of the one- and two-electron reduced
density matrix. Further, we discuss the relativistic real-time electron
dynamics mean-field methods with an emphasis on Density-Functional Theory and
Gaussian basis, starting from the four-component (Dirac) picture and continue
to the two-component (Pauli) picture, where we introduce various flavours of
modern exact two-component (X2C) Hamiltonians for real-time electron dynamics.
We also overview several numerical techniques for real-time propagation and
signal processing in quantum electron dynamics. We close this Chapter by
listing selected applications of real-time electron dynamics to frequency-
resolved and time-resolved spectroscopies.

###### keywords:

real-time, electron dynamics, Liouville-von Neumann equation, reduced density
matrix, density functional theory, noncollinearity, relativistic theory, Dirac
Hamiltonian, X2C Hamiltonian, Fourier transformation, absorption, circular
dichroism, nonlinear spectroscopy, pump-probe spectroscopy, time-resolved
spectroscopy

††journal: Comprehensive Computational Chemistry

url]www.respectprogram.org/michalrepisky

[inst1]organization=Hylleraas Centre for Quantum Molecular Sciences,
Department of Chemistry, UiT The Arctic University of Norway, city=Tromsø,
postcode=9037, country=Norway

[inst2]organization=Department of Physics, College of Science, Northeastern
University, city=Boston, postcode=Massachusetts 02115, country=USA

[inst3]organization=Max Planck Institute for the Structure and Dynamics of
Matter, Center for Free Electron Laser Science, city=Hamburg, postcode=22761,
country=Germany

[inst4]organization=Department of Physical and Theoretical Chemistry, Faculty
of Natural Sciences, Comenius University, city=Bratislava, postcode=84104,
country=Slovakia

##  1 Objectives box

The principal objectives of this Chapter are:

  * 1.

Introduction to the fundamental theory leading to the relativistic
particle–field interaction Hamiltonian.

  * 2.

Discussion of the equations-of-motion for exact-state wave function in terms
of the one-electron and two-electron reduced density matrix.

  * 3.

Introduction to the relativistic four-component real-time electron dynamics
mean-field methods with an emphasis on Density-Functional Theory and Gaussian
basis.

  * 4.

Detailed overview of various exact two-component (X2C) transformations towards
the relativistic two-component real-time electron dynamics.

  * 5.

Overview of numerical techniques for real-time propagation and signal
processing in quantum electron dynamics.

  * 6.

Selected application of real-time electron dynamics to frequency-resolved and
time-resolved spectroscopies.

##  2 Introduction

The rapid advancement of laser technology in the past decades allows us to
probe matter on spatiotemporal scales that approach the characteristic time
and length scales of the electron, opening the field of attosecond science [1,
2]. This development has forced quantum chemists to shift their attention from
the time-independent to the time-dependent Schrödinger or Dirac equation.
Real-time electronic structure theory thus describes the explicit time-
evolution of the wave function or the electron density driven by non-
equilibrium condition of the Hamiltonian under external perturbation(s). All
physical quantities of a molecular system are then extracted non-
perturbatively from the time-varying part of the wave function or the electron
density. Due to non-perturbative nature, real-time methods represent the most
straightforward approach to dynamical property calculations and enable the use
of external perturbations of arbitrary strength, shape and duration, capturing
in general both linear and non-linear effects within a wide spectral window
from a single run [2, 3, 4]. This distinguishes real-time electronic structure
theory from response theory where all physical quantities are obtained in the
frequency domain using perturbation expansion [5].

Historically, an early work on explicit time-propagation of electronic wave
function dates back to 1990 when Cederbaum and coworkers developed the
nonrelativistic multiconfigurational time-dependent Hartree (MCTDH) method
[6]. Shortly after, Micha and Runge developed a real-time time-dependent
Hartree–Fock (RT-TDHF) approach that couples electronic and nuclear motions
[7], whereas Theilhaber [8], and Yabana and Bertsch [9], introduced the first-
ever real-time time-dependent density functional theory (RT-TDDFT) combining
the local density approximation with real-space grid methodology. In the
condensed matter physics community, these pioneering works led to several
implementations of the time-propagation formalism using either localized basis
sets [10, 11], plane waves [12, 13, 14], or real-space grids [15].
Advancements in computing power and numerical algorithms have enabled
performing large-scale RT-TDDFT simulations even on periodic solids [16, 17,
18, 19, 20]. In the quantum chemistry community, the first nonrelativistic RT-
TDDFT implementation based on popular Gaussian-type atomic orbitals was
pioneered by Isborn and coworkers in 2007 [21] and later adopted by several
other groups [22, 23]. The extension of RT-TDDFT to the relativistic four-
component (4c) realm was presented by Repisky and coworkers in 2015 [24, 25].
As shown by these authors, significant gains in computer time was obtained by
transforming the parent 4c RT-TDDFT to an exact two-component (X2C) form [26],
although the accuracy of reference 4c results was achieved only after
inclusion of the two-electron and exchange–correlation picture-change
corrections [27, 28]. Beyond DFT, there has been growing interest in explicit
time propagation of correlated methods such as multiconfigurational self-
consistent-field [29, 30, 31, 32], configuration interaction [33, 34, 35, 36,
37, 38, 39], algebraic diagrammatic construction [40, 41, 42, 43, 44, 45, 46],
density matrix renormalization group[47], Møller-Plesset [48], and coupled
cluster [4, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61] theories.

At the non-relativistic level of theory, a plethora of applications of real-
time methods has been presented including UV/Vis absorption spectroscopy [9,
62, 22], excited-state absorption [63], photoionization [64, 65], X-ray
absorption [66, 67], chiroptical spectroscopies [68, 69, 70, 71, 72], non-
linear optical properties [73, 74], spin and magnetization dynamics [75, 76],
molecular conductance [77], pump-probe spectroscopy [78], photoinduced
electric currents [79], plasmon resonances [23], singlet–triplet transitions
[80], and magnetic circular dichroism [81]. This list of applications is by no
means exhaustive and a reader interested in a more thorough exploration of the
use of non-relativistic real-time methods is referred to the recent review [3]
and references therein.

The advent of soft X-ray free electron laser pulses with subfemtosecond
temporal widths have opened new ways to investigate time-resolved dynamics
involving inner-shell electrons. A prerequisite for reliable quantum-chemical
modeling of these processes is the inclusion of relativistic effects, defined
as differences between the exact Dirac (four-component) description of matter
and an approximate Schrödinger (one-component) description. This requirement
stems from the fact that the inner-shell orbitals involved in X-ray
absorption/emission processes are most affected by relativity, manifestations
of which are frequency shifts of spectral lines due to the scalar (SC)
relativistic effects as well as spectral fine structure splitting arising from
the spin-orbit (SO) coupling [25, 28]. The relativistic effects are
significant even in light (third row) elements [25] and increase in importance
for heavier elements [28], highlighting the need for relativistic description
across the Periodic Table. Therefore, the most accurate way to perform real-
time simulations is the use of full four-component (4c) Dirac formalism where
both the SC and SO relativistic effects are included variationally. The first
4c extension of the time-propagation formalism was presented by Repisky and
coworkers at the RT-TDDFT level [24, 25], involving the program package
ReSpect [82]. Recently, De Santis and coworkers reported a similar 4c RT-TDDFT
implementation in the BERTHA code [83]. While advancements in computing power
and numerical algorithms have enabled performing fairly large 4c real-time
electron dynamics simulations [24, 25, 26, 27, 84, 83], there is still
interest in developing approximate two-component (2c) methods that maintain
the accuracy of the parent 4c method at a fraction of its computational cost.
In this respect, the X2C Hamiltonian has gained wide popularity in quantum
chemistry community as it reduces the original 4c problem by half at the
expense of only a few simple algebraic manipulations [85, 86, 87, 88, 89]. As
shown independently by Konecny [26] and Goings [90], the central idea of X2C
transformation can be extended to the real-time electron dynamics framework,
provided the X2C decoupling matrix satisfies an adiabatic approximation [26].
However, both real-time X2C implementations utilize a crude one-electron X2C
(1eX2C) Hamiltonian model, which typically leads to absolute errors for core
spinor energies of heavier elements of the order of tens of Hartree [91]. As
shown by Knecht and coworkers [91], accuracy of X2C Hamiltonians severely
depends on the two-electron and exchange-correlation picture-change correction
models employed, and can vary as much as 5-6 orders of magnitude for core-
shell energies. As a remedy, the authors introduced two simple yet
computationally efficient and numerically accurate X2C Hamiltonian models,
dubbed as amfX2C and e(xtended)amfX2C, to correct both SC and SO two-electron
and exchange-correlation picture-change effects using simple atomic mean-field
quantities, achieving a consistent $\approx 10^{-5}$ Hartree/atom accuracy
[91]. The theoretical extension and numerical assessment of (e)amfX2C
Hamiltonian models was recently performed for conventional and time-resolved
TDDFT by Repisky and coworkers [27, 28]. In addition to the previous works,
this Chapter also provides an in-depth discussion on the transformation of the
original 4c equation-of-motion to its 2c form, particularly focusing on the
modern exact two-component (X2C) formalism.

At the relativistic level, the real-time applications are scarcer due to fewer
computer programs providing such functionality. These programs include ReSpect
[24, 82], Gaussian [90], Chronus Quantum [92], PyBerthaRT [83], BDF [93], and
FHI-aims [11]. The molecular properties addressed at the relativistic level
include UV/Vis absorption spectroscopy [24, 90], X-ray absorption [25, 94,
93], non-linear optical properties [26], chiroptical spectroscopies [84], high
harmonic generation [83], and pump-probe spectroscopy [27].

Before closing this Section, let us emphasize that we restrict ourselves to
(i) the Born–Oppenheimer approximation and therefore coupled electron–nuclear
dynamics is not considered here; (ii) the semiclassical approximation where
electronic degrees of freedom are described quantum mechanically, while
electromagnetic fields are treated classically. In next Sections, we introduce
the fundamental theory behind the relativistic particle–field interaction
Hamiltonian, and discuss the equation-of-motion for exact-state wave function
in terms of the one-electron and two-electron reduced density matrix. Later,
we dive into the relativistic four-component real-time electron dynamics mean-
field methods with an emphasis on Density-Functional Theory and Gaussian
basis, followed by a detailed overview of various exact two-component (X2C)
transformation models within the time domain. Finally, we offer a brief
overview of numerical techniques for real-time propagation and signal
processing, and close this Chapter by listing selected applications in
relativistic quantum electron dynamics.

##  3 Relativistic Particle–Field Interaction Hamiltonians

For the theoretical description of spectroscopic processes, quantum chemistry
commonly employs a semiclassical theory. In this framework, the molecules are
described quantum-mechanically, whereas the electromagnetic (photon) field is
treated classically. This assumption is justified in the limit of large photon
numbers – to be specific, when the photon density exceeds one per cubic
wavelength (for a discussion, see Ref. [95]). If this is not the case, one may
need to quantize the photon field as well and work within the framework of
quantum electrodynamics [96]. To provide an illustrative example, let us
consider a laser pulse with intensity $10^{14}\mathrm{W/cm^{2}}$ and
wavelength $\lambda=1064\mathrm{nm}$. The number of photons per cubic
wavelength is then given by ($\hbar$ and $c$ denote the reduced Planck
constant and the speed of light, respectively):

$\frac{\text{energy
flux}}{\hbar\omega}\frac{V}{c}=\frac{10^{14}\mathrm{W/cm^{2}}}{1.86\times
10^{-19}\mathrm{J}}\frac{1.064^{3}\times 10^{-12}\mathrm{cm^{3}}}{3\times
10^{10}\mathrm{cm/s}}\approx 2\times 10^{10},$ (1)

which is obviously much greater than one. Therefore, this semiclassical
theoretical framework is appropriate for absorption and emission processes,
and we rely on this framework throughout this Chapter.

Before we actually dive into particle–field interactions, let us first
consider a $N$-electron system alone, _i.e._ in the absence of any
electromagnetic (photon) field. In this case, the system is governed by the
relativistic electronic Hamiltonian

$\hat{H}=\sum_{i}^{N}\hat{h}^{\mathrm{D}}_{i}+\frac{1}{2}\sum_{i\neq
j}^{N}\hat{g}_{ij}.$ (2)

Here, $\hat{h}^{\mathrm{D}}_{i}$ is the famous relativistic Dirac Hamiltonian
of a single electron $i$, while $\hat{g}_{ij}$ is the interaction Hamiltonian
between electrons $i$ and $j$. A factor one half in front of $\hat{g}$
corrects for double counting of the two-electron interactions.
$\hat{h}^{\mathrm{D}}$ describes the relativistic kinetic energy of an
electron as well as its interaction energy with the electrostatic scalar
potential $\phi_{0}(\bm{r})$ due to the fixed atomic nuclei. It bears the
$4\\!\times\\!4$ matrix form [97, 98]

$\begin{split}\hat{h}^{\mathrm{D}}_{i}&=\beta^{\prime}_{i}m_{e}c^{2}+c\big{(}\bm{\alpha}_{i}\cdot\bm{p}_{i}\big{)}-e\phi_{0}(\bm{r}_{i})\mathbb{I}_{4}.\end{split}$
(3)

Here, $\bm{r}_{i}$ and $\bm{p}_{i}\\!=\\!-i\hbar\bm{\nabla}_{i}$ refer to the
position and canonical momentum of the $i$th electron, respectively.
$\mathbb{I}_{4}$ is a $4\\!\times\\!4$ identity matrix, and $-e$, $m_{e}$ and
$c$ are constants referring to the electron charge, electron mass and the
speed of light in vacuum. When compared to the original expression of Dirac
[97, 98], $\hat{h}^{\mathrm{D}}$ utilizes the reduced rest mass energy
$\beta^{\prime}m_{e}c^{2}$ with $\beta^{\prime}\equiv\beta-\mathbb{I}_{4}$ to
align the relativistic and non-relativistic energy scales. $\beta$ is one of
four new $4\\!\times\\!4$ matrix variables

$\displaystyle\beta=\begin{pmatrix}\mathbb{I}_{2}&0_{2}\\\
0_{2}&-\mathbb{I}_{2}\end{pmatrix};\qquad\bm{\alpha}=\begin{pmatrix}0_{2}&\bm{\sigma}\\\
\bm{\sigma}&0_{2}\end{pmatrix},$ (4)

introduced by Dirac to formulate relativistic quantum-mechanical equations of
motion for spin-1/2 particles that are linear in space and time [97, 98].
These variables fulfill the anti-commutation relations

$[\alpha_{k},\beta]_{+}=0_{4};\qquad[\alpha_{k},\alpha_{l}]_{+}=2\delta_{kl}\mathbb{I}_{4};\qquad
k,l\in x,y,z,$ (5)

and are customarily written in terms of the two-component Pauli spin matrices

$\displaystyle\sigma_{x}=\begin{pmatrix}0&1\\\
1&0\end{pmatrix};\quad\sigma_{y}=\begin{pmatrix}0&-\mathrm{i}\\\
\mathrm{i}&0\end{pmatrix};\quad\sigma_{z}=\begin{pmatrix}1&0\\\
0&-1\end{pmatrix}.$ (6)

For further reading on the properties and physical interpretation of the Dirac
Hamiltonian, the reader is referred to several excellent quantum chemistry
textbooks [99, 100, 101].

Now, let us subject the $N$-electron system to a classical electromagnetic
radiation characterized by the fundamental electromagnetic field vectors: the
electric field $\bm{E}\equiv\bm{E}(\bm{r},t)$ and the magnetic field
$\bm{B}\equiv\bm{B}(\bm{r},t)$. These vectors satisfy the microscopic
Maxwell’s equations, which are the basic equations of motion of
electromagnetism where charged particles appear as sources. Here, we apply the
perturbation theory viewpoint of quantum electrodynamics: to first order it is
assumed that the particles of whose motion is being studied do not affect the
radiation field, which thus appear as a ”driving field” [96]. Therefore,
assuming that the sources of the radiation field are sufficiently remote from
a molecule of interest, the $\bm{E}$ and $\bm{B}$ fields are source- and
divergence-free, and conveniently described in terms of the scalar potential
$\phi\equiv\phi(\bm{r},t)$ and the vector potential
$\bm{A}\equiv\bm{A}(\bm{r},t)$, satisfying [102]:

$\begin{split}\bm{E}(\bm{r},t)&=-\bm{\nabla}\phi(\bm{r},t)-\frac{\partial\bm{A}(\bm{r},t)}{\partial
t};\qquad\bm{\nabla}\cdot\bm{E}(\bm{r},t)=0;\\\
\bm{B}(\bm{r},t)&=\bm{\nabla}\times\bm{A}(\bm{r},t);\qquad\qquad\qquad~{}\bm{\nabla}\cdot\bm{B}(\bm{r},t)=0.\end{split}$
(7)

In fact, both electromagnetic potentials enter the Dirac Hamiltonian and
describe the coupling of an electron to the classical electromagnetic field as
[98]

$\begin{split}\hat{h}^{\mathrm{D}}_{i}(t)&=\beta^{\prime}_{i}m_{e}c^{2}+c\big{(}\bm{\alpha}_{i}\cdot\bm{p}_{i}\big{)}-e\phi_{0}(\bm{r}_{i})\mathbb{I}_{4}-e\phi(\bm{r}_{i},t)\mathbb{I}_{4}+ec\big{(}\bm{\alpha}_{i}\cdot\bm{A}(\bm{r}_{i},t)\big{)}.\end{split}$
(8)

When compared to the field-free Dirac Hamiltonian in Eq. (3), the scalar
electrostatic potential due to the nuclei $\phi_{0}(\bm{r})$ is substituted by
the time-dependent potential:
$\phi_{0}(\bm{r})\rightarrow\phi_{0}(\bm{r})+\phi(\bm{r},t)$, and the
canonical momentum of an electron $\bm{p}$ is substituted by the mechanical
momentum: $\bm{p}\rightarrow\bm{p}+e\bm{A}(\bm{r},t)$. The latter substitution
is known in literature as the principle of minimal electromagnetic coupling
substitution [103].

To gain insights into the physical interpretation of the matter–field
interaction, let us consider the expectation value of the relativistic one-
electron interaction Hamiltonian given by last two terms in Eq. (8)

$\begin{split}\int\psi^{\dagger}(\bm{r},t)\Big{[}-e\phi(\bm{r},t)\mathbb{I}_{4}+ec\big{(}\bm{\alpha}\cdot\bm{A}(\bm{r},t)\big{)}\Big{]}\psi(\bm{r},t)d^{3}\bm{r}=\\\
=\int\Big{[}\rho(\bm{r},t)\phi(\bm{r},t)-\bm{j}(\bm{r},t)\cdot\bm{A}(\bm{r},t)\Big{]}d^{3}\bm{r}.\end{split}$
(9)

Assuming multiplicative potentials $\phi$ and $\bm{A}$, the second equation
reveals that the scalar potential is coupled to the electron charge density
$\rho$ – _i.e._ , the charge of the electron times its probability
distribution

$\rho(\bm{r},t)=-e\psi^{\dagger}(\bm{r},t)\mathbb{I}_{4}\psi(\bm{r},t),$ (10)

whereas the vector potential is coupled to the electron current density
$\bm{j}$ – _i.e._ the charge of the electron times its velocity distribution

$\bm{j}(\bm{r},t)=-e\psi^{\dagger}(\bm{r},t)c\bm{\alpha}\psi(\bm{r},t).$ (11)

In order to write the interaction Hamiltonian in its explicit form, we need to
know analytical expressions for both potentials $\phi$ and $\bm{A}$. By the
use of Maxwell’s equations for the source-free field, it can be shown that
these potentials satisfy [102]

$\begin{split}\nabla^{2}\phi+\frac{\partial}{\partial
t}(\bm{\nabla}\cdot\bm{A})=0,\\\
\nabla^{2}\bm{A}-\frac{1}{c^{2}}\frac{\partial^{2}\bm{A}}{\partial
t^{2}}-\bm{\nabla}\Big{(}\bm{\nabla}\cdot\bm{A}+\frac{1}{c^{2}}\frac{\partial\phi}{\partial
t}\Big{)}=0.\end{split}$ (12)

However, there exists a certain arbitrariness in the definition of the
potentials, in that it is possible to shift them by the transformation

$\phi\rightarrow\phi^{\prime}=\phi-\frac{\partial\chi}{\partial
t};\qquad\bm{A}\rightarrow\bm{A}^{\prime}=\bm{A}+\bm{\nabla}\chi,$ (13)

where $\chi\equiv\chi(\bm{r},t)$ is an arbitrary scalar function of space and
time coordinates called a gauge function. Since the physics, _i.e._ the force
law and Maxwell’s equations, is sensitive only to the electric field $\bm{E}$
and the magnetic field $\bm{B}$, the transformation of potentials, called a
gauge transformation, does not affect it. This is known in physics as gauge
invariance and may be readily verified by inserting two pairs of potentials
$(\phi,\bm{A})$ and $(\phi^{\prime},\bm{A}^{\prime})$ into the expression in
Eq. (7). In addition, gauge invariance may be exploited to simplify Eq. (12),
and this means also the interaction Hamiltonian.

In quantum chemistry, the gauge freedom is fixed by choosing the so-called
Coulomb gauge defined by the condition [103]

$\bm{\nabla}\cdot\bm{A}(\bm{r},t)=0.$ (14)

With this condition and the fact that the electric field is divergence-free in
free space (7), the scalar potential is a constant, _i.e._
$\phi(\bm{r},t)=\phi$, and may be taken as zero to satisfy $|\phi|\rightarrow
0$ at spatial infinity. In this case, the equations of motion for
electromagnetic potentials (12) simplify to

$\begin{split}\nabla^{2}\phi=\phi=0,\\\
\Big{(}\nabla^{2}-\frac{1}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}\Big{)}\bm{A}(\bm{r},t)=0.\end{split}$ (15)

The wave equation for the vector potential is identical in form in many
problems of wave motion, with a real solution in the form of a monochromatic,
linearly polarized electromagnetic plane-wave [102]

$\bm{A}(\bm{r},t)=\bm{A}_{0}\cos(\bm{k}\cdot\bm{r}-\omega t),$ (16)

where $\bm{A}_{0}$ is a constant real vector called amplitude factor. The
argument of the cosine function is called the _phase_ of $\bm{A}$ and is given
in terms of the wave vector $\bm{k}$ (characterizing the direction of wave
propagation) and the angular frequency $\omega$. Note that the phase sometimes
contains a phase constant, manipulation of which the cosine function can be
converted to a sine function. By substitution of the solution Eq. (16) back
into the wave equation Eq. (15), we find that the magnitude of $\bm{k}$ obeys
$|\bm{k}|=k=\omega/c$. In addition, noting that the angular frequency is
$\omega=2\pi\nu=2\pi c/\lambda$ with the frequency $\nu$ and wavelength
$\lambda$, we also find $|\bm{k}|=k=2\pi/\lambda$.

To summarize, time evolution of a $N$-electron system subjected to a classical
electromagnetic field is governed by the electronic Hamiltonian

$\hat{H}(t)=\sum_{i}^{N}\hat{h}^{\mathrm{D}}_{i}(t)+\frac{1}{2}\sum_{i\neq
j}^{N}\hat{g}_{ij},$ (17)

where its one-electron part given in the Dirac’s relativistic formalism as

$\begin{split}\hat{h}^{\mathrm{D}}_{i}(t)&=\beta^{\prime}_{i}m_{e}c^{2}+c\big{(}\bm{\alpha}_{i}\cdot\bm{p}_{i}\big{)}-e\phi_{0}(\bm{r}_{i})\mathbb{I}_{4}+\hat{h}^{\mathrm{(v)}}(\bm{r}_{i},t),\end{split}$
(18)

contains also the electron–field interaction Hamiltonian characterized in the
Coulomb gauge entirely by the vector potential

$\begin{split}\hat{h}^{\mathrm{(v)}}(\bm{r}_{i},t)&=ec\bm{\alpha}_{i}\cdot\bm{A}(\bm{r}_{i},t)=ec\big{(}\bm{\alpha}_{i}\cdot\bm{A}_{0})\cos(\bm{k}\cdot\bm{r}_{i}-\omega
t).\end{split}$ (19)

In the literature, $\hat{h}^{\mathrm{(v)}}$ is known as the one-electron
interaction Hamiltonian in _velocity representation_ which we shall label with
the superscript (v).

Note that the spatial phase of $\hat{h}^{\mathrm{(v)}}$ can be simplified by
considering that wavelengths of electromagnetic waves in the ultraviolet or
visible range are very large compared with the spatial extent of typical
molecular systems under study. To provide an illustrative example, let us
consider a laser pulse with wavelength $\lambda=1064\mathrm{nm}$ applied to a
molecule of size $|\bm{r}|=r=10\mathrm{\AA}$. Hence,

$\bm{k}\cdot\bm{r}\leq kr=\frac{2\pi}{1064\mathrm{nm}}1\mathrm{nm}\ll 1,$ (20)

which implies that the spatial phase of an oscillating electromagnetic wave
can be approximated by a constant over the length scale of a molecule (or more
precisely over the mean-value of an electronic position), _i.e._

$\exp{[i(\bm{k}\cdot\bm{r})]}=1+i(\bm{k}\cdot\bm{r})-\frac{1}{2}(\bm{k}\cdot\bm{r})^{2}+...\approx
1.$ (21)

Independence of the electromagnetic wave on a spatial coordinate is known as
_dipole_ (or _long-wavelength_) approximation, which brings the velocity
interaction Hamiltonian in Eq. (19) into a particularly simple form labelled
here as (vd)

$\begin{split}\hat{h}^{\mathrm{(v)}}(\bm{r}_{i},t)\approx\hat{h}^{\mathrm{(vd)}}_{i}(t)=ec\bm{\alpha}_{i}\cdot\bm{A}(t)=\frac{ec}{2}\big{(}\bm{\alpha}_{i}\cdot\bm{A}_{0}\big{)}\big{[}\exp(-i\omega
t)+\exp(i\omega t)\big{]}.\end{split}$ (22)

Here, we used $\cos(x)=[\exp(ix)+\exp(-ix)]/2$. However, special care has to
be taken for short wavelengths used for instance in hard X-ray spectroscopy
where the dipole approximation may not be adequate. In particular, this is
true for heavy-element K-edge X-ray absorption spectroscopy [104, 105]. By
including higher-order powers of $\bm{k}\cdot\bm{r}$ in the expansion, one
gets multipolar contributions known as electric-quadrupole, magnetic-dipole,
etc., and there exist techniques to include these contributions into quantum-
chemical calculations [104, 105, 106, 107].

Before we close this section, let us mention that there exists a unitary
transformation of the wave function which yields an altered form of the
interaction Hamiltonian that may be more useful for practical calculations.
Let us start from the time-dependent Schrödinger/Dirac equation with the
electronic Hamiltonian $\hat{H}$ containing the velocity-dipole interaction
Hamiltonian ($\hat{h}^{\mathrm{(vd)}}$) given by Eq. (22):

$\left(i\hbar\frac{\partial}{\partial
t}-\hat{H}\right)\Psi=0;\qquad\hat{H}\equiv\hat{H}(t)=\sum_{i}^{N}\big{[}\hat{h}^{\mathrm{D}}_{i}+\hat{h}^{\mathrm{(vd)}}_{i}(t)\big{]}+\frac{1}{2}\sum_{i\neq
j}^{N}\hat{g}_{ij}.$ (23)

The wave function $\Psi\\!\equiv\\!\Psi(t)$ can undergo a unitary (gauge)
transformation with a freely chosen function $\Lambda\equiv\Lambda(t)$

$\Psi=\exp(-i\Lambda)\Psi^{\prime}.$ (24)

The new wave function $\Psi^{\prime}\\!\equiv\\!\Psi^{\prime}(t)$ is as
physically meaningful as the old one, provided

$\left(i\hbar\frac{\partial}{\partial
t}-\hat{H}^{\prime}\right)\Psi^{\prime}=0;\qquad\hat{H}^{\prime}\equiv\hat{H}^{\prime}(t)=\exp(i\Lambda)\hat{H}(t)\exp(-i\Lambda)-\hbar\frac{\partial\Lambda}{\partial
t}.$ (25)

Now, by selecting $\Lambda$ as

$\Lambda(t)=\frac{e}{\hbar}\sum_{i}^{N}\mathbb{I}_{4}\bm{r}_{i}\cdot\bm{A}(t),$
(26)

one replaces the velocity-dipole interaction Hamiltonian in the original
electronic Hamiltonian $\hat{H}$ by a new interaction Hamiltonian
$\hat{h}^{\mathrm{(ld)}}$ in the so-called _length-dipole_ representation (ld)
in the new electronic Hamiltonian $\hat{H}^{\prime}$:

$H^{\prime}(t)=\sum_{i}^{N}\big{[}\hat{h}^{\mathrm{D}}_{i}+\hat{h}^{\mathrm{(ld)}}_{i}(t)\big{]}+\frac{1}{2}\sum_{i\neq
j}^{N}\hat{g}_{ij};\qquad\hat{h}^{\mathrm{(ld)}}_{i}(t)=e\mathbb{I}_{4}\bm{r}_{i}\cdot\bm{E}(t).$
(27)

Physically, $\hat{h}^{\mathrm{(ld)}}$ couples the classical electric field
$\bm{E}(t)$ defined as

$\bm{E}(t)=-\frac{\partial}{\partial t}\bm{A}(t),$ (28)

to the quantum-mechanical system characterized by the sum of the electric
dipole moment operators of individual electrons
($\bm{\mu}_{i}=-e\bm{r}_{i}\mathbb{I}_{4}$). This gauge transformation was
first discussed at the nonrelativistic level by Göppert-Mayer in 1931 [108]
and therefore it is often named after her. An additional reading on quantum-
mechanical gauge invariance and general unitary transformations for atoms and
molecules in interactions with radiation can be found in Ref. [109].

##  4 Equations-of-motion for exact-state wave function

###  4.1 Time-dependent Schrödinger equation

In the most general case, the time evolution of a quantum-mechanical system is
governed by the time-dependent equation-of-motion

$i\hbar\frac{\partial\Psi(t)}{\partial t}=\hat{H}(t)\Psi(t),$ (29)

where $\hat{H}(t)$ is the Hamiltonian operator that is explicitly dependent on
time via external electromagnetic fields. If we are interested in the response
of molecules subjected to ultrafast laser pulses and similar processes
occurring at atto- or femtosecond time scales, we can study electron dynamics
decoupled from the nuclear motion. However, molecular vibrations and nuclear
relaxation occur on a time scale of 10–100 fs, and in principle should be
included in the computation. Nevertheless, performing electron dynamics
simulations with fixed nuclear configuration at this time scale is still
beneficial for improving the spectral resolution and aids analyzing electron
excitations without the effect of nuclear dynamics. In such a case, $\Psi(t)$
in Eq. (29) is the many-electron wave function depending on the position
$\bm{r}_{i}$ and spin of all electrons, and the Hamiltonian $\hat{H}$ is the
many-electron Hamiltonian defined in Eq. (17) containing the electron kinetic
operator, Coulomb interactions between electrons and nuclei with fixed
positions, and interactions between the system and external electromagnetic
fields discussed in more detail in Section 3.

Eq. (29) also remains valid in relativistic case, provided additional
approximations are assumed. The Hamiltonain $\hat{H}$ needs to be treated as a
multicomponent operator acting on multicomponent wave functions to reflect the
fact that in relativistic theory, electron spin and orbital degrees of freedom
interact with each other via the spin–orbit coupling terms. However, in a
truly relativistic picture, we would need to consider multiple time variables
associated with each electron’s frame of reference. Such effects arising from
the relative time are always neglected when studying molecular systems, and
Eq. (29) thus assumes the absolute time approximation, which leads to a single
time variable $t$. For further discussion on the relativistic theory of many
electrons, see Ref. [101].

Response of the system to external time-dependent electromagnetic fields can
be studied by solving Eq. (29). This can be achieved by using the formalism of
response theory [110, 5], in case the external fields are weak and can be
regarded as small perturbations to the system compared to the intrinsic
unperturbed Hamiltonian. Alternatively, the equation can be solved numerically
by propagating the wave function in real time, which facilitates studying
processes that involve arbitrarily strong fields.

###  4.2 Reduced density matrices

Since the many-electron wave function is a complicated object that depends on
the spatial coordinates of each electron, for the forthcoming discussion, it
will be more convenient to work in the formalism of reduced density matrices
(RDMs) [111]. In the time domain, we can define the one-electron and two-
electron RDMs, respectively, as

$D(\bm{r}_{1};\bm{r}^{\prime}_{1};t)=N\int\Psi(\bm{r}_{1},x_{2},\ldots,x_{N},t)\Psi^{\dagger}(\bm{r}^{\prime}_{1},x_{2},\ldots,x_{N},t)dx_{2}\ldots
dx_{N},$ (30)

and

$\Gamma(\bm{r}_{1},\bm{r}_{2};\bm{r}^{\prime}_{1},\bm{r}^{\prime}_{2};t)=N(N-1)\int\Psi(\bm{r}_{1},\bm{r}_{2},x_{3},\ldots,x_{N},t)\\\
\times\Psi^{\dagger}(\bm{r}^{\prime}_{1},\bm{r}^{\prime}_{2},x_{3},\ldots,x_{N},t)dx_{2}\ldots
dx_{N},$ (31)

where $N$ is the number of electrons. We note here, that whereas $\bm{r}_{i}$
represents spatial coordinates in three-dimensional space,
$x_{i}\equiv(\bm{r}_{i},\tau_{i})$ denotes both the position $\bm{r}_{i}$
_and_ the spin $\tau_{i}$ of the $i$-th electron, and the integration
symbolically also labels the summation over the spin degrees of freedom in
addition to the integration over the spatial variables. In the relativistic
theory with SOC, it is convenient to keep the indices associated with
$\bm{r}_{1}$ and $\bm{r}_{2}$ free. As a consequence, $D$ and $\Gamma$ are
still multicomponent tensors, for instance, in case of the Dirac theory, $D$
and $\Gamma$ have the dimensions of $4\times 4$ and $4\times 4\times 4\times
4$, respectively. Hence, the scalar electron charge density is obtained as

$\rho(\bm{r},t)=-e\operatorname{Tr}D(\bm{r};\bm{r};,t),$ (32)

where $\operatorname{Tr}$ indicates the trace over the bispinor components.
Likewise, for the four-component current density, it follows that

$\bm{j}(\bm{r},t)=-ec\operatorname{Tr}\bm{\alpha}D(\bm{r};\bm{r};,t).$ (33)

Eqs. (32) and (33) generalize the one-electron definitions of the charge and
current densities in Eqs. (10) and (11) for many-electron wave functions,
since they are agnostic to the method that was used to calculate the one-
electron RDM.

Exact time propagation determined by Eq. (29) can equivalently be formulated
in the language of RDMs, which avoids the use of the cumbersome many-electron
wave function. Let us assume that we have a set of orthonormal spin-orbitals
$\varphi_{p}(\bm{r})$. The one- and two- RDMs matrices in the spin-orbital
basis then become

$D_{pq}(t)=\int d^{3}\bm{r}_{1}\int
d^{3}\bm{r}^{\prime}_{1}\varphi^{\dagger}_{p}(\bm{r}_{1})D(\bm{r}_{1};\bm{r}^{\prime}_{1};t)\varphi_{q}(\bm{r}^{\prime}_{1}),$
(34)

and

$\Gamma_{pqrs}(t)=\int d^{3}\bm{r}_{1}\ldots
d^{3}\bm{r}^{\prime}_{2}\varphi^{\dagger}_{p}(\bm{r}_{1})\varphi^{\dagger}_{r}(\bm{r}_{2})\Gamma(\bm{r}_{1},\bm{r}_{2};\bm{r}^{\prime}_{1},\bm{r}^{\prime}_{2};t)\varphi_{q}(\bm{r}^{\prime}_{1})\varphi_{s}(\bm{r}^{\prime}_{2}),$
(35)

respectively. The time-dependent one-electron RDM can be obtained by solving
the equation of motion of Liouville-von Neumann (LvN) type [112]

$i\hbar\frac{\partial}{\partial
t}\mathbf{D}(t)=\left[\mathbf{h}(t),\mathbf{D}(t)\right]+\frac{1}{2}\operatorname{Tr}_{1}\left[\mathbf{G},\mathbf{\Gamma}(t)\right],$
(36)

where $[,]$ denotes the commutator, $\mathbf{h}(t)$ and $\mathbf{G}$ are
matrices of one- and antisymmetrized two-electron integrals

$\displaystyle
h_{pq}(t)\equiv\int\varphi^{\dagger}_{p}(\bm{r})\hat{h}^{\text{D}}(t)\varphi_{q}(\bm{r})d^{3}\bm{r},$
(37) $\displaystyle
G_{pqrs}\equiv\mathcal{I}_{pqrs}-\mathcal{I}_{psrq};\quad\mathcal{I}_{pqrs}\equiv\iint\varphi^{\dagger}_{p}(\bm{r}_{1})\varphi_{q}(\bm{r}_{1})r^{-1}_{12}\varphi^{\dagger}_{r}(\bm{r}_{2})\varphi_{s}(\bm{r}_{2})d^{3}\bm{r}_{1}d^{3}\bm{r}_{2},$
(38)

and

$\displaystyle(\mathbf{G\Gamma})_{pqrs}$ $\displaystyle\equiv
G_{perf}\Gamma_{eqfs},$ (39)
$\displaystyle(\operatorname{Tr}_{1}\mathbf{X})_{pq}$ $\displaystyle\equiv
X_{pqrr}$ (40)

for any two-electron matrix $\mathbf{X}$. Upon inspecting Eq. (36), we can see
that the exact time evolution of the one-electron RDM also depends on the two-
electron RDM $\mathbf{\Gamma}(t)$, which is also not known. Likewise, we could
proceed by writing the equation of motion for the two-electron RDM. However,
in general, the equation of motion for the $N$-electron RDM will contain the
RDM of the order $N+1$, leading to an infinite hierarchy of coupled equations
for RDMs, mirroring the same situation that occurs in the theory of Green’s
functions [113]. Solving the resulting system of equations is impractical,
hence, approximations to the higher-order second term that decouple the
equations are sought. In the following sections, we will describe the LvN
equation for the one-electron RDM where the second term containing
$\mathbf{\Gamma}(t)$ is approximated in the mean-field manner using only the
one-electron RDM in the framework of time-dependent Hartree–Fock theory and
density functional theory.

###  4.3 Time-reversal symmetry

One of the most important properties of quantum-mechanical equations of motion
(and all microscopic laws) is their symmetry with respect to the reversal of
time. Let us use the shorthand notation for the many-electron wave function
$\Psi(t)\equiv\Psi(x_{1},\ldots,x_{N},t)$. Replacing $t\rightarrow-t$ in Eq.
(29) gives

$-i\hbar\frac{\partial\Psi(-t)}{\partial t}=\hat{H}(-t)\Psi(-t).$ (41)

This equation differs from the original one in two ways. First, the
Hamiltonian is expressed in the inverted time $-t$. Second, there is an extra
minus sign on the left hand side of the equation. Let us assume we have an
_antiunitary_ operator $\mathcal{K}$ that is unitary
($\mathcal{K}^{\dagger}\mathcal{K}=\mathbb{I}$) and antilinear

$\mathcal{K}i=-i\mathcal{K}.$ (42)

Letting this operator act from the left on the Eq. (41), and denoting

$\displaystyle\bar{\Psi}(t)$ $\displaystyle:=\mathcal{K}\Psi(-t),$ (43)
$\displaystyle\bar{H}(t)$
$\displaystyle:=\mathcal{K}\hat{H}(-t)\mathcal{K}^{\dagger},$ (44)

we obtain

$i\hbar\frac{\partial\bar{\Psi}(t)}{\partial t}=\bar{H}(t)\bar{\Psi}(t).$ (45)

In principle, this is a new equation of motion with a new solution, however,
if we can assume that the Hamiltonian satisfies $\bar{H}(t)=\hat{H}(t)$,
_i.e._

$\mathcal{K}\hat{H}(t)=\hat{H}(-t)\mathcal{K},$ (46)

then Eqs. (29) and (45) represent the _same_ equation, for which we obtained a
pair of solutions $\Psi(t)$ and $\bar{\Psi}(t)$. The condition in Eq. (46) is
known as time-reversal symmetry (TRS).

Due to the requirement in Eq. (42), the operator $\mathcal{K}$ must _at least_
contain complex conjugation. This is a sufficient condition for scalar wave
functions in nonrelativistic theory, where Eq. (46) reduces to
$\hat{H}^{*}(t)=\hat{H}(-t)$ and additionally the condition that the time-
independent Hamiltonian is real-valued. However, for spinor wave functions and
multicomponent relativistic theories, $\mathcal{K}$ can have a more
complicated matrix form. For instance, in case of the Dirac four-component
one-electron Hamiltonian, the operator $\mathcal{K}$ takes the form [100, 114,
115]

$\mathcal{K}=-i\begin{pmatrix}\sigma_{y}&0_{2}\\\
0_{2}&\sigma_{y}\end{pmatrix}\mathcal{K}_{0},$ (47)

where $\mathcal{K}_{0}$ denotes the complex conjugation, and $\sigma_{y}$ is
the $y$-th $2\times 2$ Pauli matrix.

We conclude this section by noting that the condition in Eq. (46) is satisfied
for nonrelativistic as well as relativistic Hamiltonians. Neither internal
electromagnetic interactions nor spin–orbit coupling terms break TRS, _i.e._
the Hamiltonain consists of symmetric operators
($\mathcal{K}\hat{A}\mathcal{K}^{\dagger}=\hat{A}$) or bilinear products of
antisymmetric ($\mathcal{K}\hat{A}\mathcal{K}^{\dagger}=-\hat{A}$) operators,
such as $\bm{\sigma}\cdot\bm{p}$, that are again symmetric. However, if
external fields are introduced, they can break the TRS, for instance, an
electric field with the time dependence given by an odd function of $t$. More
importantly, it is often discussed in the literature that the presence of a
magnetic field breaks TRS. This is only true if the magnetic field $\bm{B}$ is
considered as external and does not change its orientation upon time reversal,
_i.e._ $\mathcal{K}$ only acts on the electronic degrees of freedom. In such
situations, terms like $\bm{B}\cdot\hat{\bm{S}}$, where $S$ denotes the
electron spin, become antisymmetric with respect to the time reversal, because
the operator $\mathcal{K}$ only acts on $\hat{\bm{S}}$
($\mathcal{K}\hat{\bm{S}}\mathcal{K}^{\dagger}=-\hat{\bm{S}}$) and not on
$\bm{B}$.

##  5 Equations-of-motion for approximate-state wave functions

The previous section dealt with exact state theory. In practical calculations,
model quantum chemistries are used to treat systems containing many particles.
The theory presented here focuses on both time-dependent Hartree–Fock (TDHF)
theory and time-dependent density functional theory (TDDFT) in the time-
dependent Kohn–Sham (TDKS) framework. From the practical point of view, both
TDHF and KS TDDFT are mean-field theories solving equations for one-electron
molecular orbitals. Therefore, we use the term time-dependent self-consistent
field (TDSCF) when referring to both methods together. In the following text,
we sketch the derivation of working equations for TDHF and TDKS theories.
Since the final form of the equations is the same for both methods, the rest
of this chapter concerning propagators, evaluation of molecular properties,
and analysis applies equally to both of them.

###  5.1 Time-dependent Hartree–Fock theory

The main idea of the TDHF method is to approximate the many-electron time-
dependent wave function $\Psi(x_{1},x_{2},\ldots;t)$ by a single Slater
determinant built from time-dependent molecular spin-orbitals (MO)
$\varphi_{i}(x,t)$, where we grouped the electron’s spatial and spin degrees
of freedom into a single variable $x\equiv(\bm{r},\tau)$. Hence, the ansatz
reads

$\displaystyle\Psi(x_{1},x_{2},\ldots;t)$
$\displaystyle=\frac{1}{\sqrt{N!}}\begin{vmatrix}\varphi_{1}(x_{1},t)&\varphi_{1}(x_{2},t)&\cdots&\varphi_{1}(x_{N},t)\\\
\varphi_{2}(x_{1},t)&\varphi_{2}(x_{2},t)&\cdots&\varphi_{2}(x_{N},t)\\\
\vdots&\vdots&\ddots&\vdots\\\
\varphi_{N}(x_{1},t)&\varphi_{N}(x_{2},t)&\cdots&\varphi_{N}(x_{N},t)\\\
\end{vmatrix}$
$\displaystyle=\frac{1}{\sqrt{N!}}\sum_{\\{P\\}}(-1)^{|P|}\left|\varphi_{P(1)}(x_{1},t)\varphi_{P(2)}(x_{2},t)\ldots\varphi_{P(N)}(x_{N},t)\right|,$
(48)

where $P$ denotes a permutation of indices, $P(i)$ is the new index after
permutation, and $\sum_{\\{P\\}}$ is the sum over all possible permutations of
MO indices. The prefactor $(-1)^{|P|}$ is the sign $\pm 1$ of the permutation
based on the permutation length $|P|$. This ansatz uses complex spin-orbitals
instead of real scalar orbitals and facilitates a direct extension of the
nonrelativistic HF theory into the relativistic domain. Furthermore, in the
TDHF theory, we assume that the many-electron wave function retains the form
of a single Slater-determinant during the entire time evolution.

The working equations of TDHF can be derived using the time-dependent
variational principle. Several functionals to be minimized have been
formulated, such as the Dirac–Frenkel functional [116, 5]

$I^{\mathrm{DF}}=\int dt\Braket{\Psi(t)}{i\hbar\frac{\partial}{\partial
t}-\hat{H}}{\Psi(t)},$ (49)

or the McLachlan functional

$I^{\mathrm{ML}}(t)=\Braket{\left(i\hbar\frac{\partial}{\partial
t}-\hat{H}\right)\Psi(t)}{\left(i\hbar\frac{\partial}{\partial
t}-\hat{H}\right)\Psi(t)},$ (50)

where we used the braket notation $\Braket{\ldots}$ to indicate the
integration over degrees of freedom (spin and spatial) of all electrons. These
functionals can be used to derive the final form of the TDHF equations for MOs
[117]

$i\hbar\frac{\partial}{\partial
t}\varphi_{i}(\bm{r},t)=\hat{F}_{\text{HF}}\left[\\{\varphi_{j}(\bm{r},t)\\}\right]\\!(\bm{r},t)~{}\varphi_{i}(\bm{r},t),$
(51)

where $\hat{F}_{\text{HF}}$ is the Fock operator known from time-independent
HF theory. Here, $\hat{F}_{\text{HF}}$ contains the Coulomb interaction of the
electron with the mean-field of other electrons, the Fock exchange operator,
and the one-electron Dirac Hamiltonian $\hat{h}^{\text{D}}(t)$ that includes
the interaction with external time-dependent electromagnetic fields. As a
consequence, in the four-component Dirac theory, $\hat{F}_{\text{HF}}$ is a
$4\times 4$ operator acting on bispinor orbitals $\phi_{i}$. The presence of
the explicitly time-dependent external fields in the one-electron part of
$\hat{F}_{\text{HF}}$ and the dependence of the mean-field and exchange terms
on the spin-orbitals, which are now time-dependent, represents the most
distinct difference of the Fock operator in the TDHF theory from its time-
independent counterpart.

###  5.2 Time-dependent Kohn–Sham DFT

Analogously to the static case, the idea of time-dependent density functional
theory (TDDFT)[118] is to replace the many-electron wave function of $3N$
spatial variables and time with the simpler object — electron density
$\rho(\bm{r},t)$. In the nonrelativistic framework, the theoretical
foundations of TDDFT are provided by two theorems. The first one, the
Runge–Gross theorem [119], is a time-dependent analogue of the Hohenberg–Kohn
theorem connecting the time-dependent external potential and time-dependent
density. The second one, the Van Leeuwen theorem [120], connects the real
system with a fictitious system with different interaction potential.
Application of these theorems allows for introducing a fictitious KS system of
non-interacting electrons for which the many-electron wave function is a
single Slater determinant built of one-particle functions called KS orbitals.
The final TDKS equations are similar to Eq. (51),

$i\hbar\frac{\partial}{\partial
t}\varphi_{i}(\bm{r},t)=\hat{F}_{\text{KS}}\left[\\{\varphi_{j}(\bm{r},t)\\}\right]\\!(\bm{r},t)~{}\varphi_{i}(\bm{r},t),$
(52)

except that the Fock operator $\hat{F}_{\text{KS}}$ also contains the
exchange–correlation (XC) potential derived from the approximation to the XC
energy functional instead of the exact HF exchange term. This XC term links
the fictitious KS system to the studied real system. In hybrid DFT [121],
functionals allow for a fraction of the HF exchange contribution to also enter
$\hat{F}_{\text{KS}}$, bringing the very important element of (exact)
antisymmetry of the many-electron wave function to DFT.

Approximating the XC functional in the time domain is more challenging than in
time-independent theory. In principle, TDDFT requires the development and use
of special time-dependent XC potentials that may generally depend on the
density in previous times. However, a widespread practice is to simply use
potentials from time-independent DFT, with the time variable only entering via
the time-dependence of the density (and its gradient). This local-in-time
approximation is known as the adiabatic approximation in TDDFT [122, 123]. The
term adiabatic approximation is actually used to label a combination of two
approximations: firstly the adiabatic approximation itself [118] and secondly
the approximations that were used in the construction of the time independent
XC functional [124]. The adiabatic approximation is valid when a system
remains in its instantaneous eigenstate for slowly varying perturbations that
act on it [118] and is widely used in TDDFT due to the lack of accepted time
(memory) dependent functionals. The memory effects were also shown to be
negligible in the context of nonlinear processes and strong-field excitations
studied in non-perturbative electron dynamics [125]. However, non-adiabatic
effects in the XC functional become important for high-frequency oscillations
[125], double and charge-transfer excitations [126]. Extending the XC
potential beyond the adiabatic approximation while still exploiting the local
gradient expansion can be achieved if the current density is used as a central
variable [127, 128]. More general framework for time-dependent functionals
with memory in TDDFT introduces viscoelastic stresses known in hydrodynamics
for the electron liquid [129, 126] or formulates TDDFT in a comoving
Lagrangian reference frame [130].

In a similar manner, the extension of DFT to the relativistic domain [131,
132] also makes use of non-relativistic XC potentials that take relativistic
densities as input. Linear-response TDDFT has been extended to the
relativistic approximate two-component framework and applied to calculate
absorption spectra of solids [133], however, proper theoretical foundations
that incorporate both effects of time-dependent fields as well as relativity
and generalize the Runge–Gross and Van Leeuwen theorems to the relativistic
domain do not exist. Despite this, the relativistic real-time and linear-
response TDDFT has been applied to study a number of molecular properties, as
will be discussed in more detail in Section 8.

###  5.3 Liouville–von Neumann equation in four-component framework

As discussed in Section 4.2, it is often more practical to work in the
formalism of density matrices. This is especially the case for theories that
express the many-electron wave function as a single Slater determinant, such
as TDHF and TDKS, where the two-electron RDM is not needed, and the entire
information about the time evolution of a quantum state of the many-electron
system is encoded in the one-electron RDM.

Let us express the time-dependent spin-orbitals $\varphi_{i}(\bm{r},t)$
appearing in Eqs. (51) and (52) using a set of $n$ static _orthonormal_
functions $\\{X(\bm{r})\\}$. Then

$\varphi_{i}(\bm{r},t)=\sum_{\mu}^{n}X_{\mu}(\bm{r})C_{\mu i}(t),$ (53)

where $C_{\mu i}(t)$ are the complex-valued expansion coefficients. For
purposes of this Chapter, $\\{X(\bm{r})\\}$ shall refer to orthonormal atomic
orbitals (AOs). For cases where the wave function $\Psi(t)$ is a Slater
determinant, the one-electron RDM from Eq. (30) can be expressed through the
occupied (occ) spin orbitals as

$D(\bm{r};\bm{r}^{\prime};t)=\sum_{i}^{\text{occ}}\varphi_{i}(\bm{r},t)\varphi^{\dagger}_{i}(\bm{r}^{\prime},t).$
(54)

Inserting Eqs. (53) and (56) into Eq. (34) gives the following orthonormal AO
representation of the RDM

$D_{\mu\nu}(t)=\sum_{i}^{\text{occ}}C_{\mu i}(t)C^{\dagger}_{\nu i}(t).$ (55)

Introducing the matrices $\mathbf{D}(t)$ and $\mathbf{C}(t)$ with elements
$D_{\mu\nu}(t)$ and $C_{\mu i}$, respectively, we can write

$\mathbf{D}(t)=\mathbf{C}(t)\mathbf{C}^{\dagger}(t).$ (56)

Taking the time derivative of this equation and using the time-dependent
equations for spin-orbitals (Eqs. (51) or (52)) in combination with the
expansion in Eq. (53), we obtain the Liouville–von Neumann (LvN) equation of
motion (EOM) for the RDM

$i\hbar\frac{\partial\mathbf{D}(t)}{\partial
t}=[\mathbf{F}(t),\mathbf{D}(t)],$ (57)

where we dropped the HF and KS labels on the Fock matrix $\mathbf{F}(t)$. We
note, that in the HF theory, this equation coincides with the general Eq.
(36), since the two-electron RDM is approximated as

$\Gamma_{\mu\nu\kappa\lambda}(t)=D_{\mu\nu}(t)D_{\kappa\lambda}(t)-D_{\mu\lambda}(t)D_{\kappa\nu}(t).$
(58)

Using this factorization of the two-electron RDM in Eq. (36) gives rise to
both the mean-field Coulomb as well as the exact exchange terms that
complement the one-electron Hamiltonian in the Fock matrix of HF theory. The
advantage of solving the LvN equation over the respective equations for spin
orbitals is that the RDM is gauge invariant with respect to orbital rotations
$\varphi^{\prime}_{p}\rightarrow\varphi^{\prime}_{q}V_{qp}$, _i.e._ unitary
transformations $\mathbf{V}$ that do not mix the occupied and virtual spin
orbitals. This gauge freedom of the orbitals was utilized in the work of Jia
_et al._ [134] to allow for much larger time steps used in real-time
simulations based on solving the EOM for orbitals in the parallel transport
gauge.

Real-time methods are based on directly solving Eq. (57) in the time domain by
numerically propagating the RDM (see Section 6). Since the Fock matrix is
Hermitian, the time evolution must be unitary. However, Eq. (57) is sometimes
augmented by an extra term to model the relaxation of the system to the
equilibrium (ground) state $\mathbf{D}_{eq}$ with an empirical rate of
relaxation matrix $\gamma$. The LvN equation then reads

$i\hbar\frac{\partial\mathbf{D}(t)}{\partial
t}=[\mathbf{F}(t),\mathbf{D}(t)]-i\hbar\mathbf{\gamma}\big{(}\mathbf{D}(t)-\mathbf{D}_{eq}\big{)}.$
(59)

At the level of theory presented here, the matrix $\gamma$ is phenomenological
and is commonly approximated by a single parameter, referred to as a damping
parameter. In this case, the LvN equation can be solved without the damping
parameter and its application is postponed to a post-processing step (see
discussion below). We note, that if the damping term is included in the LvN
equation, the time propagation is no longer unitary and the energy of the
system is not conserved even in the absence of external field(s).

Within a finite time window, the solution of the LvN equation Eq. (57) reduces
to the evaluation of the time-dependent Fock matrix at discrete time steps,
and to the propagation of the density matrix in time. Here, we briefly outline
the main features of the Fock matrix evaluation, assuming the full four-
component level of theory. By following the previous discussion, the Fock
matrix in Eq. (57) is given in a set of $n$ _orthonormal_ AOs
$\\{X(\bm{r})\\}$,

$F_{\mu\nu}(t)=\Braket{X_{\mu}(\bm{r})}{\hat{F}\left[\\{\varphi(\bm{r},t)\\}\right]\\!(\bm{r},t)}{X_{\nu}(\bm{r})}.$
(60)

Of particular interest in this Chapter are applications where molecular
systems are irradiated by classical time-dependent electric field(s). In this
case, the 4c Fock matrix can easily be derived from the electronic Hamiltonian
in the length-dipole representation (see Eq. (27) in Section 3) [24, 27]

$\displaystyle
F^{\text{4c}}_{\mu\nu}(t)=F^{\mathrm{4c}}_{\mu\nu}[\bm{\mathcal{E}},\bm{\mathcal{F}}](t)$
$\displaystyle=h^{\mathrm{D}}_{\mu\nu}+\sum_{\kappa\lambda}^{n}G^{\text{4c}}_{\mu\nu,\kappa\lambda}D^{\text{4c}}_{\lambda\kappa}(t,\bm{\mathcal{E}},\bm{\mathcal{F}})$
(61) $\displaystyle+\sum_{u\in 0,x,y,z}\int
v^{xc}_{u}\\!\left[\bm{\rho}^{\text{4c}}(\bm{r},t,\bm{\mathcal{E}},\bm{\mathcal{F}})\right]\Omega_{u,\mu\nu}^{\text{4c}}(\bm{r})\,d^{3}\bm{r}$
$\displaystyle-\sum_{u\in
x,y,z}P^{\mathrm{4c}}_{u,\mu\nu}\mathcal{E}_{u}(t)-\sum_{u\in
x,y,z}P^{\mathrm{4c}}_{u,\mu\nu}\mathcal{F}_{u}(t).$

The right-hand side includes the matrix representation of the one-electron
Dirac operator, the two-electron (2e) Coulomb interaction operator, the
exchange–correlation (xc) operator, and the particle–field interaction
operators. For generality we involve two time-dependent electric fields
$\bm{\mathcal{E}}(t)$ and $\bm{\mathcal{F}}(t)$ which are coupled to the
molecular system via the electric dipole moment operator matrix
($\mathbf{P}^{\mathrm{4c}}_{u}$).

Computationally most demanding is the 2e contribution as it requires the
evaluation of generalized anti-symmetrized electron repulsion integrals (ERIs)
[82]

$G^{\text{4c}}_{\mu\nu,\kappa\lambda}=\mathcal{I}^{\text{4c}}_{\mu\nu,\kappa\lambda}-\zeta\mathcal{I}^{\text{4c}}_{\mu\lambda,\kappa\nu};\qquad\mathcal{I}^{\text{4c}}_{\mu\nu,\kappa\lambda}=\iint\Omega_{0,\mu\nu}^{\text{4c}}(\bm{r}_{1})r_{12}^{-1}\Omega_{0,\kappa\lambda}^{\text{4c}}(\bm{r}_{2})d^{3}\bm{r}_{1}d^{3}\bm{r}_{2},$
(62)

in terms of 4c _charge_ distribution functions

$\Omega^{\text{4c}}_{0,\mu\nu}(\bm{r})=X_{\mu}^{\dagger}(\bm{r})X_{\nu}(\bm{r}).$
(63)

Here, each 4c basis function
$X_{\mu}(\bm{r})\equiv\\{X^{\text{L}}_{\mu}(\bm{r})\oplus
X^{\text{S}}_{\mu}(\bm{r})\\}$ consists of the direct product of the large 2c
function $X^{\text{L}}_{\mu}(\bm{r})$ and the small 2c function
$X^{\text{S}}_{\mu}(\bm{r})$, related to each other to the lowest order in
$c^{-1}$ by the restricted kinetically balanced (RKB) relation [135]:
$X^{\text{S}}_{\mu}\simeq(\bm{\sigma}\cdot\bm{p}X^{\text{L}}_{\mu})$. Obvious
computational cost and complexity of 4c ERIs arise from the presence of the
canonical momentum operator ($\bm{p}$) as well as the Pauli spin operator
($\bm{\sigma}$) in the expression for the small-component basis. Therefore, as
discussed in Ref. [82], a single 4c ERI requires even in the most compact
formalism of real quaternions the simultaneous evaluation and processing of 25
times more real scalar integrals than the simpler 1c or 2c cases. This ratio
further increases when RKB is substituted by the restricted magnetically
balanced (RMB) relation [136, 137], which is needed for handling interactions
with magnetic fields and requires the ERI evaluation formalism to be based on
complex quaternions [82].

In addition to the charge distribution function
$\bm{\Omega}_{0}^{\text{4c}}(\bm{r})$ used in Eq. (62), one can define three
_spin_ distribution functions along the Cartesian directions

$\Omega^{\text{4c}}_{k,\mu\nu}(\bm{r})=X_{\mu}^{\dagger}(\bm{r})\Sigma_{k}X_{\nu}(\bm{r});\qquad\Sigma_{k}=\begin{pmatrix}\sigma_{k}&0_{2}\\\
0_{2}&\sigma_{k}\end{pmatrix};\qquad k\in x,y,z,$ (64)

in terms of which the 4c electron charge density ($\rho_{0}^{\text{4c}}$) as
well as the electron spin densities
($\rho_{x}^{\text{4c}},\rho_{y}^{\text{4c}},\rho_{z}^{\text{4c}}$) have a
particularly simple form

$\rho_{k}^{\text{4c}}=\rho_{k}^{\text{4c}}(\bm{r},t)=\sum_{\mu\nu}^{n}\Omega_{k,\mu\nu}^{\text{4c}}(\bm{r})D^{\text{4c}}_{\nu\mu}(t);\qquad
k\in 0,x,y,z,$ (65)

where $\Sigma_{k}$ is the Dirac spin operator. Note that all current
_noncollinear_ extensions of nonrelativistic xc functionals employ those four
densities (alongside of their gradients) as basic variables [138, 139, 140,
141, 142]. In the relativistic 2c and 4c theory, the use of a noncollinear
formalism is necessary since the spatial and spin degrees of freedom are no
longer independent and are coupled by the spin-orbit interaction. This
coupling results in a lack of rotational invariance of the xc energy if the
energy is calculated collinearly through the $z$ spin-component only [140]. A
common way to circumvent this variance problem is to formulate the
nonrelativistic exchange–correlation functionals noncollinearly. Therefore, we
utilize in our real-time TDSCF implementation the noncollinear variables of
Scalmani and Frisch [141] and evaluate the noncollinear xc potential
$v_{k}^{xc}$ in Eq. (61) within a generalized gradient approximation as

$v^{xc}_{k}\\!\left[\bm{\rho}^{\text{4c}}(t)\right]=\frac{\partial\varepsilon^{xc}}{\partial\rho_{k}^{\text{4c}}(t)}-\left(\bm{\nabla}\cdot\frac{\partial\varepsilon^{xc}}{\partial\bm{\nabla}\rho_{k}^{\text{4c}}(t)}\right);\qquad
k\in 0,x,y,z.$ (66)

Here, $\varepsilon^{xc}$ and $\bm{\rho}^{\text{4c}}$ refer to a
nonrelativistic xc energy density and an electron density vector consisting of
the electron charge and spin densities (together with their gradients). For
further details on our noncollinearity implementation, the reader is referred
to Refs. [82, 142].

###  5.4 Reduction of the Liouville–von Neumann equation to the exact two-
component (X2C) form

While full four-component (4c) relativistic real-time electron dynamics
simulations are nowadays feasible [24, 25, 26, 27, 84, 83], there is interest
in developing approximate methods enabling these simulations to be performed
more efficiently at the two-component (2c) level while maintaining the
accuracy of the parent 4c regime. Therefore, we shall discuss the
transformation of the original 4c Liouville-von Neumann (LvN) equation to its
2c form, with a particular focus on the modern exact two-component (X2C)
formalism.

The X2C Hamiltonian model has gained wide popularity in recent years as it
reduces the original 4c problem by half while requiring only a few simple
algebraic manipulations [85, 86, 87, 88, 89]. However, accuracy of this
Hamiltonian strongly depends on the two-electron (2e) and exchange–correlation
(xc) picture-change correction models employed [91] and can vary as much as
5-6 orders of magnitude for core-shell energies. Since the pioneering X2C RT-
TDDFT implementations [26, 90] utilize a crude one-electron X2C (1eX2C)
Hamiltonian model where the picture-change corrections are entirely neglected,
the inner-shell spinors (and their energies) substantially differ from the
reference 4c results [91]. Therefore, our focus here is to provide theoretical
insights into three numerically accurate X2C Hamiltonian models [91, 27, 28],
dubbed as amfX2C, eamfX2C and mmfX2C that enable accounting for the two-
electron and exchange-correlation picture-change effects.

By following the matrix-algebraic approach of X2C, let us assume that at an
arbitrary time $t$ there exists a unitary transformation matrix
$\mathbf{U}(t)$ that block-diagonalizes/decouples the 4c Fock matrix

$\mathbf{F}^{\text{4c}}(t)\rightarrow\tilde{\mathbf{F}}^{\text{4c}}(t)=\mathbf{U}^{\dagger}(t)\mathbf{F}^{\text{4c}}(t)\mathbf{U}(t)=\left(\begin{array}[c]{cc}\mathbf{\tilde{F}}^{\mathrm{LL}}(t)&\bm{0}_{2}\\\
\bm{0}_{2}&\mathbf{\tilde{F}}^{\mathrm{SS}}(t)\\\ \end{array}\right).$ (67)

Note that: (i) we use tildes to indicate all transformed quantities; (ii)
$\mathbf{F}^{\text{4c}}(t)$ and $\mathbf{U}(t)$ also depend on the electric
field $\bm{\mathcal{E}}(t)$ and $\bm{\mathcal{F}}(t)$, but for clarity of
presentation this dependence is omitted now. Under the X2C transformation, the
parent 4c EOM for MO coefficients becomes

$i\hbar\frac{\partial\bm{\tilde{C}}^{\mathrm{4c}}_{i}(t)}{\partial
t}=\mathbf{\tilde{F}}^{\mathrm{4c}}(t)\bm{\tilde{C}}^{\mathrm{4c}}_{i}(t)+i\hbar\left(\frac{\partial\mathbf{U}^{\dagger}(t)}{\partial
t}\right)\mathbf{U}(t)\bm{\tilde{C}}^{\mathrm{4c}}_{i}(t),$ (68)

where

$\bm{\tilde{C}}^{\mathrm{4c}}_{i}(t)=\mathbf{U}^{\dagger}(t)\bm{C}^{\mathrm{4c}}_{i}(t).$
(69)

A similar relation also holds for the X2C transformed LvN equation

$i\hbar\frac{\partial\mathbf{\tilde{D}}^{\mathrm{4c}}(t)}{\partial
t}=\left[\mathbf{\tilde{F}}^{\mathrm{4c}}(t),\mathbf{\tilde{D}}^{\mathrm{4c}}(t)\right]+i\hbar\left[\left(\frac{\partial\mathbf{U}^{\dagger}(t)}{\partial
t}\right)\mathbf{U}(t),\mathbf{\tilde{D}}^{\mathrm{4c}}(t)\right],$ (70)

with the density matrix

$\tilde{\mathbf{D}}^{\text{4c}}(t)=\sum_{i}^{\text{occ}}\bm{\tilde{C}}^{\mathrm{4c}}_{i}(t)\left(\bm{\tilde{C}}^{\mathrm{4c}}_{i}(t)\right)^{\dagger}=\mathbf{U}^{\dagger}(t)\mathbf{D}^{\text{4c}}(t)\mathbf{U}(t).$
(71)

The right hand side of Eqs. (68) and (70) involves the matrix product
$\mathbf{\dot{U}}^{\dagger}(t)\mathbf{U}(t)$ which has nonzero off-diagonal
blocks that prevent expressing these equations in the complete decoupled
(block-diagonal) form. However, as discussed in Ref. [28] for the case of a
single electric field $\bm{\mathcal{E}}(t)$, the matrix values of
$\mathbf{\dot{U}}^{\dagger}(t)$ are of the order
$\mathcal{O}(|\bm{\mathcal{E}}|\omega c^{-1})$, and therefore become
negligibly small within a weak-field limit ($|\bm{\mathcal{E}}|\\!\ll\\!1$)
and a dipole approximation ($r\omega c^{-1}\\!\ll\\!1$). As a result, the X2C
transformation matrix remains approximately constant in time, _i.e._
$\mathbf{U}(t)\approx\mathbf{U}$, and Eqs. (68) and (70) reduce to the simple
form

$\begin{split}i\hbar\frac{\partial\bm{\tilde{C}}^{\mathrm{4c}}_{i}(t)}{\partial
t}=\mathbf{\tilde{F}}^{\mathrm{4c}}(t)\bm{\tilde{C}}^{\mathrm{4c}}_{i}(t);\qquad
i\hbar\frac{\partial\mathbf{\tilde{D}}^{\mathrm{4c}}(t)}{\partial
t}=\left[\mathbf{\tilde{F}}^{\mathrm{4c}}(t),\mathbf{\tilde{D}}^{\mathrm{4c}}(t)\right].\end{split}$
(72)

This time-independence of matrix $\mathbf{U}$ is generally denoted as the
_adiabatic X2C transformation_ [26, 28].

The best possible transformation matrix $\mathbf{U}$ can be obtained from a
so-called mmfX2C approach [143]. In this approach, $\mathbf{U}$ is obtained _a
posteriori_ from converged 4c SCF HF/KS solutions (MO coefficients) applying
for instance the one-step X2C transformation of Ilias and Saue [89]. From the
real-time dynamics point of view, these solutions are associated with the
initial simulation time $t_{0}$. An important observation is that at $t_{0}$
the 4c occupied positive-energy MO coefficients $\bm{C}^{\text{4c}}_{i}$ as
well as the 4c density matrix $\mathbf{D}^{\text{4c}}$ can be expressed in
terms of their 2c counterparts,

$\begin{split}\bm{C}^{\text{4c}}_{i}(t_{0})&=\mathbf{U}\bm{\tilde{C}}^{\text{2c}}_{i}(t_{0})\quad\Rightarrow\quad\big{[}C^{\text{4c}}\big{]}^{\text{X}}_{\mu
i}=\sum_{\nu}\big{[}U\big{]}^{\text{XL}}_{\mu\nu}\big{[}\tilde{C}^{\text{2c}}\big{]}_{\nu
i}\\\\[5.69046pt]
\mathbf{D}^{\text{4c}}(t_{0})&=\mathbf{U}\mathbf{\tilde{D}}^{\text{2c}}(t_{0})\mathbf{U}^{\dagger}\quad\Rightarrow\quad\big{[}D^{\text{4c}}\big{]}^{\text{XY}}_{\mu\nu}=\sum_{\kappa\lambda}\big{[}U\big{]}^{\text{XL}}_{\mu\kappa}\big{[}\tilde{D}^{\text{2c}}\big{]}_{\kappa\lambda}\big{[}U^{\dagger}\big{]}^{\text{LY}}_{\lambda\nu}.\end{split}$
(73)

Here, X and Y refer to the large-component (L) and small-component (S) subset
of orthonormal AO basis. Within the adiabatic X2C transformation it is assumed
that the relation (73) remain valid also at an arbitrary future time
$t>t_{0}$, and therefore 4c real-time dynamic results can be obtained just
from the solution of simple 2c EOMs

$i\hbar\frac{\partial\bm{\tilde{C}}^{\mathrm{2c}}_{i}(t)}{\partial
t}=\mathbf{\tilde{F}}^{\mathrm{2c}}(t)\bm{\tilde{C}}^{\mathrm{2c}}_{i}(t);\qquad
i\hbar\frac{\partial\mathbf{\tilde{D}}^{\mathrm{2c}}(t)}{\partial
t}=\left[\mathbf{\tilde{F}}^{\mathrm{2c}}(t),\mathbf{\tilde{D}}^{\mathrm{2c}}(t)\right].$
(74)

However, as shown by Knecht and coworkers for static SCF case [91], the
correctly transformed 2c Fock matrix $\mathbf{\tilde{F}}^{\text{2c}}$ involves
a so-called picture-change transformation of density matrix, overlap
distribution matrix, and one- and two-electron integrals. Repisky and
coworkers extended this observation to the time domain and derive [27, 28]:

$\displaystyle\tilde{F}^{\mathrm{2c}}_{\mu\nu}(t)=\left[\mathbf{U}^{\dagger}\mathbf{F}^{\text{4c}}(t)\mathbf{U}\right]^{\text{LL}}_{\mu\nu}$
$\displaystyle=\tilde{h}^{\mathrm{2c}}_{\mu\nu}+\sum_{\kappa\lambda}\tilde{G}^{\text{2c}}_{\mu\nu,\kappa\lambda}\tilde{D}^{\text{2c}}_{\lambda\kappa}(t,\bm{\mathcal{E}},\bm{\mathcal{F}})$
(75) $\displaystyle+\sum_{u\in 0,x,y,z}\int
v^{xc}_{u}\\!\left[\bm{\tilde{\rho}}^{\text{2c}}(\bm{r},t,\bm{\mathcal{E}},\bm{\mathcal{F}})\right]\tilde{\Omega}_{u,\mu\nu}^{\text{2c}}(\bm{r})\,d^{3}\bm{r}$
$\displaystyle-\sum_{u\in
x,y,z}\tilde{P}^{\mathrm{2c}}_{u,\mu\nu}\mathcal{E}_{u}(t)-\sum_{u\in
x,y,z}\tilde{P}^{\mathrm{2c}}_{u,\mu\nu}\mathcal{F}_{u}(t).$

There are two important points to note here: (i) all transformed quantities
are marked with tilde; (ii) the presence of the picture-change transformed
charge distribution matrix ($\bm{\tilde{\Omega}}^{\text{2c}}$) in both 2e and
xc interaction terms makes the evaluation of
$\mathbf{\tilde{F}}^{\mathrm{2c}}$ computationally more demanding than the
original 4c Fock matrix.

Therefore, it is desirable to seek for an approximation that enables us to
carry out electron dynamics simulations in 2c mode such that they are
computationally efficient and reproduce the reference 4c results as closely as
possible. Keeping this in mind, one can compare Eq. (75) with an approximate
and computationally efficient form of the Fock matrix built with untransformed
(without the tilde) two-electron integrals $\mathbf{G}^{\text{2c}}$ and
overlap distribution matrix $\bm{\Omega}^{\text{2c}}$; that is

$\displaystyle F^{\mathrm{2c}}_{\mu\nu}(t)$
$\displaystyle=\tilde{h}^{\mathrm{2c}}_{\mu\nu}+\sum_{\kappa\lambda}G^{\text{2c}}_{\mu\nu,\kappa\lambda}\tilde{D}^{\text{2c}}_{\lambda\kappa}(t,\bm{\mathcal{E}},\bm{\mathcal{F}})$
(76) $\displaystyle+\sum_{u\in 0,x,y,z}\int
v^{xc}_{u}\\!\left[\bm{\rho}^{\text{2c}}(\bm{r},t,\bm{\mathcal{E}},\bm{\mathcal{F}})\right]\Omega_{u,\mu\nu}^{\text{2c}}(\bm{r})\,d^{3}\bm{r}$
$\displaystyle-\sum_{u\in
x,y,z}\tilde{P}^{\mathrm{2c}}_{u,\mu\nu}\mathcal{E}_{u}(t)-\sum_{u\in
x,y,z}\tilde{P}^{\mathrm{2c}}_{u,\mu\nu}\mathcal{F}_{u}(t).$

Here, it is important to emphasize that $\bm{\rho}^{\text{2c}}$ also remains
untransformed in the sense that an untransformed $\bm{\Omega}_{u}^{\text{2c}}$
is used but with the correctly transformed density matrix
$\mathbf{\tilde{D}}^{\text{2c}}$. We immediately find that the difference
between these two Fock matrices expresses the picture-change corrections
associated with the two-electron integrals and the xc contribution

$\Delta\tilde{F}^{\text{2c}}_{\mu\nu}(t)=\tilde{F}^{\text{2c}}_{\mu\nu}(t)-F^{\text{2c}}_{\mu\nu}(t)=\sum_{\kappa\lambda}\Delta\tilde{G}^{\text{2c}}_{\mu\nu,\kappa\lambda}\tilde{D}^{\text{2c}}_{\lambda\kappa}(t)+\Delta\tilde{F}^{\text{2c,xc}}_{\mu\nu}(t),$
(77)

where

$\begin{split}\Delta\tilde{G}^{\text{2c}}_{\mu\nu,\kappa\lambda}&=\tilde{G}^{\text{2c}}_{\mu\nu,\kappa\lambda}-G^{\text{2c}}_{\mu\nu,\kappa\lambda},\\\\[5.69046pt]
\Delta\tilde{F}^{\text{2c,xc}}_{\mu\nu}(t)&=\int
v_{k}^{xc}\\!\left[\bm{\tilde{\rho}}^{\text{2c}}(\bm{r},t)\right]\tilde{\Omega}_{k,\mu\nu}^{\text{2c}}(\bm{r})\,d^{3}\bm{r}-\int
v_{k}^{xc}\\!\left[\bm{\rho}^{\text{2c}}(\bm{r},t)\right]\Omega_{k,\mu\nu}^{\text{2c}}(\bm{r})\,d^{3}\bm{r}.\end{split}$
(78)

Here, we dropped the dependence on $\bm{\mathcal{E}}$ and $\bm{\mathcal{F}}$
for clarity. The central idea of X2C real-time electron dynamics is the
solution of the 2c LvN equation (74) with the Fock matrix

$\displaystyle\tilde{F}^{\mathrm{2c}}_{\mu\nu}(t)$
$\displaystyle=\tilde{h}^{\mathrm{2c}}_{\mu\nu}+\Delta\tilde{F}^{\text{2c}}_{\mu\nu}(t)+\sum_{\kappa\lambda}G^{\text{2c}}_{\mu\nu,\kappa\lambda}\tilde{D}^{\text{2c}}_{\lambda\kappa}(t,\bm{\mathcal{E}},\bm{\mathcal{F}})$
(79) $\displaystyle+\sum_{u\in 0,x,y,z}\int
v^{xc}_{u}\\!\left[\bm{\rho}^{\text{2c}}(\bm{r},t,\bm{\mathcal{E}},\bm{\mathcal{F}})\right]\Omega_{u,\mu\nu}^{\text{2c}}(\bm{r})\,d^{3}\bm{r}$
$\displaystyle-\sum_{u\in
x,y,z}\tilde{P}^{\mathrm{2c}}_{u,\mu\nu}\mathcal{E}_{u}(t)-\sum_{u\in
x,y,z}\tilde{P}^{\mathrm{2c}}_{u,\mu\nu}\mathcal{F}_{u}(t),$

where $\Delta\mathbf{\tilde{F}}^{\mathrm{2c}}(t)$ accounts for the picture-
change corrections associated with the 2e integrals and the xc contribution.
Note that $\mathbf{\tilde{F}}^{\mathrm{2c}}(t)$ in Eqs. (75) and (79) are
equal, and all differences between various flavours of X2C are due to
approximations in $\Delta\mathbf{\tilde{F}}^{\mathrm{2c}}(t)$.

In the simplest but least accurate case, dubbed one-electron X2C (1eX2C),
$\Delta\mathbf{\tilde{F}}^{\mathrm{2c}}(t)$ in Eq. (79) is completely
discarded, while the decoupling matrix $\mathbf{U}$ is obtained simply from
the parent one-electron Dirac Hamiltonian. This approach was employed in
pioneering X2C RT-TDDFT implementations [26, 90]. Due to its simplicity the
1eX2C Hamiltonian still remains very popular, but caution is needed when
applying this model beyond valence electric properties as shown for instance
in Ref. [28].

In the second model, coined as molecular mean-field X2C (mmfX2C),
$\Delta\mathbf{\tilde{F}}^{\mathrm{2c}}(t)$ in Eq. (79) is approximated by a
_static_ model $\Delta\mathbf{\tilde{F}}^{\mathrm{2c}}$, which is evaluated
according to Eqs. (77) and (78) only _once_ using the converged 4c molecular
self-consistent field solutions [27, 28]. Similarly, $\mathbf{U}$ is
determined from the same 4c solutions. For theoretical and numerical
justification of the static approximation used in the real-time and response
mmfX2C theory, readers are referred to the original publication [28]. Due to
the late X2C transformation (post-SCF), the mmfX2C approach was found as most
accurate among all X2C Hamiltonian models [28, 91]; though, the price for this
accuracy is the implementation and execution of 4c molecular SCFs.

In line with the idea of Knecht _et al._ [91] on the amfX2C Hamiltonian for
time-independent Hartree-Fock and Kohn-Sham mean-field theories, one may
exploit the local atomic nature of the static picture-change correction matrix
$\Delta\mathbf{\tilde{F}}^{\text{2c}}$ discussed in the previous paragraph. In
the third model, dubbed as atomic mean-field X2C (amfX2C) [27, 28, 91],
$\Delta\mathbf{\tilde{F}}^{\mathrm{2c}}(t)$ in Eq. (79) is approximated by a
_static_ model $\Delta\mathbf{\tilde{F}}^{\mathrm{amfX2C}}_{\bigoplus}$
obtained by a superposition of converged atomic quantities rather than the
converged molecular one, _i.e._

$\displaystyle\Delta\mathbf{\tilde{F}}^{\text{2c}}(t)\approx\Delta\mathbf{\tilde{F}}^{\text{amfX2C}}_{\bigoplus}=\bigoplus_{K=1}^{M}\Delta\mathbf{\tilde{F}}^{\text{2c}}_{K}[\mathbf{\tilde{D}}_{K}^{\text{2c}}].$
(80)

Here, $K$ runs over all atoms in an $M$-atomic system.

The main advantage of the amfX2C approach is that it introduces picture-change
corrections to both spin-independent and spin-dependent parts of the two-
electron and xc interaction just from simple atomic quantities. On the other
hand, the fact that $\Delta\mathbf{\tilde{F}}^{\text{2c}}_{\bigoplus}$ has
only atomic diagonal blocks means that, for instance, the off-diagonal
electron-nucleus contribution will not cancel out with the direct electron-
electron contribution at long distances from the atomic centers. This becomes
problematic in solid-state calculations, where the exact cancellation of these
contributions is essential at long distances. In fact, this motivated Knecht
and coworkers [91] to introduce our last X2C Hamiltonian model, called
extended amfX2C (eamfX2C). The generalization of eamfX2C to the time domain
was recently discussed by Konecny and coworkers [28], and it requires to
approximate $\Delta\mathbf{\tilde{F}}^{\mathrm{2c}}(t)$ in Eq. (79) by a
_static_ model $\Delta\mathbf{\tilde{F}}^{\mathrm{eamfX2C}}_{\bigoplus}$
obtained from the time-independent version of equation (77)

$\Delta\tilde{F}^{\text{eamfX2C}}_{\bigoplus,\mu\nu}=\tilde{F}^{\text{2c}}_{\bigoplus,\mu\nu}-F^{\text{2c}}_{\bigoplus,\mu\nu}=\sum_{\kappa\lambda}\Delta\tilde{G}^{\text{2c}}_{\mu\nu,\kappa\lambda}\tilde{D}^{\text{2c}}_{\bigoplus,\lambda\kappa}+\Delta\tilde{F}^{\text{2c,xc}}_{\bigoplus,\mu\nu},$
(81)

with elements on the right hand side given in Eq. (78). The picture-change
corrections associated with the two-electron integrals and the xc contribution
involve the 2c density matrix $\mathbf{\tilde{D}}^{\text{2c}}_{\bigoplus}$
obtained from a superposition of converged 4c atomic density matrices
$\mathbf{D}^{\text{4c}}_{K}$, i.e.

$\displaystyle\mathbf{\tilde{D}}^{\text{2c}}_{\bigoplus}=\bigoplus_{K=1}^{M}\left[\mathbf{U}_{K}^{\dagger}\mathbf{D}^{\text{4c}}_{K}\mathbf{U}_{K}\right]^{\text{LL}}.$
(82)

Here, $K$ runs over all atoms in an $M$-atomic system.

##  6 Real-time propagation

###  6.1 Evolution operator

The solution of the TDHF and TDKS equations for an arbitrary time $t$ can be
written in the compact form by defining the _evolution operator_
$U(t,t^{\prime})$ that propagates the state from time $t^{\prime}$ to time $t$
as

$\varphi_{i}(\bm{r},t)=U(t,t^{\prime})\varphi_{i}(\bm{r},t^{\prime}),$ (83)

or, equivalently, in the matrix form

$\mathbf{C}(t)=\mathbf{U}(t,t^{\prime})\mathbf{C}(t^{\prime}).$ (84)

We can also use the same evolution operator to obtain the solution of the LvN
equation in the language of the RDM as

$\mathbf{D}(t)=\mathbf{U}(t,t^{\prime})\mathbf{D}(t^{\prime})\mathbf{U}^{\dagger}(t,t^{\prime}).$
(85)

It is required that the evolution operator is unitary, _i.e._
$\mathbf{U}^{\dagger}(t,t^{\prime})\mathbf{U}(t,t^{\prime})=\mathbb{I}$, so
that the time evolution preserves the norm of the $\varphi_{i}$ as well as the
idempotence and trace of the density matrix. It follows from the definition
that

$\displaystyle\mathbf{U}(t,t)$ $\displaystyle=\mathbb{I},$ (86)
$\displaystyle\mathbf{U}(t_{3},t_{1})$
$\displaystyle=\mathbf{U}(t_{3},t_{2})\mathbf{U}(t_{2},t_{1}),$ (87)
$\displaystyle\mathbf{U}^{-1}(t_{1},t_{2})$
$\displaystyle=\mathbf{U}(t_{2},t_{1}).$ (88)

The last property is related to the TRS described in Section 4.3 and only
holds when no external magnetic fields are present.

Inserting the definition of $\mathbf{U}$ in Eq. (85) into the LvN equation
recasts the problem of time propagation into determining $\mathbf{U}$ by
solving

$i\hbar\frac{\partial}{\partial
t}\mathbf{U}(t,t^{\prime})=\mathbf{F}(t)\mathbf{U}(t,t^{\prime}).$ (89)

It is possible to write a closed-form solution of this equation in the form of
the Dyson series as

$\mathbf{U}(t,t^{\prime})=\sum_{n=0}^{\infty}\frac{(-i/\hbar)^{n}}{n!}\int_{t^{\prime}}^{t}dt_{1}\ldots\int_{t^{\prime}}^{t}dt_{n}\mathcal{T}\big{\\{}\mathbf{F}(t_{1})\ldots\mathbf{F}(t_{n})\big{\\}},$
(90)

where $\mathcal{T}$ represents the time-ordering of the product such, that the
leftmost term has the latest time, and each following term is applied at an
earlier time than the one before it. The time-ordering is necessary since
$\left[\mathbf{F}(t_{1}),\mathbf{F}(t_{2})\right]\neq 0$. Except for the time-
ordering, this series represents the expansion of the exponential function and
thus is often written in the short-hand form

$\mathbf{U}(t,t^{\prime})=\mathcal{T}\exp\left[-\frac{i}{\hbar}\int_{t^{\prime}}^{t}\mathbf{F}(\tau)d\tau\right].$
(91)

This expression, albeit formally exact, requires truncation of the series in
numerical implementations. Such truncation inevitably leads to the loss of the
unitary property of $\mathbf{U}$, and consequently the idempotence and trace
of the density matrix, which can result in numerically unstable time
propagation [144, 145, 146].

###  6.2 Magnus expansion

As an alternative to the Dyson expansion, the evolution operator can be
written as a true exponential function that does not require the time
ordering, _i.e._ in the form of the exponent of the infinite series as

$\mathbf{U}(t,t^{\prime})=e^{\mathbf{A}(t,t^{\prime})},$ (92)

where

$\mathbf{A}(t,t^{\prime})=\sum_{n=1}^{\infty}\mathbf{A}_{n}(t,t^{\prime}).$
(93)

This form was proposed by Magnus in 1954 with the first terms given by [147,
117]

$\displaystyle\mathbf{A}_{1}(t,t^{\prime})=$
$\displaystyle\frac{1}{i\hbar}\int_{t^{\prime}}^{t}dt_{1}\mathbf{F}(t_{1}),$
(94) $\displaystyle\mathbf{A}_{2}(t,t^{\prime})=$
$\displaystyle-\frac{1}{2}\left(\frac{1}{i\hbar}\right)^{2}\int_{t^{\prime}}^{t}dt_{2}\int_{t^{\prime}}^{t_{2}}dt_{1}\big{[}\mathbf{F}(t_{1}),\mathbf{F}(t_{2})\big{]},$
(95) $\displaystyle\mathbf{A}_{3}(t,t^{\prime})=$
$\displaystyle-\frac{1}{6}\left(\frac{1}{i\hbar}\right)^{3}\int_{t^{\prime}}^{t}dt_{3}\int_{t^{\prime}}^{t_{3}}dt_{2}\int_{t^{\prime}}^{t_{2}}dt_{1}\left(\Big{[}\mathbf{F}(t_{1}),\big{[}\mathbf{F}(t_{2}),\mathbf{F}(t_{3})\big{]}\Big{]}\right.$
$\displaystyle+\left.\Big{[}\big{[}\mathbf{F}(t_{1}),\mathbf{F}(t_{2})\big{]},\mathbf{F}(t_{3})\Big{]}\right).$
(96)

Both Dyson and Magnus expansions are equivalent if their respective series are
considered in the infinite limit, however, the Magnus expansion has an
advantage in the numerical implementations as it retains the unitary property
of the evolution operator even if truncated after any number of terms. Higher-
order propagators based on the Magnus expansion require the evaluation of
multiple commutators of the Hamiltonian (Fock) matrix, _e.g._ in Eqs. (95) and
(96). A commutator-free version based on the Magnus series can be obtained by
writing the evolution operator as an (infinite) product of exponentials [148,
149, 150] which leads to a powerful approach for deriving higher order
commutator-free exponential time propagators.

###  6.3 Approximate evolution

Real-time simulations typically start from an initial state defined at $t=0$
by $\mathbf{D}(t\\!=\\!0)$. This initial state is in most cases obtained from
a converged ground-state optimization procedure, though starting the evolution
from approximate excited states is also possible [151]. Once the evolution
operator is determined, the orbitals and RDM at arbitrary times can be
calculated by applying Eqs. (83) and (85). For instance, $\mathbf{D}(t)$ is
obtained as follows

$\mathbf{D}(t)=\mathbf{U}(t,0)\mathbf{D}(0)\mathbf{U}^{\dagger}(t,0).$ (97)

In practice, however, the global propagator $\mathbf{U}(t,0)$ is not known –
this is the case even if the infinite series (Dyson or Magnus) from the
previous section are truncated. Numerical implementations require that the
time is discretized into a finite number of time steps $N$ of size $\Delta t$,
and the evolution operator is factored using Eq. (87) as

$\mathbf{U}(N\Delta t,0)=\prod_{i=0}^{N-1}\mathbf{U}\big{(}(i+1)\Delta
t,i\Delta t\big{)}.$ (98)

Thus, the propagation over one time step is achieved by the application of
$\mathbf{U}(t+\Delta t,t)$ on the density matrix or orbitals, where $t\equiv
i\Delta t$. This short-time propagation allows us to approximate the integrals
in $\mathbf{U}(t+\Delta t,t)$

The most commonly used approximation for $\mathbf{U}$ is the _midpoint Magnus_
propagator [144, 152, 153, 154]

$\mathbf{U}(t+\Delta
t,t)\approx\exp\left[\frac{1}{i\hbar}\mathbf{F}\left(t+\frac{\Delta
t}{2}\right)\Delta t\right],$ (99)

where only the first term $\mathbf{A}_{1}$ in Eq. (94) of the Magnus expansion
is considered, and the integral $\int_{t}^{t+\Delta t}dt$ is approximated
using the midpoint quadrature. This integrator is of the second order as it is
correct to $\mathcal{O}(\Delta t^{2})$. Eq. (99) is also the basis of the
modified midpoint unitary transformation (MMUT) method in the literature [152,
155], where the time step is modified to $2\Delta t$. A fourth-order Magnus
propagator can be constructed by taking the first two terms in the Magnus
series $\bm{A}_{1}$ and $\bm{A}_{2}$ and approximating the integrals using a
two-point Gaussian quadrature [144, 93]. A comparison of the performance of
the second- and fourth-order Magnus propagators with the predictor–corrector
scheme was presented in the recent study of Müller, Sharma, and Sierka [156]
based on their efficient implementation of RT-TDDFT. The matrix exponential in
Eq. (99) can be evaluated directly by diagonalizing the Fock matrix and
constructing the exponent from its eigenvalues. However, techniques that
circumvent the expensive diagonalization step by using a series of matrix
multiplications, such as the Baker-Campbell-Hausdorff formula [22] or the
Chebyshev expansion [157], lower the computational cost as well as improve
parallelization.

Integrators that are not based on the Magnus series, such as the Runge-Kutta
Method [158, 153] or the Crank–Nicholson propagator [159, 73, 160, 161, 20]

$\mathbf{U}(t+\Delta t,t)\approx\frac{1+i\mathbf{F}\left(t+\frac{\Delta
t}{2}\right)\frac{\Delta t}{2\hbar}}{1-i\mathbf{F}\left(t+\frac{\Delta
t}{2}\right)\frac{\Delta t}{2\hbar}},$ (100)

can also be used. However, unlike the midpoint Magnus or the Crank–Nicholson,
the Runge-Kutta integrator does not preserve the unitary property of
$\mathbf{U}$ and can lead to instabilities in the time evolution, since
neither the electron number nor the total energy are strictly conserved during
a time propagation that is not unitary [144, 162]. On another note, the
exponential Runge-Kutta method [163] was successfully used to solve equations
of motion in the time-dependent coupled-cluster theory [164] which is based on
the exponential coupled-cluster parametrization of the wave function [165,
166].

The issue of stable time propagation requires even more attention in the TDHF
and TDKS theories. Due to the mean-field, XC, and HF exchange terms, the Fock
matrix in the LvN equation depends on the density matrix $\mathbf{D}(t)$ that
is not known at the time $t$. Hence, the LvN equation with this Fock matrix is
nonlinear and must be solved self-consistently, _i.e._ the Fock matrix is
constructed from the density matrix from the previous iterations (referred to
as microiterations in the context of time domain methodology). Unless the loop
over microiterations is introduced, this implicit time-dependence of the Fock
matrix on the unknown density matrix results in the inability to express the
future time midpoint Fock matrix $\mathbf{F}\left(t+\frac{\Delta t}{2}\right)$
that appears in most approximations for the evolution operator. This issue is
commonly mitigated by using predictor–corrector or extrapolation–interpolation
schemes [77, 22, 24, 160, 20], where the unknown midpoint Fock matrix is first
constructed from the previous time step using the linear extrapolation

$\mathbf{F}\left(t+\frac{\Delta
t}{2}\right)=2\mathbf{F}(t)-\mathbf{F}\left(t-\frac{\Delta t}{2}\right).$
(101)

Once the $\mathbf{D}(t+\Delta t)$ is obtained by applying the evolution
operator $\mathbf{U}(t+\Delta t,t)$, a new Fock matrix at $t+\Delta t$ is
formed. The two Fock matrices at $t$ and $t+\Delta t$ are then linearly
interpolated to create the updated midpoint Fock matrix

$\mathbf{F}\left(t+\frac{\Delta
t}{2}\right)=\frac{1}{2}\mathbf{F}(t)+\frac{1}{2}\mathbf{F}\left(t+\Delta
t\right),$ (102)

and the time propagation restarts from the initial time $t$. This process is
repeated until the self-consistence is reached. A thorough comparison of
various propagation schemes in the context of nonrelativistic TDKS can be
found in the study of Pueyo _et al._ [167]. A detailed analysis is provided in
the chapter by Ye _et al._ [93]. Simulations of X-ray absorption spectra are
particularly sensitive to the size of the time step as the core excitations
appearing in the high energy region typically involve rapid oscillations of
the wave functions that require a very small time step to be described
properly [25, 94, 168, 93]. For such studies, Ye _et al._ [93] presented a
relativistic X2C approach based on the fourth-order commutator-free Magnus
propagator that adaptively chooses optimal time step and simulation time. The
discussion of the approximate time evolution is expanded further in this
volume in the chapter by L. Ye, H. Wang, Y. Zhang, Y. Xiao and W. Liu entitled
“Real-Time Time-Dependent Density Functional Theories with Large Time Step and
Short Simulation Time”.

##  7 Signal processing

RT-TDSCF solves the EOM by direct propagation in the time domain. This allows
direct simulation of time-resolved experiments and to obtain the entire
spectral information from a single real-time calculation. On the other hand,
physical quantities of experimental interest are often defined in the
frequency domain, which creates a demand for techniques that efficiently
extract the frequency-domain quantity from the simulations that are carried
out with finite numerical accuracy, time step length, and simulation time.

A time-dependent property $f(t)$ can be translated to the frequency domain as
the Fourier integral

$\tilde{f}(\omega)=\int^{\infty}_{-\infty}f(t)e^{i\omega t}dt.$ (103)

However, due to the presence of periodic oscillations in $f(t)$ that originate
in the quantum mechanical evolution containing the excitation energies, the
integral in Eq. (103) leads to $\delta$-functions in the frequency domain.
Such stick spectra are difficult to describe in numerical simulations, hence,
the Fourier integral is replaced with the Laplace transform

$\tilde{f}(\omega)=\int^{\infty}_{0}f(t)e^{i\omega t-\gamma t}dt,$ (104)

where we assumed that there is no response of the studied system before the
perturbation is applied at time $t=0$, and we introduced a phenomenological
damping parameter $\gamma>0$. This damping parameter accounts for the fact
that practical real-time simulations are performed with a finite simulation
time, and truncating the periodically oscillating signal that is not damped
results in undesirable features in the spectrum. To this end, $\gamma$ is
empirically set to a value inversely proportional to the simulation time to
ensure that the oscillations diminish before the simulation is terminated.
Analytic evaluation of Eq. (104) for periodic signals leads to a series of
Lorentzian peaks with finite width in the spectrum.

Signal processing techniques serve as a means to extract the frequency-
dependent molecular properties $\tilde{f}(\omega)$ from the simulation results
$f(t)$, _i.e._ to approximate the integrals in Eqs. (103) and (104). In this
section we present an overview of approaches used in the context of RT-TDSCF,
starting with the simplest one, the Discrete Fourier transform, that is
instructive in order to explain the relationship between the time and
frequency domains, and then moving to the more sophisticated methods that
achieve better resolution from a shorter simulation. In general, spectra with
a high density of states pose a bigger challenge to the signal processing
methods.

###  7.1 Discrete Fourier transform

In practice, RT-TDSCF simulations are performed in a series of discrete time
steps $t_{j}$ for which the induced dipole moment is calculated from a trace
of the dipole moment matrix and the time-dependent density matrix

$\bm{\mu}^{\mathrm{ind}}(t_{j})=\operatorname{Tr}[\mathbf{P}\mathbf{D}(t_{j})]-\bm{\mu}^{\mathrm{static}},$
(105)

where the static dipole moment is calculated as
$\bm{\mu}^{\mathrm{static}}=\operatorname{Tr}[\mathbf{P}\mathbf{D}_{0}]$. The
most straightforward approximation to Eq. (104) is the discrete Fourier
transform

$\tilde{f}_{k}=\sum_{j=0}^{n-1}\Delta t\,f_{j}e^{2\pi i\frac{jk}{n}-\gamma
j\Delta t}.$ (106)

Here, $k=0,1,\ldots,n-1$ where $n$ is the number of time steps and
$\omega_{k}=2\pi k/(n\Delta t)$ is the $k$-th frequency point. The
coefficients $f_{j}\equiv\bm{\mu}^{\mathrm{ind}}(t_{j})$ and
$\tilde{f}_{k}\equiv\bm{\mu}^{\mathrm{ind}}(\omega_{k})$ represent the
components of the induced dipole moment in time and frequency domains,
respectively.

###  7.2 Relationship between the time and frequency domains

The frequency-domain results are obtained by discrete Fourier transform of the
time-domain results. If we perform a time-domain simulation that consists of
$n$ steps of length $\Delta t$, the Fourier transform yields a frequency-
domain interval of length

$\Omega=\frac{2\pi}{\Delta t}.$ (107)

Since the number of points in both domains is the same, the resolution in the
frequency domain is

$\Delta\omega=\frac{2\pi}{n\Delta t}.$ (108)

This relationship tells us that in order to increase the resolution in
calculated spectra we need to increase the total simulation length $n\Delta t$
by increasing the number of time steps (which makes the simulation more time
consuming) or increasing the size of the time step (which puts extra demands
on the solver). However, because the frequency-domain interval depends
inversely on the time-step length, see Eq. (107), in order to describe high-
frequency parts of spectra, such as in X-ray spectroscopies, shorter time
steps are required. Therefore, a balance between the resolution, frequency
range and computational cost has to be achieved by choosing suitable
simulation parameters. Eq. (108) represents the major limitation for obtaining
spectra with high resolution when using the discrete Fourier transform.

###  7.3 Padé approximants

The use of the Padé approximation as a signal processing technique [169, 170]
was introduced to RT-TDSCF by Bruner, LaMaster, and Lopata [171] and quickly
gained popularity due to its advantage compared to the widespread Fourier
transform. In the Padé approximants, the expression for the Fourier components
$\tilde{f}(\omega)$ (_e.g._ induced dipole moment in frequency domain)

$\tilde{f}(\omega)=\sum_{j=0}^{M}f(t_{j})\Delta te^{i\omega j\Delta
t}e^{-\gamma j\Delta t},$ (109)

is understood as a power series

$\tilde{f}(\omega)=\sum_{j=0}^{M}c_{j}(z_{\omega})^{j},$ (110)

where $z_{\omega}=e^{i\omega\Delta t}$ and $c_{j}=f(t_{j})\Delta te^{-\gamma
j\Delta t}$ that can be approximated as a division of two other power series
using the Padé ansatz

$\tilde{f}(z)=\frac{\sum_{k=0}^{N}a_{k}z^{k}}{\sum_{k^{\prime}=0}^{N}b_{k^{\prime}}z^{k^{\prime}}},$
(111)

where $N=M/2$. The comparison of Eqs. (110) and (111) leads to a system of
equations, or a matrix equation, for the coefficients $b$

$\mathbf{b}=\mathbf{G}^{-1}\mathbf{d},$ (112)

where $G_{km}=c_{N-m+k}$ and $d_{k}=-c_{N+k}$. The system is overdetermined,
leading to a customary choice of setting $b_{0}=1$. The knowledge of the
$b$-coefficients can then be used to determine the $a$-coefficients from
$a_{k}=\sum_{m=0}^{k}b_{m}c_{k-m}$. The $a$ and $b$ coefficients are
subsequently used to approximate the Fourier transform $\tilde{f}(\omega)$
with the advantage that $a$ and $b$ are not functions of frequency. Hence, the
frequency can be chosen at will without the limitation of Eq. (108), which
allows the spectrum to be evaluated with arbitrary frequency resolution. In
practice, the Padé approximation can suffer from numerical instabilities. To
mitigate this, the original work [171] suggested to combine it with a MO-based
decomposition [24], where each occupied–virtual MO pair that contributes to
the net frequency-dependent polarizability is transformed into the frequency
domain using its own Padé approximation, _i.e._ the coefficients $a_{k}$ and
$b_{k}$ become different for each MO pair. A spectrum constructed from an
individual MO pair is typically sparser and thus less prone to defective
behaviour. The final spectrum is then a sum of spectra over all MO pairs.

###  7.4 Compressed sensing

Compressed sensing is a technique based on the observation that a small number
of points in time domain is sufficient to sample a frequency domain signal
that is sparse, _i.e._ when many Fourier-domain coefficients $\tilde{f}_{k}$
in Eq. (106) are near zero. The application of compressed sensing in RT-TDDFT
was first explored by Andrade, Sanders, and Aspuru-Guzik [172] in the context
of electronic and nuclear dynamics for the calculation of vibrational and
optical spectra. The compressed sensing method recasts the problem of finding
the Fourier coefficients into solving a system of linear equations

$\mathbf{A}\tilde{\bm{f}}=\bm{f},$ (113)

where $f_{j}\equiv(\bm{f})_{j}$ and $\tilde{f}_{k}\equiv(\tilde{\bm{f}})_{k}$
are vectors with components in time and frequency domains, respectively, and
$\mathbf{A}$ is the matrix containing the complex phase factors. This system
allows for a different number of time and frequency points
$j=0,1,\ldots,n_{t}-1$ and $k=0,1,\ldots,n_{\omega}-1$. For a small number of
time points $n_{t}<n_{\omega}$ the system is underdetermined with infinitely
many solutions $\tilde{\bm{f}}$. The sparse solution with the largest number
of zero coefficients is then obtained by finding $\tilde{\bm{f}}$ that
minimize the norm $|\tilde{f}|$ while satisfying

$\left|\mathbf{A}\tilde{\bm{f}}-\bm{f}\right|<\eta,$ (114)

where $\eta\ll 1$ accounts for a certain amount of numerical noise in the
signal. Even though the benefits of compressed sensing are expected to be
lower for dense signals, for systems with low density-of-states, significant
savings can be obtained by reconstructing the spectra from shorter time
simulations [172].

###  7.5 Filter diagonalization

Time signals often take the form of a sum of damped oscillations, _i.e._

$f(t)=\sum_{m}d_{m}e^{-i\omega_{m}t},$ (115)

where the frequencies $\omega_{m}$ can be considered complex to account for
the damping factors. In an ideal case (with no numerical noise), this is also
the form obtained from the quantum mechanical time propagation discussed
Section 6, with $\omega_{m}$ and $d_{m}$ corresponding to excitation energies
and oscillator strengths, respectively. Determining the values of $\omega_{m}$
and $d_{m}$ from the known signal $f(t)$ is referred to as the _harmonic
inversion_ problem [173], and it was the connection to the quantum mechanical
evolution that lead to the formulation of harmonic inversion as an eigenvalue
problem [174, 175, 176]. The signal $f(t)$ in Eq. (115) can be considered a
correlation function

$f(t)=\braket{\psi_{0}}{e^{-i\hat{\Omega}t}}{\psi_{0}},$ (116)

some unknown Hamiltonian $\hat{\Omega}$ and an initial state $\psi_{0}$. The
frequencies $\omega_{m}$ are the eigenvalues of $\hat{\Omega}$ and are
obtained by solving the generalized eigenvalue equation

$\mathbf{U}\bm{b}_{m}=u_{m}\mathbf{S}\bm{b}_{m},$ (117)

where $\hat{U}=e^{-i\hat{\Omega}\Delta t}$, $u_{m}=e^{-i\omega_{m}\Delta t}$,
and $\Delta t$ is the time step. This equation can be formulated in Krylov
basis constructed by a consecutive application of the operator $\hat{U}$ on a
reference vector $\bm{v}_{0}$ as $\bm{v}_{j}:=\hat{U}^{j}\bm{v}_{0}$. In such
a basis, the matrices $\mathbf{U}$ and $\mathbf{S}$ take simple forms of
$U_{jj^{\prime}}=f_{j+j^{\prime}+1}$ and $S_{jj^{\prime}}=f_{j+j^{\prime}}$,
respectively, obtained from the time signal as $f_{j}\equiv f(j\Delta t)$. The
coefficients $d_{m}$ are calculated using the eigenvectors $\bm{b}_{m}$ as
$\sqrt{d_{m}}=\bm{b}_{m}^{T}\bm{f}$. Note, that this method assumes that the
time signal has the form of Eq. (115), but it is not necessary that the signal
was generated by an actual propagation with a quantum mechanical Hamiltonian,
nor do we need to know the explicit forms of $\hat{\Omega}$ or $\psi_{0}$.

Even though Eq. (117) provides a formally exact solution to the harmonic
inversion problem (up to the dimension of the Krylov vector space), it suffers
from the cubic scaling $\mathcal{O}(n^{3})$ of the diagonalization procedure
with the number of time steps compared to $\mathcal{O}(n\log n)$ scaling of
the Discrete Fourier Transform. The filter diagonalization method circumvents
this problem by transforming the matrices in Eq. (117) from the Krylov basis
$\bm{v}_{j}$ into the Fourier basis

$\bm{w}_{k}=\sum_{j=0}^{n-1}e^{ij\Delta t\xi_{k}}\bm{v}_{j},$ (118)

where the frequencies $\xi_{k}$ can be chosen to form an equidistant grid.
Since the basis vectors $\bm{w}_{k}$ are localized in the frequency domain,
the eigenvectors of the operator $\hat{U}$ can be expressed using a small
number $k=1,\ldots,n_{\text{win}}\ll n$ of $\bm{w}_{k}$. Nonnegligible
contributions to the $m$-th eigenvector arise only from basis vectors for
which $\xi_{k}\approx\omega_{m}$. Thus, matrices $\mathbf{U}$ and $\mathbf{S}$
expressed in this basis exhibit large diagonal and diminishing off-diagonal
terms. This enables defining a small spectral window
$\left[\omega_{\text{min}},\omega_{\text{max}}\right]$ of $n_{\text{win}}$
frequencies and diagonalizing the matrix $\mathbf{U}$ in Eq. (117) in the
Fourier subspace spanned by $\bm{w}_{k}$. The main disadvantage of the filter
diagonalization method is the assumption that the time signal takes the form
of Eq. (115). Even though the method can efficiently circumvent the
uncertainty relation in Eq. (108) and provide high resolution spectra of
sparsely distributed dominant peaks, practical real-time simulations are
hampered by numerical noise, which complicates the use of the filter
diagonalization method for extracting spectral information in highly dense
regions.

##  8 Calculation of molecular properties using real-time methods

In this section we review some of the areas where relativistic RT-TDSCF
methods have been employed to combine the advantages of both the relativistic
and real-time treatments. Moreover, our aim is to explain the physical context
of the molecular properties, to show how to construct a computational protocol
for obtaining these properties, and what aspects of the calculations to pay
attention to. First we explore the calculation of linear and non-linear
response properties where real-time propagation is an alternative to
perturbation theory. Then we focus on non-equlibrium spectroscopies, for which
real-time methods are the only viable approach. We focus on heavy-element
systems and X-ray spectroscopies, where relativistic effects are paramount.
Such applications are the primary motivation behind the development of
relativistic real-time methods.

###  8.1 Linear response properties

Linear response properties include some of the most commonly measured
spectroscopies such as electron absorption spectroscopy (EAS), including X-ray
absorption spectroscopy (XAS), and chiroptical spectroscopies such as electron
circular dichroism (ECD) and optical rotatory dispersion (ORD). [177]
Therefore, they are both experimentally relevant as well as provide a good
introduction to the workflow and analysis of real-time simulations.

In this section, we consider a molecular system perturbed by a single external
field $\bm{\mathcal{E}}(t)$ interacting with the molecule within a dipole
approximation which induces a time-dependent dipole response in the molecule.
The response is a physical quantity $R(t)$ calculated as an expectation value
of its operator $\hat{R}$ from the time-dependent wave function,
$R(t)=\braket{\Psi(t)}{\hat{R}}{\Psi(t)}$. In cases when $\hat{R}$ is a one-
electron operator, such as the electric dipole operator, its expectation value
can be calculated as the trace that contains the time-dependent one-electron
RDM $R(t)=\operatorname{Tr}[\mathbf{D}(t)\mathbf{R}]$. Restricting ourselves
to electric and magnetic dipole moment operators for the perturbation and
response, different combinations lead to different linear spectroscopies.
First, we follow the induced electric dipole resulting from an electric dipole
perturbation, leading to EAS spectrum and the frequency-dependent index of
refraction.

#### Electron absorption spectroscopy

Electron absorption spectrum at all frequencies from UV/Vis to X-ray is
determined by the complex frequency-dependent polarizability tensor
$\bm{\alpha}(\omega)$. The polarizability tensor connects the induced electric
dipole moment to an applied electric field, which in the frequency domain
reads

$\mu^{\mathrm{ind}}_{u}(\omega)=\alpha_{uv}(\omega)\mathcal{E}_{v}(\omega)+\ldots,$
(119)

where $\mathcal{E}_{v}(\omega)$ is the Fourier transform of the external
electric field $\bm{\mathcal{E}}(t)=\mathcal{E}\bm{n}F(t)$ defined by its
amplitude $\mathcal{E}$, directional unit vector $\bm{n}$, and time dependence
$F(t)$. The external field couples to the molecular system via the electric
dipole operator, resulting in its appearance in the Fock matrix as the term

$V^{\mathrm{ext}}(t)=-\mathcal{E}F(t)\bm{n}\cdot\mathbf{P}$ (120)

where $\mathbf{P}$ is a matrix representation of the electric dipole moment
operator. By connecting the molecular induced dipole moment, i.e. the
polarization of a bulk material, to the applied external field, the
polarizability tensor determines the complex index of refraction whose real
part is the standard index of refraction while the imaginary part corresponds
to the attenuation coefficient describing the absorption of light and
appearing in the Beer–Lambert law. However, the more common way of expressing
the absorption spectrum is via the photoabsorption cross-section tensor

$\bm{\sigma}(\omega)=\frac{4\pi\omega}{c}\Im\left[\bm{\alpha}(\omega)\right],$
(121)

where $\Im$ denotes the imaginary part, and $c$ is the speed of light. The
absorption spectrum is then the dipole strength function obtained from the
rotational average of the tensor $\bm{\sigma}$

$S(\omega)=\frac{1}{3}\left[\operatorname{Tr}\bm{\sigma}(\omega)\right],$
(122)

where Tr is the trace over the Cartesian components.

A calculation of the absorption spectrum defined in Eq. 122 from a real-time
simulation then proceeds in the following steps:

  1. 1.

Obtain the reference ground-state density matrix $\mathbf{D}_{0}$ by solving
the time-independent SCF equation.

  2. 2.

Perturb the ground state to obtain the initial state $\mathbf{D}(t_{0})$. This
is usually performed by a short “kick” in the time domain that corresponds to
a broadband pulse in the frequency domain, thus exciting all molecular
transitions. A pure form of such a pulse is the Dirac $\delta$ function
$\bm{\mathcal{E}}(t)=\mathcal{E}\bm{n}\delta(t-t_{0})$ which in practical
simulations can be represented numerically by a narrow Gaussian function or
rectangle, or by an analytic expression

$\mathbf{D}(t_{0})=e^{i\mathcal{E}\bm{n}\cdot\mathbf{P}/\hbar}\mathbf{D}_{0}e^{-i\mathcal{E}\bm{n}\cdot\mathbf{P}/\hbar}.$
(123)

which represents an infinitesimally short time evolution by
$\mathbf{U}(t_{0}+\varepsilon,t_{0}-\varepsilon)$ driven by the
$\delta(t-t_{0})$ field in the limit $\varepsilon\rightarrow 0$ [24].

  3. 3.

Propagate the density matrix $\mathbf{D}(t_{0})$ in time for $n$ time steps of
length $\Delta t$ while recording the induced dipole moment
$\bm{\mu}(t)=\operatorname{Tr}[\mathbf{D}(t)\mathbf{P}]$ at each time step.

  4. 4.

Transform the induced dipole moment to the frequency domain using some of the
techniques discussed in Section 7, i.e. calculate

$\bm{\mu}(\omega)=\int_{t_{0}}^{\infty}dt\bm{\mu}(t)e^{i\omega t-\gamma t},$
(124)

where the damping term $e^{-\gamma t}$ is introduced to resolve the problem
that arises when periodic signals are truncated in numerical simulations with
finite time length.

A graphical summary of these steps is shown in Figure 1 for the case of the
EAS spectrum of the mercury atom (SVWN5 functional [178, 179], uncontracted
Dyall’s VDZ basis [180]) calculated from a four-component RT-TDDFT simulation.
Besides illustrating the workflow of an EAS calculation from RT-TDDFT, the
figure also demonstrates the importance of including relativistic effects in
such simulations by capturing the formally forbidden singlet–triplet
transition. Even though non-relativistic approaches to these transitions have
been presented employing for example a spin-dependent perturbation, [80] it is
only in relativistic theories that include spin–orbit coupling that
singlet–triplet transitions appear in the spectra naturally from first
principles and with correct intensities. Therefore, several works on
relativistic RT-TDSCF for electron absorption spectroscopy have focused on
describing singlet–triplet transitions in the spectra. [24, 90, 83]

Figure 1: Calculated absorption spectrum of the Hg atom using the relativistic
four-component RT-TDDFT. The electronic ground state is perturbed by a
$\delta$-type pulse, the induced dipole moment is recorded in time and
transformed to the frequency domain to yield the spectrum. A formally
forbidden singlet–triplet transition is present due to the relativistic level
of theory. Reprinted with permission from Ref. [24]. Copyright 2015, American
Chemical Society.

#### X-ray absorption

X-ray absorption spectroscopy (XAS) is a subset of electron absorption
spectroscopy where high-frequency X-ray radiation is absorbed in molecules
while exciting electrons from core orbitals. Therefore, in XAS the same
physical quantities are evaluated. However, relativistic effects, both scalar
relativistic effects manifesting as shifts of spectral lines, as well as
spin–orbit interaction causing the splitting of spectral lines, are more
pronounced in XAS necessitating the use of computational methods based on
relativistic Hamiltonians. These effects are observable even in light (3rd
row) elements [25], highlighting the need for a relativistic description also
in these cases. The computational protocol for calculating XAS is the same as
presented in the previous paragraph for EAS in the UV/Vis frequency range with
two important caveats: (i) X-ray absorption occurs at higher frequencies so
that the settings of the simulation such as time step and the number of time
steps have to be adjusted in order to reach the desired frequencies with
sufficient resolution and numerical accuracy; (ii) in simulations using finite
atom-centered basis sets, a broadband $\delta$-type pulse excites all
molecular modes including excitations from valence orbitals to high-lying
above-ionization virtuals that may fall into the XAS frequency range, but are
non-physical relicts of an improper description of continuum states, and thus
have to be eliminated either in post-processing or during the application of
the external field. [25, 168]

#### Chiroptical spectroscopies

Optical activity and circular dichroism are effects arising when chiral matter
interacts with polarized light. Chiral molecules possess a different complex
index of refraction for right- and left-handed circularly polarized (CP)
light. The real part determines the different refraction of CP light and also
the rotation of the plane of polarization of linearly polarized (LP) light,
while the imaginary part determines the difference in absorption of CP light
and the induced ellipticity of LP light. [181] At the molecular level, the
property underpinning these processes is the electric dipole–magnetic dipole
tensor $\mathbf{\beta}$ (also known as Rosenfeld tensor), that also connects
to the first order the induced electric dipole moment
$\bm{\mu}^{\mathrm{ind}}$ to the time derivative of the external magnetic
field $\bm{B}$ as well as the induced magnetic dipole moment
$\bm{m}^{\mathrm{ind}}$ to the time derivative of the external electric field
$\mathbf{E}$

$\displaystyle\mu^{\mathrm{ind}}_{i}(\omega)$ $\displaystyle=$
$\displaystyle\beta_{ij}(\omega)\dot{B}_{j}(\omega),$ (125) $\displaystyle
m^{\mathrm{ind}}_{i}(\omega)$ $\displaystyle=-$
$\displaystyle\beta_{ji}(\omega)\dot{E}_{j}(\omega).$ (126)

Note that we have restricted ourselves to isotropic samples where a
quadrupolar contribution that is non-zero for a single molecule vanishes after
averaging over molecular orientations. RT-TDSCF calculations of chiroptical
properties are based on Eq. (126) rather than on the direct simulation of an
interaction of molecules with circularly polarized light. The calculation
proceeds analogously to the computational protocol outlined here for electron
absorption spectroscopy: a molecule in its ground state $\mathbf{D}_{0}$ is
perturbed by an external electric field in the form of a $\delta$-pulse and
the induced _magnetic_ dipole moment

$\bm{m}^{\mathrm{ind}}(t)=\operatorname{Tr}[\mathbf{M}\mathbf{D}(t)]-\bm{m}^{\mathrm{static}},$
(127)

is evaluated in the course of the simulation. In Eq. (127),
$\bm{m}^{\mathrm{static}}=\operatorname{Tr}[\mathbf{M}\mathbf{D}_{0}]$ is the
static magnetic dipole moment and

$\mathbf{M}^{\mathrm{4c}}_{\mu\nu}=-\frac{1}{4c}\begin{pmatrix}\bm{0}&\braket{X_{\mu}}{(\bm{r}_{\mathrm{g}}\times\bm{\sigma})(\bm{\sigma}\cdot\bm{p})}{X_{\nu}}\\\
\braket{X_{\mu}}{(\bm{\sigma}\cdot\bm{p})(\bm{r}_{\mathrm{g}}\times\bm{\sigma})}{X_{\nu}}&\bm{0}\end{pmatrix},$
(128)

is the matrix representation of the magnetic dipole moment operator in the RKB
basis with $\bm{r}_{\mathrm{g}}=\bm{r}-\bm{R}_{\mathrm{g}}$ standing for the
electronic position operator relative to a fixed gauge $\bm{R}_{\mathrm{g}}$.
The induced magnetic dipole moment is transformed to the frequency domain and
used to calculate the Rosenfeld tensor via

$\beta_{ji}(\omega)=-i\frac{m^{\mathrm{ind}}_{i}(\omega)}{\mathcal{E}},$ (129)

where $\mathcal{E}$ is again the amplitude of the perturbing external field.
Chiroptical properties are notoriously sensitive to different parameters of a
calculations such as the choice of functional, basis set, conformation of the
molecule, solvent effects etc. [182]. Using relativistic RT-TDDFT it was shown
[84] on a series of model molecules – analogs of dimethyloxirane with the
oxygen atom replaced with heavier homologues (S, Se, Te, Po, Lv), that
relativity alone can change the sign of the spectral function, i.e. the factor
discriminating between the enantiomers. An example of such a spectrum is shown
in Figure 2 for dimethylpolonirane (PBE functional [178, 183, 184],
uncontracted Dyall’s aug-cVDZ basis [185, 186] for Po, and uncontracted
Dunning’s aug-cc-pVDZ [187, 188] for light elements). Therefore, relativistic
real-time methods should be an important tool in practical calculations of
chiroptical spectra, especially of molecules containing heavy elements.

Figure 2: ECD spectrum of polonium analogue of dimethyloxirane calculated
using relativistic 4c as well as 2c RT-TDDFT compared to a non-relativistc 1c
result. The inclusion of relativistic effects changes the sign of the spectrum
in the highlighted area demonstrating the need of proper treatment of
relativity when addressing chiroptical spectroscopic properties. Adapted from
Ref. 84 with permission. Copyright 2018, AIP.

###  8.2 Nonlinear optical properties

For a weak external field, the spectra resulting from the real-time
propagation will be equivalent to the results obtained using response theory.
However, in stronger fields, real-time simulations contain corrections of
higher orders. This is seen from comparing the perturbation expansion for the
induced dipole moment, schematically

$\bm{\mu}^{\mathrm{ind}}=\bm{\alpha}^{\mathrm{PT}}\bm{\mathcal{E}}+\bm{\beta}^{\mathrm{PT}}\bm{\mathcal{E}}^{2}+\bm{\gamma}^{\mathrm{PT}}\bm{\mathcal{E}}^{3}+\ldots,$
(130)

with the way the induced dipole moment from a real-time simulation is
processed, again schematically

$\bm{\mu}^{\mathrm{ind}}=\bm{\alpha}^{\mathrm{RT}}\bm{\mathcal{E}}=\left[\bm{\alpha}^{\mathrm{PT}}+\bm{\beta}^{\mathrm{PT}}\bm{\mathcal{E}}+\bm{\gamma}^{\mathrm{PT}}\bm{\mathcal{E}}^{2}+\ldots\right]\bm{\mathcal{E}}.$
(131)

In Eqs. (130) and (131) the indices PT and RT refer to perturbation theory and
real-time, respectively, and the molecular properties correspond to
polarizability ($\bm{\alpha}$), first hyperpolarizability ($\bm{\beta}$), and
second hyperpolarizability ($\bm{\gamma}$). While this feature of real-time
methods enables the study of strong-field effects in spectra, the properties
of higher orders are incorporated in the non-perturbative
$\bm{\mu}^{\mathrm{ind}}$ or $\bm{\alpha}^{\mathrm{RT}}$ and are not readily
available for further analysis. However, in some applications, it is desirable
to know the values of higher-order responses individually. This is also
possible to achieve using real-time methods by combining simulations with
various field strengths.

To show how a method for obtaining nonlinear responses from real-time
simulations can work, let us examine more closely the Taylor expansion of the
time-dependent induced dipole moment

$\mu_{i}(t)=\mu_{ij}^{(1)}(t)\mathcal{E}_{j}+\mu_{ijk}^{(2)}(t)\mathcal{E}_{j}\mathcal{E}_{k}+\mu_{ijkl}^{(3)}(t)\mathcal{E}_{j}\mathcal{E}_{k}\mathcal{E}_{l}+\ldots$
(132)

where $\mathcal{E}_{j}$ combines the amplitude and direction of the external
field, i.e. $\bm{\mathcal{E}}=\mathcal{E}\bm{n}$, and we defined the $n$-th
order contributions $\mu^{(n)}$ to the induced dipole moment. These
contributions are convolutions of the time-dependent (hyper)polarizability
tensors with the time dependence of the external field(s)

$\displaystyle\mu_{ij}^{(1)}(t)$
$\displaystyle=\int\\!dt_{1}\,\alpha_{ij}(t-t_{1})F(t_{1}),$ (133a)
$\displaystyle\mu_{ijk}^{(2)}(t)$
$\displaystyle=\frac{1}{2!}\int\\!dt_{1}\int\\!dt_{2}\,\beta_{ijk}(t-t_{1},t-t_{2})F(t_{1})F(t_{2}),$
(133b) $\displaystyle\mu_{ijkl}^{(3)}(t)$
$\displaystyle=\frac{1}{3!}\int\\!dt_{1}\int\\!dt_{2}\int\\!dt_{3}\,\gamma_{ijkl}(t-t_{1},t-t_{2},t-t_{3})F(t_{1})F(t_{2})F(t_{3}).$
(133c)

Again, the experimentally relevant quantities are the frequency-dependent
(hyper)polarizability tensors

$\displaystyle\alpha_{ij}(\omega)=$
$\displaystyle\int\\!dt_{1}\,\alpha_{ij}(t_{1})\,e^{-i\omega t_{1}},$ (134a)
$\displaystyle\beta_{ijk}(\omega_{1},\omega_{2})=$
$\displaystyle\int\\!dt_{1}\int\\!dt_{2}\,\beta_{ijk}(t_{1},t_{2})\,e^{-i\omega_{1}t_{1}}e^{-i\omega_{2}t_{2}},$
(134b) $\displaystyle\gamma_{ijkl}(\omega_{1},\omega_{2},\omega_{3})=$
$\displaystyle\int\\!dt_{1}\int\\!dt_{2}\int\\!dt_{3}\,\gamma_{ijkl}(t_{1},t_{2},t_{3})\,e^{-i\omega_{1}t_{1}}e^{-i\omega_{2}t_{2}}e^{-i\omega_{3}t_{3}}.$
(134c)

If we choose a harmonic external field,
$V^{\mathrm{ext}}(t)=\mathcal{E}\cos(\omega t)\bm{n}\cdot\mathbf{P}$, the
integrals in Eqs. (133) and (134) can be simplified to obtain expressions
relating $\mu^{(n)}$ to specific nonlinear optical (NLO) properties

$\displaystyle\mu_{ij}^{(1)}(t)$
$\displaystyle=\alpha_{ij}(-\omega;\omega)\cos(\omega t),$ (135a)
$\displaystyle\mu_{ijk}^{(2)}(t)$
$\displaystyle=\frac{1}{4}\left[\beta_{ijk}(-2\omega;\omega,\omega)\cos(2\omega
t)+\beta_{ijk}(0;\omega,-\omega)\right],$ (135b)
$\displaystyle\mu_{ijkl}^{(3)}(t)$
$\displaystyle=\frac{1}{24}\left[\gamma_{ijkl}(-3\omega;\omega,\omega,\omega)\cos(3\omega
t)+3\bar{\gamma}_{ijkl}(-\omega;\omega,\omega,-\omega)\cos(\omega t)\right].$
(135c)

The frequency-dependent molecular property tensors in equations (135) are the
dipole polarizability $\alpha_{ij}(-\omega;\omega)$, and higher-order
properties governing processes involving several photons, namely, the second
harmonic generation (SHG) coefficient $\beta_{ijk}(-2\omega;\omega,\omega)$,
the optical rectification (OR) coefficient $\beta_{ijk}(0;\omega,-\omega)$,
the third harmonic generation (THG) coefficient
$\gamma_{ijkl}(-3\omega;\omega,\omega,\omega)$ and the averaged degenerate
four-wave mixing (DFWM) coefficient
$\bar{\gamma}_{ijkl}(-\omega;\omega,\omega,-\omega)$. [189, 73, 74]

The workflow of the procedure for evaluating NLO properties from real-time
simulations [74, 26] is as follows

  1. 1.

Starting from a converged ground-state SCF, perform several real-time
simulations employing a cosine-shaped external field with different amplitudes
of the field, for example $\mathcal{E}_{1}=\mathcal{E}$,
$\mathcal{E}_{2}=2\mathcal{E}$, $\mathcal{E}_{3}=-\mathcal{E}$ and
$\mathcal{E}_{4}=-2\mathcal{E}$. Note that to improve the stability of time
evolution and smoothness of extracted responses, the cosine function is
multiplied with a linear envelope $\omega t/(2\pi)$ in the first period. [74]
Different envelopes with improved performance has also been suggested. [190]

  2. 2.

Calculate $\mu^{(n)}$ as derivatives of induced dipole moment

$\mu_{ij}^{(1)}(t)=\left.\frac{\partial\mu_{i}(t)}{\partial\mathcal{E}_{j}}\right|_{\bm{\mathcal{E}}=0},\quad\mu_{ijk}^{(2)}(t)=\left.\frac{1}{2}\frac{\partial^{2}\mu_{i}(t)}{\partial\mathcal{E}_{j}\partial\mathcal{E}_{k}}\right|_{\bm{\mathcal{E}}=0},\quad\mu_{ijkl}^{(3)}(t)=\left.\frac{1}{6}\frac{\partial^{2}\mu_{i}(t)}{\partial\mathcal{E}_{j}\partial\mathcal{E}_{k}\partial\mathcal{E}_{l}}\right|_{\bm{\mathcal{E}}=0},$
(136)

by means of numerical differentiation – a finite field method in each time
step. For example, the first- and second-order responses can be calculated
from simulations emplying fields from step 1) with precision of the order
$\mathcal{E}^{4}$ via

$\displaystyle\mu_{ij}^{(1)}(t)$
$\displaystyle=\frac{8\left[\mu_{i}(t,\mathcal{E}_{j})-\mu_{i}(t,-\mathcal{E}_{j})\right]-\left[\mu_{i}(t,2\mathcal{E}_{j})-\mu_{i}(t,-2\mathcal{E}_{j})\right]}{12\mathcal{E}_{j}},$
(137a) $\displaystyle\mu_{ijj}^{(2)}(t)$
$\displaystyle=\frac{16\left[\mu_{i}(t,\mathcal{E}_{j})+\mu_{i}(t,-\mathcal{E}_{j})\right]-\left[\mu_{i}(t,2\mathcal{E}_{j})+\mu_{i}(t,-2\mathcal{E}_{j})\right]}{24\mathcal{E}_{j}^{2}}.$
(137b)

  3. 3.

Fit the obtained $n$-th order induced dipole moment contributions to
analytical expressions in Eqs. (135) to evaluate numerical values of the NLO
properties.

An example of such a real-time finite field procedure is depicted in Figure 3
for the second-order response $\mu_{xxx}^{(2)}(t)$ of W(CO)5py, py = pyridine
at the 1eX2C level of theory (B3LYP functional [178, 179, 191, 192, 193],
Dyall’s uncontracted valence DZ basis set [180] for W, uncontracted aug-cc-
pVDZ basis [188] for the light elements). The fitting was used to determine
the second harmonic generation and optical rectification coefficients
$\beta^{\mathrm{SHG}}_{xxx}$ and $\beta^{\mathrm{OR}}_{xxx}$, respectively.
The figure is based on data underpinning Ref. [26] where it was shown that the
inclusion of relativistic effects contributed to about 35% of the final value,
highlighting the importance of a relativistic treatment in NLO applications
where heavy metal-containing compounds are of interest due to the favourable
electronic properties of the metallic centre [194, 195].

Figure 3: Non-linear optical properties from RT-TDDFT: fitting of second-order
induced dipole moment contribution $\mu_{xxx}^{(2)}(t)$ of W(CO)5py obtained
as a numerical second derivative in each time step to analytical expression
$\mu_{xxx}^{(2)}(t)=\tfrac{1}{4}\left[\beta^{\mathrm{SHG}}_{xxx}\cos(2\omega
t)+\beta^{\mathrm{OR}}_{xxx}\right]$ from which elements of the
hyperpolarizability tensor can be extracted. This figure is generated by the
authors based on data underpinning Ref. [26].

A special category of non-linear phenomena is high harmonic generation (HHG)
during which photons from a strong laser recombine into fewer photons of
higher energy via an interaction with a material. A HHG spectrum thus contains
peaks corresponding to multiples of the frequency of the laser field. HHG has
practical importance both as a spectroscopy technique as well as a means of
generating coherent high-frequency radiation. HHG presents a challenge for
theoretical modelling due to the necessity of using strong fields – several
orders of magnitude stronger than the applications discussed so far, thus
requiring a stable propagation, and due to the requirements on basis sets that
need to be able to describe electrons oscillating far from nuclei. The first
relativistic RT-TDDFT calculations of HHG were presented by De Santis et al.
using the PyBerthaRT program for Au2, capturing harmonics up to the 13th
order. [83]

###  8.3 Non-equilibrium spectroscopies

So far we have discussed molecular properties where formulations in terms of
perturbation theory exist, and are usually the preferred mode of calculation.
However, real-time methods are particularly well suited for the simulation of
experiments where the use of response theory would be too cumbersome. Such is
the case of non-equilibrium spectroscopies where more than one laser pulses
are used to drive the molecule. In the so-called pump–probe or transient
absorption (TA) spectroscopies, the first pulse (the pump) is used to excite
the molecule into a non-equlibrium state while the second pulse (the probe)
then measures the response of the driven molecule. By varying the time delay
between the pump and the probe it is possible to follow quantum dynamics of
electrons in molecules in real time. While a response theory-based description
for pump–probe experiments exists in the form of non-equlibrium response
theory, the ability of real-time methods to tailor the pulse shape to match
the experiments and handle strong fields offers a distinct advantage over
perturbative techniques.

In the case of transient absorption spectroscopies, two external pulses are
used in the simulation, the pump $\bm{\mathcal{E}}(t)$ and the probe
$\bm{\mathcal{F}}(t)$ as introduced in the Fock matrix in Eq. (61). The pump
first excites the molecule to a non-stationary excited state. This perturbed
state then evolves in time and its evolution is probed by the second pulse
applied after a time delay $\tau$. As an example, let us consider a set-up
with the pump pulse taking the form

$\bm{\mathcal{E}}(t)=\bm{n}\mathcal{E}(t)=\bm{n}\mathcal{E}\cos^{2}\left(\pi\frac{t-t_{0}}{T}\right)\sin(\omega_{0}t+\phi),$
(138)

with amplitude $\mathcal{E}$, polarization direction $\bm{n}$, and the pulse
shape defined by the carrier frequency $\omega_{0}$ of a sine wave, $\cos^{2}$
envelope, carrier–envelope phase (CEP) $\phi$, and time duration $T$. The
carrier frequency is usually tuned to an excitation energy of the molecule
which then becomes the primary excited state in the superposition state
created by the pump. However, even with relatively large amplitudes, the
ground state remains the most populated one.

For the probe, we use a broadband $\delta$-function pulse

$\bm{\mathcal{F}}(t)=\bm{m}\mathcal{F}(t)=\bm{m}\mathcal{F}_{0}\delta(t-(T+\tau)),$
(139)

applied at time $\tau$ after the pump pulse, that similarly to the case of
linear spectroscopies induces a time-dependent dipole moment that can be
processed to yield an absorption spectrum. However, since the initial state
now corresponds to the superposition state instead of pure ground state, the
spectrum contains imprints of quantum dynamics of the non-equilibrium state.
The probe pulse can be applied while the pump is still active, overlapping
regime, or after the pump has been turned off, non-overlapping regime.

The final TA spectra are obtained within the RT-TDDFT framework from the
differential induced dipole moment

$\Delta\mu^{\mathrm{TAS}}_{uv}(t)=\operatorname{Tr}\left\\{\mathbf{P}_{u}\left[\mathbf{D}_{v}^{\mathrm{pp}}(t)-\mathbf{D}_{v}^{\mathrm{p}}(t)\right]\right\\}\equiv\mu^{\mathrm{ind,pp}}_{uv}(t)-\mu^{\mathrm{ind,p}}_{uv}(t),\qquad
u,v\in\\{x,y,z\\},$ (140)

where $\mu^{\mathrm{ind}}_{uv}(t)$ denotes the expectation value of the dipole
operator. The computation of TA spectra involves performing two simulations
for recording the dipole moment at each time step; these simulations and their
quantities are denoted by p and pp subscripts, indicating that the real-time
propagation used pump-only pulse and pump together with the probe pulse,
respectively.

As an example, let us consider the TA spectrum of thiophene (PBE0 [178, 183,
184, 196] functional modified to contain 40% of Hartree–Fock exchange,
uncontracted aug-cc-pVXZ, X = T (S), D (C,H), basis [187, 188, 197]) depicted
in Figure 4. Here, the pump carrier frequency was set to correspond to the
first excitation energy while the X-ray absorption at sulfur L_2,3-edges was
investigated after the application of the probe. The spin–orbit splitting of
the sulfur 2p orbitals is visible in the spectrum which is thus correctly
described only by relativistic methods, in this case using the 4c
Dirac–Coulomb Hamiltonian (4c) and the amfX2C Hamiltonian (2c). The pump–probe
delay $\tau$ adds an extra degree of freedom to the TA spectrum which is then
normally plotted as a heat map where alternating low- and high-intensity
signals can be observed, tracing the dynamics of the superposition state as
induced by the pump pulse. Due to the increased computational cost of
obtaining such a spectrum – several spectra need to be calculated from
simulations with different $\tau$ in order to generate such a heat map –
efficient 2c relativistic methods are mandatory for these applications.

Figure 4: Transient absorption spectra of thiophene – comparison of 1c non-
relativistic (_orange_), 2c amfX2C (_blue_) and 4c Dirac–Coulomb (_cyan_)
Hamiltonian-based RT-TDDFT. (a) TAS at $\tau=0$; (b) Variation in 1c TAS
spectra with pump–probe delay $\tau$; (c) Variation in 2c TAS spectra with
pump–probe delay $\tau$. All spectra are obtained with a damping factor
$\Gamma=0.01$ au. Figure reprinted from Ref. [27] under CC-BY-4.0 licence.

##  9 Conclusion and perspectives

Real-time methods are based on a direct integration of the quantum mechanical
equation of motion. In non-relativistic quantum chemistry, they gained
popularity in previous decades due to their ability to describe phenomena
ranging from linear response properties to interaction with strong laser
fields and time-resolved spectroscopies – areas at the forefront of
experimental research. In relativistic quantum chemistry, the pioneering
theory development, computational implementation, and first applications
arrived later and the field has yet to catch up with the breadth of the scope
of applications of its non-relativistic counterparts.

In this chapter we summarized some of the advances of relativistic real-time
methods in quantum chemistry while restricting ourselves to mean-field methods
and pure electron dynamics. It has been our ambition that our introduction
explains the fundamental principles of this methodology and inspires readers
to join in this rapidly developing and promising field.

## References

  * [1] P. B. Corkum, F. Krausz, Attosecond science, Nature Phys. 3 (6) (2007) 381–387.
  * [2] M. Nisoli, P. Decleva, F. Calegari, A. Palacios, F. Martín, Attosecond electron dynamics in molecules, Chem. Rev. 117 (16) (2017) 10760–10825.
  * [3] X. Li, N. Govind, C. Isborn, A. E. DePrince III, K. Lopata, Real-time time-dependent electronic structure theory, Chem. Rev. 120 (18) (2020) 9951–9993.
  * [4] B. S. Ofstad, E. Aurbakken, Ø. S. Schøyen, H. E. Kristiansen, S. Kvaal, T. B. Pedersen, Time-dependent coupled-cluster theory, WIREs Computational Molecular Science (2023).
  * [5] T. Helgaker, S. Coriani, P. Jørgensen, K. Kristensen, J. Olsen, K. Ruud, Recent advances in wave function-based methods of molecular-property calculations, Chem. Rev. 112 (1) (2012) 543–631.
  * [6] H.-D. Meyer, U. Manthe, L. Cederbaum, The multi-configurational time-dependent hartree approach, Chem. Phys. Lett. 165 (1) (1990) 73–78.