MSc Chemistry
Molecular Sciences
Master Thesis
Beyond the Electric Dipole Approximation
in Relativistic Simulations:
Applications in Electronic Circular
Dichroism and Convergence Difficulties of
the Multipole Expansion
by
Martin van Horn
10992413
42 EC
May 2019 – December 2019
Supervisor/Examiner: Daily supervisor:
Examiner:
Prof. dr. L. Visscher Prof. dr. T. Saue
Prof. dr. P. Gori-Giorgi
Amsterdam Institute for Molecules,
Medicines and Systems
Preface
This thesis can be regarded as a continuation of unpublished work by Nanna List, Timoth´e Melin and Trond Saue. In this unpublished work, the quantum chemical code DIRAC was augmented with a mod-ule for the relativistic simulation of X-ray spectroscopy that describes the interaction of light and matter beyond the electric dipole approximation. The first two sections of this thesis are mainly concerned with X-ray spectroscopy, while the third section is focused on electronic circular dichroism. In the first two sections, the code of List and co-workers was used for all calculations and in the third section this code is extended to describe electronic circular dichroism. The derivation and implementation of the equations in the third section were originally carried out in this thesis, while the equations in the former two sections are inspired from the aforementioned work of List and co-workers.
The trick to forgetting the big picture is to look at everything close-up.
Abstract
The main topic of this thesis are the applications of semi-classical light-matter interaction operators be-yond the electric dipole approximation, which can either be expanded in terms of multipole moments (truncated interaction operator), or utilized in its sinusoidal form (full interaction operator). Within con-temporary scientific literature, several authors have applied operators of this type within non-relativistic theory in the calculation of X-ray absorption intensities, or electronic circular dichroism. Only the quan-tum chemical code DIRAC is capable of describing the semi-classical light-matter interaction operator beyond the electric dipole approximation within a relativistic framework, although this code lacks the capacity to describe electronic circular dichroism. Furthermore, the truncated interaction operator intro-duced divergences of the oscillator strengths, while the full interaction operator required additional basis set testing. Therefore, the aim of this thesis is threefold: firstly, the divergences are studied that are in-troduced by the truncated interaction operator; secondly, the basis set requirements of the full interaction operator are tested; and thirdly, the aforementioned implementation is extended to describe electronic circular dichroism. From the various tests with the truncated interaction operator could be concluded that the transition moments in fact did not diverge, but rather converge extremely slowly, suggesting that this kind of operator is impractical to use for the calculation of X-ray absorption intensities. Moreover, it was found that the transition moments of the full interaction operator were stable with respect to basis set augmentation, suggesting that this operator has equivalent basis set requirements as the energy. From a test calculation on hydrogen peroxide could be concluded that the ECD implementation was functional, as the calculations reproduced the reference values from the dipole limit, although further testing needed to be performed to be certain about this conclusion.
Contents
Abstract i
1 Introduction 1
2 Theory 4
2.1 Relativistic Interaction Operators . . . 4
2.1.1 Dirac Hamiltonian . . . 4
2.1.2 Effective Interaction Operator . . . 6
2.2 Semi-Classical Light-Matter Interaction Operator . . . 7
2.2.1 Classical Electromagnetic Waves . . . 7
2.2.2 Full Interaction Operator . . . 10
2.2.3 Truncated Interaction Operator: Generalized Length Representation . . . 12
2.2.4 Truncated Interaction Operator: Generalized Velocity Representation . . . 15
2.3 Absorption Cross Section . . . 16
2.3.1 Full Interaction Operator: Oriented Absorption Cross Section . . . 16
2.3.2 Full Interaction Operator: Isotropic Absorption Cross Section . . . 17
2.3.3 Truncated Interaction Operator: Generalized Length Representation . . . 19
2.3.4 Truncated Interaction Operator: Generalized Velocity Representation . . . 20
2.4 Orbitals and Basis Sets Suitable for X-ray Transitions . . . 21
2.4.1 Completeness-Optimized Basis Sets . . . 21
2.4.2 Improved Virtual Orbital . . . 22
3 Computational Details 24 3.1 Truncated Interaction Operator . . . 24
3.1.1 Generalized Length Representation . . . 24
3.1.2 Generalized Velocity Representation . . . 26
3.2 Full Interaction Operator . . . 27
3.2.1 Basis Set Test . . . 27
3.2.2 Electronic Circular Dichroism . . . 28
4 Results 29 4.1 Truncated Interaction Operator . . . 29
4.1.1 Generalized Length Representation . . . 29
4.1.2 Generalized Velocity Representation . . . 32
4.2 Full Interaction Operator . . . 37
4.2.1 Basis Set Test . . . 37
4.2.2 Electronic Circular Dichroism . . . 41
6 Outlook 44 7 Acknowledgements 45 8 Bibliography 46 Appendices 51 A Response Theory 51 B Property Integrals 52
1
Introduction
Spectroscopy is a ubiquitous experimental technique that has a wide variety of applications in chemistry, ranging from the analysis of molecular properties to the characterization of unknown compounds.1Each spectroscopic method is based on the interaction of light with matter, although these methods may differ with respect to the range of wavelengths used to excite molecular systems. At short wavelengths, sam-ples are analysed with X-ray absorption spectroscopy (XAS), which is a collective term for a multitude of spectroscopic techniques characterized by the nature of the excitations it aims to study. For example, X-ray absorption near-edge structure (XANES) describes excitations close to the ionization threshold, whereas extended X-ray absorption fine structure (EXAFS) gathers information from electrons that are scattered by high energy photons.2,3 In the lower X-ray regime, core electrons are excited to low-lying bound states, giving rise to the prepeaks. Compared to other spectroscopic methods, XAS provides more local information on molecular properties due to the limited spatial extension of core orbitals. The results of these measurements can be used in conjunction with quantum chemical calculations to iden-tify the spin state, the local coordination environment, or the oxidation states of metallic complexes.4–8 Therefore, improving the theory that describes the absorption of electromagnetic radiation is essential as to gain information about these properties.9
Classically, electromagnetic waves in vacuum can be modeled as linear combinations of sinusoidal functions that obey the homogeneous Maxwell equations.10Within a quantum mechanical framework, observable properties are represented by Hermitian operators and their eigenvalues correspond to mea-surement outcomes, which are probabilistic in nature. Classical mechanics can be connected to the quan-tum realm by promoting the observable properties to Hermitian operators, yielding the semi-classical matter-light full interaction operator in the case of electromagnetic waves.
This theoretical description of absorption phenomena can be simplified by employing the electric dipole (ED) approximation, which is derived by expanding the semi-classical matter-light interaction operator in terms of the wave vector. Subsequently, this expansion is truncated at zeroth order, resulting in the ED approximation. At this level of approximation, the interaction operator describes the pertur-bation induced by a uniform electric field with a negligible magnetic component. Ignoring the magnetic field at zeroth order is justified by the homogeneous Maxwell equations, which prescribe that the mag-netic field is at least linear in the wave vector. It can be shown that the transition dipole moment can either be calculated as matrix elements of the momentum or position operator, referred to as the length and velocity representation, respectively.11–13 Provided that the extent of the system under considera-tion is small compared to the wavelength of light, the electric dipole approximaconsidera-tion can yield reasonable accuracy for the calculation of absorption intensities.
However, in the X-ray regime, more thorough theory ought to be used, due to the short wavelength of radiation.14,15At the length scale of the system under consideration, the spatial variation of the electro-magnetic wave cannot be neglected, giving rise to dipole-forbidden transitions. Transitions of this type, also referred to as multipole transitions, have different symmetries and dependencies on the incident angle of the electromagnetic radiation. This property was exploited by Drager et al. in their analysis of the multipole character of X-ray transitions in iron compounds.14 By irradiating crystals of these
compounds at various incident angles and comparing the results with the mathematical expressions of multipole moments, they could prove the presence of these transitions experimentally.
The mathematical expressions of multipole moments emerge from the semi-classical full interaction operator by truncating its Taylor expansion at higher orders than zero. As a result, the transition moments can be decomposed in contributions from individual multipole moments, thus providing a language to understand these quantities. The truncated interaction operator can be expressed in two equivalent rep-resentations: the generalized length or velocity representation, in analogy with the electric dipole ap-proximation. The generalized velocity representation follows directly from the Taylor-expansion of the semi-classical light-matter full interaction operator and it can be converted to the length representation by introducing the multipolar gauge.16–22By employing this gauge, the electromagnetic potentials are expanded in terms of their respective fields.
The theoretical description of multipole moments has been met with some complications, as trun-cation above zeroth order induces a gauge-origin dependency. Upon varying this origin, each multipole moment introduces additional terms that depend on lower-ordered multipoles. Berry and co-workers proposed an algorithm that determines a gauge-origin by minimizing these additional terms.6In most cases this origin is then positioned at places with high transition density, e.g. at the nucleus where the transition occurs. However, it has been shown by Atkins et al. that for dipole-forbidden transitions, this algorithm possibly locates the gauge-origin far from the nuclei.23 Furthermore, this scheme has difficulties with metallic clusters and delocalized transitions. In a seminal paper by Bernadotte and co-workers, the problem of gauge-origin dependency was alleviated by expanding the oscillator strength in orders of the wave vector, whereas formerly, the transition moments were expanded in terms of the wave vector.24This alteration can be justified by noting that the oscillator strength is an experimentally observable quantity, whereas the transition moment is merely an auxiliary quantity. It has been shown by Lestrange et al. that this method can possibly yield negative oscillator strengths, which was ascribed to an unbalanced choice of basis set by Sørensen et al.25,26
Alternatively, the semi-classical light-matter interaction operator can be utilized without truncation. Within the contemporary scientific literature, several attempts have been made to implement the full semi-classical light-matter interaction operator in this form. List et al. and Khamesian et al. imple-mented the full interaction operator in quantum chemical software, while utilizing a Gaussian basis
set.27–29 By doing so, integrals over this operator could be expressed as a Fourier transformation of
Gaussian functions, thereby yielding simple and analytic expressions. Sørensen et al. proceeded slightly differently, as they computed these integrals numerically, utilizing Gauss-Hermite quadrature.30 These attempts to go beyond the electric dipole approximation (BED) were carried out with non-relativistic quantum chemical software. Typical X-ray transitions involve core orbitals, which are susceptible to relativistic effects, due their short distance from the nuclei. Therefore, a relativistic theory is required to describe transitions of this kind properly.
In recent unpublished work, List and co-workers implemented both the full and truncated interac-tion operators within a relativistic framework, resulting in a mathematically simple formulainterac-tion of this operator.31 The implementation of the full interaction operator yielded reasonable results for the
cal-culations of oscillator strength in the X-ray regime, as these calcal-culations gave small corrections to the electric dipole approximation. Additional basis set studies needed to be performed to investigate the basis set requirements of the full interaction operator. However, peculiar results were obtained for the same calculations carried out with the truncated interaction operator in both the generalized length and velocity representation: upon increase of the truncation order, the absolute value of the contributions to the oscillator strengths would monotonically increase, indicating divergence in the limit of infinite order. Hypothetically, the divergence of the oscillator strengths in the X-ray regime can be explained by nu-merical noise, bugs in the code, basis set error, a finite radius of convergence of the transition moments, or the incompatibility of this theory with a basis set expansion.
Besides its uses in the modeling of XAS, the full interaction operator can also be applied to electronic circular dichroism (ECD), which can be used in combination with quantum chemical calculations to determine the absolute configuration of chiral molecules.29XAS can be modeled by implementing an operator representation of a linearly polarized electromagnetic wave, while ECD can be described by the difference between the oscillator strengths of right and left handed circularly polarized light. It should be noted that the implementation of both properties require the same integrals, thus the implementation of the full interaction operator can easily be extended to describe ECD.
The aim of this thesis is threefold: firstly, the divergences are studied that are introduced by the truncated interaction operator; secondly, the basis set requirements of the full interaction operator are tested; and thirdly, the implementation of List et al. is extended to describe electronic circular dichroism. This thesis can be thought of as a continuation of the work by List and co-workers, hence aforementioned interaction operators will be described in a relativistic framework. It should be stressed that the first two objectives were carried out with existing code, whereas the implementation and the derivation of the equations concerning of ECD were originally carried out in this thesis.
2
Theory
All calculations in this thesis were performed with DIRAC, which is a quantum chemical code, designed for calculations with relativistic four component or two component methods.32The aim of this section is to describe all equations that are implemented in this code and show the principles from which they are derived. In the first subsections, the interaction operators will be derived, after which the derivation of the absorption cross sections are carried out. Furthermore, it should be noted that in the entire theory section, SI units will be used.
2.1 Relativistic Interaction Operators
2.1.1 Dirac Hamiltonian
In relativistic theory, the energy of a classical point particle is described by expression (2.1),
E =pm2c4+ π2c2+ qφ; π = p − qA, (2.1)
where c represent the speed of light and m, q, p and E the mass, charge, momentum and energy of the point particle, respectively.33 The vector potential, A, and scalar potential φ enter this expression through minimal coupling, that is:
p → p − qA; E → E + qφ. (2.2)
Note how the classical kinetic energy can be retrieved from (2.1) in the limit mcπ 1, by Taylor-expanding the square root term and truncating it at first order:
mc2 r 1 + π 2 m2c2 ≈ mc 2+ π2 2m− π4 8m3c2 + O( π6 m5c4) (2.3)
In principle, a relativistic Hamiltonian can be derived from (2.1), by promoting the observable quantities to Hermitian operators. However, in operator form, the square root term in (2.1) is only defined as an infinite Taylor series, which needs to be truncated at finite order, thereby losing Lorentz-invariance.33 Moreover, at higher orders, this expansion becomes more and more problematic, due to the non-commutativity of the momentum and the vector potential. This behaviour can be illustrated by expanding the quadratic term:
π2 = (p − qA)2 = p2− qp · A − qA · p + q2A2. (2.4) At higher orders, the amount of additional terms would get even higher, thus, in practice, the operator form of equation (2.1) is rarely used.
In Dirac’s approach, a linearization of (2.1) is proposed in terms of α and β, circumventing the square root term altogether,
where (2.5) also allows negative energy solutions. After moving qφ to the left-hand side of equations (2.5) and (2.1) and squaring both expressions, one finds that they can only be equal if αiand β obey the
algebraic rules below:
α2i = β2= 1 αiαj = −αjαi; i 6= j αiβ = −βαi m2c4+ π2c2= c2αiαjπiπj+ c3mαiβπi+ c3mβαjπj+ m2c4β2.
In this equation, the Einstein summation convention is used, implying that all components are summed for terms with repeated indices. The Dirac matrices, shown in (2.6), are one possible representation of this algebra, αk = 02 σk σk 02 ! ; k = 1, 2, 3; β = I2 02 02 −I2 ! , (2.6)
where appears the Pauli matrices, {σk}, the 2x2 identity matrix, I2 and the 2x2 zero matrix, 02. The
Pauli matrices describe the intrinsic angular momentum of a spin12 particle and are defined as:
σ1 = 0 1 1 0 ! ; σ2 = 0 −i i 0 ! ; σ3 = 1 0 0 −1 ! . (2.7)
The intrinsic angular momentum is an axial vector, which is anti-symmetric with respect to time-reversal, hence the Pauli matrices have this symmetry as well. It follows naturally from their definition in equation (2.6) that the Dirac matrices must have the same behaviour as the Pauli matrices with respect to time-reversal.
After quantizing expression (2.1), one obtains the Dirac Hamiltonian:
ˆ
hR= βmec2+ c(α · ˆp) + ec(α · A) − eφ. (2.8)
This Hamiltonian is a 4x4 matrix operator and its action is only well-defined on spinors, which are four component vector functions, analogous to the wave function in non-relativistic theory. The first and third component of a spinor can be identified with spin 12, whereas the second and fourth component describe spin −12. Furthermore, for positive energy solutions, the first two components are a factor of mc larger than the lower two components, hence both the upper and lower halves are classified as two-component spinors: the large and small component, respectively.
Interactions with external fields enter the Dirac Hamiltonian through the potentials, which appear as two linear terms in (2.8). It should be stressed that the minimal coupling procedure for the non-relativistic Pauli-Hamiltonian result in three more interaction terms than the Dirac Hamiltonian,
ˆ hN R= pˆ 2 2me + e 2me (p · A + A · p) +e 2A2 2me − eφ + e~ 2me (σ · B) (2.9)
, where the final term represents the spin-Zeeman interaction and the term quadratic in A the diamagnetic interaction.34The origin of this additional complexity is the mixing of relativistic and non-relativistic
representations. Magnetism is a purely relativistic phenomenon, since it vanishes in the non-relativistic limit, although for pragmatic reasons it is reasonable to apply it in non-relativistic theory. Why mag-netism vanished in the non-relativistic limit, can be understood by considering the Maxwell’s equations :35 ∇ · B = 0 ∇ × E + ∂B ∂t = 0 ∇ · E = 1 0 ρ ∇ × B = µ0J + 1 c2 ∂E ∂t −−−→ c→∞ ∇ · B = 0 ∇ × E +∂B ∂t = 0 ∇ · E = 1 0 ρ ∇ × B = 0.
Note that the magnetic field has both vanishing divergence and curl in the non-relativistic limit. Com-bined with the constraint that physical fields should vanish at infinite distance from the source, one can conclude that the magnetic field has to be zero everywhere, because there is no vector field to which all three conditions apply.10Therefore, electrostatics is a mathematically consistent non-relativistic limit of electrodynamics.
2.1.2 Effective Interaction Operator
We start from the Lambert-Beer law, which gives the relation between the attenuation of light, II 0, the propagation path, l, the number density of absorbing molecules, N , and the absorption cross section σ.
I I0
= e−N σl (2.10)
The latter is a measure of the absorption probability and it can be related to the microscopic transition rate, obtained from quantum mechanical calculations. In equation (2.11), this relation is shown,
I(ω)σ(ω) = ~ωwf ←i(ω); I =
1 20cE
2, (2.11)
where wf ←i(ω) denotes the microscopic transition rate, E the field strenght, 0 the permittivity of
vac-uum and σ the absorption cross section. Starting from an interaction operator of the form,
ˆ
V (t) = ˆV (ω)e−iωt+ ˆV (−ω)e+iωt; V (−ω) = ˆˆ V†(ω), (2.12) and using time-dependent perturbation theory one obtains the expression for the transition rate,
wf ←i(ω) =
2π
~2|hf | ˆV (ωf i) |ii |
2f (ω, ω
f i, γf i), (2.13)
in which f (ω, ωf i, γf i) represents a lineshape function. Therefore, an interaction operator of the form in
(2.12) needs to be found, in order to extract ˆV (ω) and determine the transition rate. Hermiticity of ˆV (t) is ensured by imposing the relation between ˆV†(ω) and ˆV (−ω). Expression (2.13) can now be plugged in (2.11) to find an expression for the absorption cross section,
σ(ω) = 2πω ~I(ω) |hf | ˆV (ωf i) |ii |2f (ω, ωf i, γf i) = πω ~ |hf | ˆT (ωf i) |ii |2f (ω, ωf i, γf i), (2.14)
where in the last equality, ˆV (ω) is substituted with −12E ˆT (ω), in order to remove the dependency on the field strength. In the following sections, various interaction operators will be defined by applying the relations shown above.
2.2 Semi-Classical Light-Matter Interaction Operator
2.2.1 Classical Electromagnetic Waves
The most general form of monochromatic electromagnetic radiation is described by equation (2.15),
E = E11sin(k · r − ωt + δ1) + E22sin(k · r − ωt + δ2); E1= E cos θ; E2 = E sin θ
(2.15) B = E1 c 2sin(k · r − ωt + δ1) − E2 c 1sin(k · r − ωt + δ2); E = q E2 1 + E22,
where c, ω, δ and E represent the speed of light, angular frequency, phase and field strength, respectively. The polarization and wave vectors of this system are described by the basis of the orthonormal unit vectors 1, 2and ˆk, the latter defined as ˆk = ckω. These vectors span a right-handed coordinate system,
that is:
1× 2 = ˆk (2.16)
Consider two specific combinations of δ1, δ2and θ:
1. δ2 = δ1+ mπ; m ∈ Z: substituting this in (2.15) yields,
E = E Esin(k · r − ωt + δ1); E = cos θ1+ (−1)msin θ2 (2.17)
B = E
cBsin(k · r − ωt + δ1); B = cos θ2− (−1)
msin θ
1 = ˆk × E,
where E, B and ˆk form a right-handed coordinate system. The physical picture of
equation (2.17) corresponds to a linearly polarized bundle of radiation at a polarization angle θ relative to 1.
2. δ2 = δ1+12π + mπ; m ∈ Z and θ = 14π: we then get,
E = √1 2E[1sin(k · r − ωt + δ1) + (−1) m 2cos(k · r − ωt + δ1)] (2.18) B = 1 c√2E[2sin(k · r − ωt + δ1) − (−1) m 1cos(k · r − ωt + δ1)],
resulting in circularly polarized light. This type of light consists of two equally contribut-ing, perpendicular linearly polarized components that are out of phase byπ2.
Figure 1. Linearly polarized light at t = 0 , δ1= 0 and δ2= π: electric component in red and magnetic in blue.
(a) Right handed circularly polarized light (b) Left handed circularly polarized light
Figure 2. Electric component of right and left handed circularly polarized light, at t = 0 and δ1= 0.
Due to its helical shape, circularly polarized light is chiral, implying that this type of radiation can not be superimposed on its reflected image. Therefore, a distinction has to be made between right and left-handed circularly polarized light. Which handedness is assigned to which type of light is a matter of convention. In this thesis, the convention is used that is recommended by the IUPAC: the time in expression (2.18) is fixed and the rotation of the electric field vector is followed, while looking against the direction of the wave vector.36The time can be conveniently fixed at −ωt + δ1 = 0, such that only
the spatial part of the trigonometric functions remains. At this fixed time, the electric field vector traces a helix in its spatial dimension with its wave vector parallel to the screw axis. When looking at the source, circularly polarized light is right-handed if the vector rotates anti-clockwise, whereas it is left-handed if it rotates in the opposite direction. In equation (2.18), the handedness is determined by the parameter m and for even values, this parameter returns left-handed polarization, while the opposite is true for odd values.
The intensity of (2.15) can be obtained by averaging the magnitude of the Poynting vector over a full oscillation,
I = hSi ; S = 1
µ0
(E × B), (2.19)
where appears the Poynting vector, S and its magnitude S.10 By plugging (2.15) into the Poynting vector, we obtain four different terms, each containing a distinct unit vector, which can be computed by using the algebra of the right handed coordinate system.
S = 1 µ0
(E × B) = 1 µ0
(E cos θ1sin(k · r − ωt + δ1) + E sin θ2sin(k · r − ωt + δ2)) (2.20)
× (E c cos θ2sin(k · r − ωt + δ1) − E c sin θ1sin(k · r − ωt + δ2)) =E 2 cµ0 ˆ
k(cos2θ sin2(k · r − ωt + δ1) + sin2θ sin2(k · r − ωt + δ2))
By averaging this expression over an oscillation, T , one obtains:
I = E 2 cµ0 1 T Z T /2 −T /2
dt[cos2θ sin2(k · r − ωt + δ1) + sin2θ sin2(k · r − ωt + δ2)] =
c0
2 E
2 (2.21)
Note how the dependency on δ1, δ2 and θ dissapears after averaging, indicating that the energy of
monochromatic radiation is independent of its polarization or phase.
In a quantum mechanical framework, the influence of electromagnetic fields enter the Hamiltonian in the form of potentials, as was pointed out in section 2.1.1. The fields from (2.15), need to be expressed as potentials, by applying Helmholtz’s theorem. According to this theorem, all vector fields can be decomposed into a divergence free (solenoidal) and curl free (irrotational) component, if the curl and divergence of this vector field go faster to zero than r12 as r → ∞.10 A vector field, F , to which this constraint applies, can be expressed as,
F = −∇U + ∇ × W , (2.22) where U (r) = 1 4π Z d3r0 D(r 0) |r − r0|; W (r) = 1 4π Z d3r0 C(r 0) |r − r0| (2.23) D(r) = ∇ · F ; C(r) = ∇ × F .
A generic electromagnetic field can be decomposed by considering the structure of the Maxwell equations, shown in (2.10). These equations imply that the magnetic field is purely solenoidal, hence it can be expressed as the curl of some vector function,
where the last equality holds, because the irrotational part of a vector field does not contribute to the curl. This vector function can be identified with the vector potential. In expression (2.24), A⊥, represents
the solenoidal part of the vector potential. The Helmholtz decomposition of the electric field is less straightforward, because this field is neither purely solenoidal, nor irrotational. We start by substituting equation (2.24) in Faraday’s law, resulting in (2.25):
∇ × E + ∂(∇ × A)
∂t = 0. (2.25)
In this expression, the irrotational part of the vector potential can not be ignored, because in general, the electric field contains an irrotational component. Electromagnetic fields are smooth functions, implying that their partial derivatives commute, so equation (2.25) can be rewritten as,
∇ × (E +∂A
∂t ) = 0, (2.26)
in which the quantity inside the brackets is purely irrotational. Therefore, this vector field can be ex-pressed as the gradient of a scalar field, φ. The Helmholtz decomposition of a generic electromagnetic field gives,
E = −∂tA⊥− ∂tAk− ∇φ; B = ∇ × A⊥, (2.27)
where Akrepresents irrotational part of the vector potential.
Electromagnetic waves, such as (2.15), do not vanish at infinity, suggesting that Helmholtz’s theorem cannot be applied to systems of this type. However, provided that the vector field under consideration is differentiable, it is still possible to decompose such fields into divergence and curl-free components, although (2.23) can not be used to construct U and W , as these integrals generally diverge if the condi-tion at infinity is not met.10In practice, φ and A are calculated for electromagnetic waves by applying (2.27) in reverse.
2.2.2 Full Interaction Operator
In the case of monochromatic electromagnetic waves in vacuum, sources lie at infinity, implying that both E and B are purely solenoidal. Furthermore, the Coulomb gauge is chosen for this derivation, thus the vector potential is solenoidal as well. The Helmholtz decomposition in (2.27) can now be rewritten as,
E = −∂tA⊥; B = ∇ × A⊥, (2.28)
where ∇φ can not contribute to the electric field, as it is purely irrotational by construction. The possi-bility remains for the scalar potential to be constant, as a constant would have a vanishing gradient. Out of convenience, the value of the scalar potential is chosen to be zero. By applying (2.28) in reverse, the potentials corresponding to (2.15) can be isolated:
A = −E1
ω 1cos(k · r − ωt + δ1) + −E2
After plugging (2.29) in the time-independent Dirac Hamiltonian, equation (2.8), one obtains the inter-action operator, ˆV (t), ˆ V (t) = −eE 2ω ((cα · )e i(k·r−ωt)+ (cα · ∗ )e−i(k·r−ωt)) (2.30) ˆ V (ω) = −eE 2ω (cα · )e
ik·r; = cos θeiδ1
1+ sin θeiδ22,
where Eulers formula was used to convert the cosines in complex exponentials, thereby creating an interaction operator of the desired form. As was indicated in section 2.1.2, the intensity in the absorption cross section can be eliminated by defining an effective interaction operator of the form ˆV = −12E ˆT (ω) and plugging it in expression (2.14). By plugging in this definition, the absorption cross section can be expressed in terms of the effective operator.
ˆ
Tf ull(ω) =
e
ω(cα · )e
i(k·r); = cos θeiδ1
1+ sin θeiδ22 (2.31)
Due to the lack of time-reversal symmetry, this form of the interaction operator is not suited for im-plementation in DIRAC, which uses a quaternion scheme. Quaternions are generalizations of complex number and a generic quaternion has the following form and multiplication rules:
q = r + s˘ı + t˘ + u˘k; ˘ı2 = ˘2= ˘k2= ˘ı˘˘k = −1, (2.32) in which r,s,t and u are real numbers and ˘ı,˘ and ˘k represent the mutually anticommuting quaternion units.37 An isomorphism exists between the quaternions and SU (2), which is the group containing the Pauli matrices. By utilizing this relation, the interaction operator from (2.8) can be expressed in terms of quaternions: ˆ V (t) = −e φ 0 0 φ ! − ec˘ı 0 Az Az 0 ! − ec˘ 0 Ay Ay 0 ! − ec˘k 0 Ax Ax 0 ! . (2.33)
In this form, the matrix representation of the interaction operator can be block diagonalized.38 Time-reversal symmetry is required to carry out this transformation, thus equation (2.31) needs to be made time-reversal symmetric by splitting it into a Hermitian and anti-Hermitian part and applying the identity,
ˆ A = 1 2( ˆA + ˆA †) +1 2( ˆA − ˆA †), (2.34)
where ˆA denotes some generic operator. Decomposing the effective interaction operator for linearly polarized light yields:
ˆ Tf ull(ω) = e ω(cα · L)e i(k·r+δ1) ˆ Tf ull,H(ω) = ωe(cα · L) cos(k · r + δ1) ˆ
Tf ull,A(ω) = ωe(icα · L) sin(k · r + δ1)
(2.35)
L= cos θ1+ (−1)msin θ2,
ˆ Tf ull(ω) = e ω(cα · C)e i(k·r+δ1) ˆ
Tf ull,H(ω) = ω√e2(cα · 1cos(k · r + δ1) − (−1)mcα · 2sin(k · r + δ1))
ˆ
Tf ull,A(ω) = ω√e2(icα · 1sin(k · r + δ1) + (−1)micα · 2cos(k · r + δ1))
(2.36) C = 1 √ 2(1+ (−1) mi 2).
The Dirac matrices, α, are time-reversal anti-symmetric, hence the anti Hermitian part of ˆT is time-reversal symmetric and the Hermitian part is anti-symmetric with respect to this symmetry operation. Therefore, a factor of i is inserted in the Hermitian part in order to be implemented in the quaternion scheme.
The effective interaction operator is not Hermitian, although this property can be imposed in the limit, k · r 1, provided the proper value of δ1is chosen. In this limit, ˆT (ω) becomes:
ˆ Tf ull(ω) = e ω(cα · )e i(k·r+δ1) −−−−→ k·r 1 ˆ Tf ull(ω) = e ω(cα · )e iδ1. (2.37)
In the generalized velocity representation δ1has to be set equal to 0, to ensure Hermiticity of the dipole
term.
2.2.3 Truncated Interaction Operator: Generalized Length Representation
Most molecular properties can be understood as modifications of multipole moments by external elec-tromagnetic fields.36Electric and magnetic multipole moments indicate the intricacy of the charge and current distribution. Potentials that are generated by these distributions can be expanded in terms of multipoles. For example, a point charge at the origin only contains a monopole term that equals the total charge. More complicated charge configurations have contributions from higher-ordered electric and magnetic multipoles, defined as Q[n]j
1···jn and m [n] i;j1···jn−1in equation (2.38) , Q[n]j 1···jn = Z d3rrj1rj2· · · rjnρ(r, t); m [n] i;j1···jn−1 = n n + 1 Z d3rrj1rj2· · · rjn−1[r × j(r, t)]i (2.38) where j(r, t) represents the current density, ρ(r, t) the charge density and the indices {jm} run over
the Cartesian components of r.10 By truncating the multipole expansion at finite order, complicated electromagnetic potentials can be simplified, without losing the essential characteristics of the system.
Moreover, the truncated multipole expansion gains accuracy if the system is observed from large distances and at zeroth order this expansion yields the electric dipole approximation. Higher terms in this expansion are called quadrupoles, octupoles, hexadexapoles, etc., each of which induces a gauge-origin dependence, if the multipole expansion is truncated.24
It is desirable to circumvent this problem, because observable quantities always have gauge invari-ance, although potentials lack this property, implying that they are not uniquely defined. By utilizing this redundancy of description, calculations involving electromagnetic interactions can be made much easier
to handle. For example, one of the gauges that is employed in this report, introduces a multipole expan-sion of the potentials by expanding them in terms of their respective fields.16,17,36The derivation of this gauge usually starts with the electromagnetic potentials being Taylor expanded, centered at the point a. In the equation below, a prime indicates a dummy variable, which is set equal to a after differentiation.
φ(r, t) = ∞ X n=0 1 n![(δ · ∇ 0)nφ(r0, t)] r0=a A(r, t) = ∞ X n=0 1 n![(δ · ∇ 0)nA(r0, t)] r0=a δ = r − a
The derivatives of the scalar potential can be written in terms of the electric field, by substituting the relation:
E = −∇φ − ∂A
∂t (2.39)
The scalar potential now be written as (2.40):
φ(r, t) = φ(a, t) − ∞ X n=1 1 n![(δ · ∇ 0)n−1(δ · E(r0, t))] r0=a − ∂ ∂t ∞ X n=1 1 n![(δ · ∇ 0)n−1(δ · A(r0, t))] r0=a . (2.40) Equation (2.40) has the form of a gauge transformation, shown in (2.41), which relates potentials that correspond to identical electromagnetic fields.
φ0→ φ − ∂χ
∂t A
0 → A + ∇χ
(2.41) By substituting these transformed potentials in expression (2.39), it can be shown that the electric field is invariant with respect to this transformation,
E = −∇φ0−∂A 0 ∂t = −∇φ − ∂A ∂t + ∇ ∂χ ∂t − ∂(∇χ) ∂t = −∇φ − ∂A ∂t, (2.42) where the last equality holds, due to the commutativity of partial derivatives. Equivalently, the gauge invariance of the magnetic field can be shown by substituting the transformed vector potential:
B = ∇ × A0= ∇ × A + ∇ × (∇χ) = ∇ × A. (2.43)
Note that the last equality holds, because the curl of any gradient vanishes. Alternatively, gauge trans-formations can be expressed as a gauge condition, which directly constrains the potentials, thus fixing the additional degree of freedom provided by gauge transformation.
Going back to our example in (2.40), the gauge function has the form shown in equation (2.44).
χ(r, t) = ∞ X n=1 1 n![(δ · ∇ 0)n−1(δ · A(r0, t))] r0=a (2.44)
By applying this gauge transformation in reverse, the third term in (2.40) can be diminished, yielding a scalar potential which is completely described in terms of the electric field. The same can be derived for the vector potential by using the relation B = ∇ × A, resulting in an expansion of this potential in terms of the magnetic field. Expression (2.45) contains the final results:
φ(r, t) = φ(a, t) − ∞ X n=0 1 (n + 1)![(δ · ∇ 0)n(δ · E(r0, t))] r0=a (2.45) A(r, t) = − ∞ X n=1 n (n + 1)![(δ · ∇ 0)n−1(δ × B(r0, t)] r0=a .
It is clear at this point that both potentials are expanded in terms of electric and magnetic fields, which on their own are gauge-invariant quantities. The gauge freedom now resides in the choice of expansion point a.36
The potentials in (2.45) are in a general form, that is the fields still have to be specified. In our case, the potentials need to describe the perturbation by linearly polarized light, depicted in (2.17),
φ(r, t) = −E 2(δ · L) ∞ X n=0 1 (n + 1)!(i
n−1(δ · k)nei(k·a−ωt+δ)+ (−i)n−1(δ · k)ne−i(k·a−ωt+δ)
) (2.46) A(r, t) = − E 2ω(δ × (k × L)) ∞ X n=1 n (n + 1)!(i
n−2(δ · k)n−1ei(k·a−ωt+δ)+ (−in−2)(δ · k)n−1e−i(k·a−ωt+δ)
).
By plugging these potentials in the linear interaction term from (2.8), one obtains,
ˆ V (t) = −(ecE 2ω α · (δ × (k × L)) ∞ X n=0 n + 1 (n + 2)!(i n−1(δ · k)n+ eE 2 (δ · L) 1 (n + 1)!i n−1(δ · k)n)ei(k·a−ωt+δ) (2.47) n + 1 (n + 2)!((−i) n−1(δ · k)n+eE 2 (δ · L) 1 (n + 1)!(−i) n−1(δ · k)n)e−i(k·a−ωt+δ) ,
which exactly has the right form to extract the frequency component of ˆV (t). After setting the expansion point, a, to zero, the effective operator, relabeled as ˆTmpol(ω), can now be written as,
ˆ
Tmpol(ω) = −ieiδQ[1]p L;p+ eiδ ∞ X n=1 in−1L;pkj1kj2· · · kjnXˆj1···jn;p; (2.48) ˆ X[n]= 1 (n + 1)!Qˆ [n+1] j1···jn,p− i ω 1 n!mˆ [n] j1···jn−1;rrjnp ˆ Q[n]j 1···jn = −er1r2· · · rjn; mˆ [n] j1···jn−1;jn= n n + 1rj1rj2· · · rjn−1(r × ˆj)jn; ˆj = −ecα, where appears the electric, ˆQ[n]j
1···jn, and magnetic, ˆm
[n]
j1···jn−1;jn, multipole operators. Note how the
with respect to the large and small components, whereas this is not the case for the electric multipole operator.
The effective interaction operator is not Hermitian, although this property can be imposed in the limit, k · r 1, provided that the proper value of δ is chosen. In this limit, ˆTmpol(ω) becomes:
ˆ
Tmpol(ω) = −ieiδQ[1]p L;p+eiδ ∞ X n=1 in−1L;pkj1kj2· · · kjnXˆj1···jn;p −−−−→ k·r 1 ˆ
Tmpol(ω) = −ieiδQ[1]p L;p.
(2.49) Therefore, the phase factor, δ, has to be set to π2, to ensure Hermiticity of the dipolar term.
In order to find the source of the diverging transition moments, it is desirable to simplify (2.48). This can be achieved by taking the non-relativistic limit, i.e. c → ∞. As pointed out in section 2.1.1, magnetism vanishes in the non-relativistic limit, hence the interaction operator becomes:
ˆ Tmpol(ω) = Q[1]p p+ ∞ X n=1 inpkj1kj2· · · kjnQˆ [n+1] j1···jn,p. (2.50)
2.2.4 Truncated Interaction Operator: Generalized Velocity Representation
In the generalized velocity representation, the truncated operator is constructed by Taylor-expanding the full-interaction operator with respect to the wave vector:
ˆ Tf ull= ∞ X n=0 kn n! dn dk0n[ e ω(cα · )e ik0·r ]|k0=0 = ∞ X n=0 e ω in n!(cα · )(k · r) n. (2.51)
This simple formulation becomes equivalent to the generalized length representation in the limit of a complete basis and both representations are related by a gauge transformation. For a given order, n, the term, (k · r)n, generates all possible multipole moments of this order by means of the trinomial expansion: (k · r)n= (kxx + kyy + kzz)n= X i,j,k i+j+k=n n i, j, k (kxx)i(kyy)j(kzz)k (2.52) n i, j, k = n! i!j!k!.
The most prominent difference between the generalized length and velocity representation is their gauge-origin dependency at finite truncation order: the former manifests itself through the choice of expansion point, whereas the latter introduces lower-ordered multipole moments upon variation of the gauge-origin. Furthermore, in the generalized length representation, the truncated interaction operator is split into a magnetic and electric contribution, while the generalized velocity representation does not distinguish these two components.
2.3 Absorption Cross Section
2.3.1 Full Interaction Operator: Oriented Absorption Cross Section
The absorption cross section is the observable quantity being measured and it is computed by taking the square modulus of transition moments derived from the effective interaction operator.
πω ~
|hf | ˆT (ωf i) |ii |2 (2.53)
For the derivations in this section, the lineshape function, f (ω, ωf i, γf i), is left out of the absorption
cross section. In exact state theories, the transition moments are calculated as off-diagonal elements of the effective interaction operator, whereas in Hartree-Fock theory, the transition moments are calculated by contracting the solution vectors and the gradient , the latter defined in (2.54):
[ET[1]]ai= − hφa| ˆT (ω) |φii . (2.54)
In this equation, {|φni} denotes an orthonormal basis of orbitals and ET[1] the gradient of the effective
interaction operator ˆT (ω).39 The gradient is a matrix stored in vector form, providing information on
the influences of the perturbing operator, with its indices ai representing a superindex, indicating the component of the vector. The solution vectors on the other hand contain information on the amplitudes of the excited states. This method for the calculation of transition moments is called response theory and a complete account of this theory is beyond the scope of this section, thus more information is given in appendix A.
For operators that are decomposed in Hermitian and anti-Hermitian parts, as in (2.36) and (2.35), the transition moments can be calculated as,
hf | ˆTf ull|ii = −i hf | i ˆTH|ii + hf | ˆTA|ii = −iE †[1]
iTHX+;n+ E
†[1]
TAX+;n, (2.55) where appears the solution vector X+;n. The interaction operators in (2.55) are symmetric with respect
to time-reversal symmetry, hence their transition moments are real. Equation (2.55) has the form a + ib, with a and b representing real numbers, thus its square modules does not contain terms of mixed Hermiticity. The absorption cross section can be computed by taking the square modulus of (2.55),
σ(ω) = πω 0~c | hf | ˆT (ωf i) |ii |2 = πω 0~c(| hf | i ˆTH(ωf i) |ii | 2+ | hf | ˆT A(ωf i) |ii |2). (2.56)
Evaluated for linearly polarized light, each of these terms become:
σH(ω) =
πe2 0~ωc
| hf | icα · Lcos(k · r) |ii |2 (2.57)
σA(ω) =
πe2 0~ωc
and for circularly polarized light:
σH(ω) =
πe2
20~ωc(| hf | icα · 1cos(k · r) |ii |
2+ | hf | icα ·
2sin(k · r) |ii |2 (2.58)
− (−1)m2 Re[hf | icα · 2sin(k · r) |ii hf | icα · 1cos(k · r) |ii∗])
σA(ω) =
πe2
20~ωc(| hf | icα · 1sin(k · r) |ii |
2+ | hf | icα ·
2cos(k · r) |ii |2
+ (−1)m2 Re[hf | icα · 1sin(k · r) |ii hf | icα · 2cos(k · r) |ii∗])
It should be noted that these transition moments are real, thus in expression 2.58, it is trivial to take the real part of the m-dependent term. The results of this expression can be used to compute the electronic circular dichroism, which is defined as the difference between absorption of left-handed and right-handed light:36
∆σ = σL− σR. (2.59)
Consider the structure in equations (2.58). By taking the difference between the absorption of left and right handed polarized light, all terms that do not involve the parameter m vanish, recalling that m is odd for right handed polarization and even for left handed polarization:
∆σ = 2πe
2
0~ωc(hf | icα · 1sin(k · r) |ii hf | icα · 2cos(k · r) |ii
∗
(2.60) − hf | icα · 2sin(k · r) |ii hf | icα · 1cos(k · r) |ii∗).
By comparing the absorption cross section of linearly polarized light in (2.57) with the ECD from ex-pression (2.60), it becomes clear that linear absorption is described by one polarization vector and there is no mixing between sine and cosine terms, whereas ECD only has mixed terms and it is constructed from two unit vectors.
2.3.2 Full Interaction Operator: Isotropic Absorption Cross Section
Until this point it was assumed that both the orientation of the molecule and the electromagnetic wave are fixed in space. The physical picture of this model corresponds to an ensemble of non-interacting molecules that are oriented in the same direction. In typical gaseous or liquid samples, the absorption is measured as a statistical average over all molecular orientations, so it suffices to take the isotropic average of the property under consideration.28In three dimensions, rotational averaging has to be carried out over the Euler angles, χ, φ and θ. For a generic function in this space, f (r), the rotational average is calculated by evaluating the integral in (2.61),
hf (r)iθ,φ,χ= R2π 0 R2π 0 Rπ 0 dθdφdχ sin θf (r) R2π 0 R2π 0 Rπ 0 dθdφdχ sin θ = 1 8π2 Z 2π 0 Z 2π 0 Z π 0 dθdφdχ sin θf (r). (2.61)
To perform the rotational averaging, two frames are needed to define the rotation. In the first frame S, the electromagnetic wave is at rest, while in the second frame, S0, the molecule is at rest.28In spherical coordinates, the transformation S → S0is described by,
er = ex0sin θ cos φ + ey0sin θ sin φ + ez0cos θ (2.62) eθ = ex0cos θ cos φ + ey0cos θ sin φ − ez0sin θ
eφ= −ex0sin φ + ey0cos φ
where θ and φ represent the rotation angles, while (ex0, ey0, ez0) and (er, eθ, eφ) denote the unit vectors in frame S0 and S, respectively. For convenience, the radial unit vector will be aligned with the wave vector, thus all dependency on χ is contained in the polarization vectors.
ek= er; L= 1= cos χeθ+ sin χeφ 2 = cos χeφ− sin χeθ. (2.63)
For linear absorption, the rotational averaging yields,
hσ(ω)iθ,φ,χ = πe
2
ω0~c
hhf | icαpL;peik·r|ii hf | icαqL;qeik·r|ii∗iθ,φ,χ (2.64)
= πe
2
ω0~c
hhL;pL;qiχhhf | icαpeik·r|ii hf | icαqeik·r|ii∗iθ,φ,
and for ECD, the rotational averaging becomes:
h∆σ(ω)iθ,φ,χ = h
2πe2
0~ωc(hf | icα · 1sin(k · r) |ii hf | icα · 2cos(k · r) |ii
∗
(2.65) − hf | icα · 2sin(k · r) |ii hf | icα · 1cos(k · r) |ii∗)iθ,φ,χ
= 2πe
2
0~ωc(h1;p2;q
iχhhf | icαpsin(k · r) |ii hf | icαqcos(k · r) |ii∗iθ,φ
− h2;p1;qiχhhf | icαpsin(k · r) |ii hf | icαqcos(k · r) |ii∗iθ,φ).
In both equations, the Einstein summation convention is used. The averaging over χ can be performed analytically in both cases, by evaluating the averaging over the polarization vectors, starting with linear absorption:
hL;pL;qiχ =
1
2(δpq− er;per;q). (2.66)
By plugging this expression in (2.64), one obtains:
hσ(ω)iθ,φ,χ = πe
2
2ω0~c(δpq
− er;per;q)hhf | icαpeik·r|ii hf | icαqeik·r|ii∗iθ,φ (2.67)
h1;p2;qiχ =
1
2(eθ;peφ;q− eφ;peθ;q). (2.68) Equation (2.68) can be greatly simplified by identifying the right hand side of the equation with the relation 1× 2 = ˆk = er. Therefore, the rotational averaging over chi becomes:
h1;p2;qiχ=
1
2ipqer;i. (2.69)
Note that this expression is anti-symmetric with respect to exchange of the indices p and q. Therefore, equation (2.69) can also be used to evaluate h2;p1;qiχ. Plugging these expressions in (2.65):
h∆σ(ω)iθ,φ,χ= 2πe
2
0~ωcipqer;i
hhf | icαpsin(k · r) |ii hf | icαqcos(k · r) |ii∗iθ,φ. (2.70)
The averaging over the angles θ and φ needs to be done numerically, which was carried out with Lebedev quadrature for our purposes.40–44
2.3.3 Truncated Interaction Operator: Generalized Length Representation
Like previous derivations, the absorption cross section is obtained by plugging the effective operator, in this case (2.48), in the absorption cross section, (2.14).
σmpol(ω) = πω 0~c ∞ X n=0 ∞ X m=0 hf | (L;pinkj1kj2· · · kjnXˆ [n] j1···jn;p(ω) |ii×hf | (L;qi mk j1kj2· · · kjmXˆ [m] j1···jm;q(ω) |ii ∗ (2.71) After gathering the terms of the same power in k and relabelling the summation indices, the contributions to σmpol(ω) at order n become:
σ[n]mpol(ω) = πω 0~cL;pL;q × n X m=0 in(−1)n−mkj1kj2· · · kjnhf | ˆX [m] j1···jm;p(ω) |ii hf | ˆX [n−m] jm+1···jn;q(ω) |ii ∗ . (2.72) In expression (2.72), the index m increases until it reaches the value n, while the index n − m decreases until it reaches zero. Therefore, every term in this summation is not unique, because for every element in the summation, there is an equivalent element that is related to the former through complex conjugation:
hf | ˆXj[m]1···jm;p(ω) |ii hf | ˆXj[n−m]m+1···jn;q(ω) |ii∗ ∗→ hf | ˆ− Xj[n−m]1···jm;p(ω) |ii hf | ˆXj[m]m+1···jn;q(ω) |ii∗. (2.73)
However, there is an exception to the above mentioned rule, as even-ordered values of expression (2.72) contain one term of the form:
hf | ˆXj[0] 1···jm;p(ω) |ii hf | ˆX [0] jm+1···jn;q(ω) |ii ∗ . (2.74)
Because the sum is alternating, the complex conjugates are added at even orders, while they are subtracted at odd orders. Due to this property, the well-known identity can be used, that relates the components of a complex number to this number itself and its conjugate:
Re (Z) = 1 2(Z + Z ∗ ); Im (Z) = 1 2i(Z − Z ∗ ). (2.75)
Therefore, the even and odd contributions to the absorption cross section can be seperated into two distinct terms: σmpol[2n+1](ω) = πω 0~cL;pL;q n+1 X m=1 (−1)mkj1kj2· · · kj2n+12 Im [hf | ˆX [n+m] j1···jn+m;p(ω) |ii hf | ˆX [n+1−m] jn+m+1···j2n+1;q(ω) |ii ∗ ] (2.76) σmpol[2n] (ω) = πω 0~cL;pL;q n X m=0 (−1)m(2 − δm0)kj1kj2· · · kj2n2 Re [hf | ˆX [n+m] j1···jn+m;p(ω) |ii hf | ˆX [n−m] jn+m+1···j2n;q(ω) |ii ∗ ].
The multipole moment operator, ˆX(ω), is time-reversal symmetric, hence its transition moments are real, while the odd-ordered contributions are composed of the imaginary parts of the transition moments, as indicated in appendix A. Therefore, only the even-ordered components of σ[n](ω) contribute to the absorption cross section.
2.3.4 Truncated Interaction Operator: Generalized Velocity Representation
Similar to the generalized length representation, the absorption cross section can be rewritten by collect-ing the terms with the same order of the wave vector.
σ(ω) = πω 0~c ∞ X n=0 n X m=0
[hf | ˆTf ull[n−m](ω) |ii hf | ˆTf ull[n+m](ω) |ii∗] (2.77)
ˆ Tf ull[n] (ω) = k n n! dn dk0n[ e ω(cα · )e ik0·r] k0=0 = e ω in n!(cα · )(k · r) n
For even values of (2.78), the transition moments derived from ˆTvel[n](ω) are anti-symmetric with respect to time-reversal, while it is symmetric for even values, implying that the transition moments are real for odd values and imaginary for even values. The summation terms in (2.77) are not unique, because they come in pairs related by complex conjugation:
hf | ˆTf ull[p] (ω) |ii hf | ˆTf ull[q] (ω) |ii∗ ∗→ hf | ˆ− Tf ull[p] (ω) |ii∗hf | ˆTf ull[q] (ω) |ii . (2.78) Because the summation in (2.77) is not alternating, each pair can be collected as,
2 Re[hf | ˆTf ull[p] (ω) |ii hf | ˆTf ull[q] (ω) |ii∗], (2.79) with the only exception being the special case where p = q. Consider the case where one of the two indices is even and the other odd. In this case, a purely real number will be multiplied with a purely
imaginary number, resulting in an imaginary number. The combination of each pair can be written as (2.79), and these contributions will vanish, because this number is purely imaginary. Compare this to the case where either p and q are both even or odd. The multiplication now results in a real number, hence only the pairs where the sum of p and q is even contribute to the absorption cross section.
σvel[2n](ω) = πω 0~c
n
X
m=0
(2 − δm0) Re [hf | ˆTf ull[n−m]|ii hf | ˆTf ull[n+m]|ii∗] (2.80)
2.4 Orbitals and Basis Sets Suitable for X-ray Transitions
2.4.1 Completeness-Optimized Basis Sets
Conventionally, Gaussian basis sets are constructed for each atom separately by finding the exponents of a fixed number of Gaussians that minimize the energy. In equation (2.81) a generic Gaussian function, Gαijk, is defined,
Gαijk = Nijkα xiyjzke−αr2, (2.81) in which appears the normalization constant, Nijkα , the exponent, α, and the Cartesian powers, i, j, k. Due to the nature of the Coulomb potential, the energy tends to be described best by a local basis set, whereas other properties can have different requirements.45For example, vacant orbitals at high energies tend to be diffuse, indicating that excitations to such orbitals presumably need a basis set containing more diffuse functions than the energy-optimized basis. The divergence of the truncated interaction operator could be an artefact of the energy optimized basis sets. Therefore, completeness-optimized basis sets were used for the calculations with the truncated interaction operator in the generalized length representation, in order to exclude this hypothesis.
A basis is said to be complete if the identity operator can be constructed out of its members, although in practice this can only be achieved approximately. Furthermore, the extent of completeness can be tested by projecting the constructed identity operator onto another Gaussian and verifying how much the norm of this transformed basis function differs from one.45 By plotting the norm against its exponent, α, this function can be made into a completeness profile, which is a curve, generally shaped as a plateau, indicating the extent of completeness by its flatness and width, shown in equation (2.82),
Y (α) =
N
X
b=1
hα| biS−1hb |αi (2.82)
where S−1indicates the inverse metric of the Gaussian test basis, |bi, and |ai the Gaussian basis function on which the approximate identity operator is projected. This idea can be made more elaborate by optimizing the above relation for a range of exponents, yielding the criterion that was proposed by Manninen and Vaara.46By using (2.83), the optimization of basis sets requires the minimization of the integrand in equation (2.83), whereas conventional basis set optimization makes use of the minimization of energy.
τI =
1
log(αmax) − log(αmin)
Z log(αmax)
log(αmin)
[1 − Y (α)]d log(α) (2.83)
In equation (2.83), αmaxand αminrepresent the upper and lower bound of the exponent range, respec-tively.
2.4.2 Improved Virtual Orbital
All calculations in this report were performed using the Hartree-Fock method, which is a variational approximation that describes the wave function as a single determinant that minimizes its corresponding energy with respect to the orbitals, while they are constrained to remain orthonormal.47To analyze the results of the basis set series, it was convenient to focus on one specific excitation. In DIRAC, excitations can be selected by specifying the indices that correspond to the orbitals involved in the excitation. The indices of the orbitals are sorted from low to high energy and it was found that the energies of the virtuals would vary significantly, resulting in a different order of these indices. Therefore, inequivalent excitations were compared, because the excitations were identified using these indices.
At the root of this problem lies the summation index of the Hartree-Fock orbital energies, which only runs over the occupied orbitals, thus overcounting a repulsion term that would vanish if the orbital is occupied.48 This concept is illustrated in the equation shown below: for an occupied orbital, j, the summation of repulsion energies contains N − 1 non-vanishing terms. One term vanishes, due to can-cellation of the Coulomb and exchange term. Virtual orbitals on the other hand, have N repulsion terms, each of which is non-vanishing. Due to this additional interaction, virtual orbital energies from a system of N electrons in fact describe the energy of an orbital contained in a system of N + 1 electrons.
a= hi| h |ii + N X j6=i (Jij − Kij) (Jjj− Kjj) = 0 (Jja− Kja) 6= 0 (2.84)
In equation (2.84), K represents the exchange interaction, J the Coulomb interaction, a the energy
of a virtual orbital and j and a represent the indices for occupied and virtual orbitals, respectively. Excitations to low-lying virtual levels can be described with decent accuracy with the usual Hartree-Fock virtuals. For higher energy levels, the virtuals are artificially increased in energy, due to the additional repulsion energy of these virtual orbital with the occupied orbitals. As a result, these orbitals become disproportionately high in energy. At high photon energies the electrons are excited to Rydberg states instead, although physically the excitations should occur to the orbitals that are artificially increased in energy. For basis set tests this can cause problems, as the energy levels of the Rydberg states shuffle significantly, thereby changing the nature of the excitation across the basis set series. An improved description of the virtual orbitals was used for the basis set test for the truncated interaction operator in the generalized length representation, in order to make the excitations unambigously accros the series.
The description of the virtuals was improved by applying the scheme described by Hunt and co-workers.48First, the neutral system is optimized and the MO coefficients are used as input for the next calculation. Secondly, these MO coefficients are used in an optimization of the vacant orbital involved in
the excitation, while the N-1 occupied orbitals are kept frozen, assuming they do not change significantly upon excitation. This step is required to yield a proper basis for the virtual orbital under consideration. Finally, the generated basis for the virtual orbital is used in an TD-HF calculation, resulting in the correct description of the excited states.
3
Computational Details
3.1 Truncated Interaction Operator
3.1.1 Generalized Length Representation
In order to assess the extent of the divergence of the truncated interaction operator in the generalized length representation, several excitations were calculated for the earth alkaline metals ranging from magnesium to radium. This series of calculations was performed non-relativistically and isotropically at the Hartree-Fock level of theory, using the DIRAC program.32This level of theory was chosen, to isolate the origin of the divergence, as DFT introduces numerical integration, unnecessarily complicating the analysis. For similar reasons, a non-relativistic Hamiltonian was used for these calculations, implying that the interaction operator reduces to (2.50).
It could be concluded that the divergences started to appear in the calculations of strontium. For this system, the oscillator strengths at higher orders would approximately stay in the same order of magnitude. the calculations with barium on the other hand, yielded oscillator strengths that diverged more substantially. Therefore, the oscillator strengths of barium were studied more thoroughly. This system was analyzed by calculating the 5px ← 1s transition at 8th truncation order. This specific
transition was selected, because it is the lowest transition of the 1s electron.
Hypothetically, the divergences of the oscillator strenghts can be explained by incorrect values of the property matrix. For barium, the elements of this matrix were tested by computing the integrals of (2.50) analytically with Mathematica and comparing these numbers with the output of DIRAC.49In an unnormalized Gaussian basis, matrix elements of this operator have the form,
Z d3rxiyjzke−ar2 = 1 8(1+(−1) i)(1+(−1)j)(1+(−1)k)a−3+i+j+k2 Γ 1 + i 2 Γ 1 + j 2 Γ 1 + k 2 , (3.1) where appears the gamma function, defined as:
Γ(z) = Z ∞
0
dttz−1e−t. (3.2)
The derivation of expression (3.1) is given in appendix B. Due to these simple analytical expression, the output of DIRAC could be tested easily for incorrect values.
An unbalanced description of the basis sets could be another explanation for the divergences. It was pointed out by Sørensen et al. that unbalanced basis sets yield negative oscillator strengths, indicating that these calculations are sensitive with respect to the choice of basis sets.26 For these calculations, completeness optimized basis sets were used, in order to exclude the possibility of artefacts introduced by energy-optimized basis sets. These basis sets were composed with Erkale, which is a software pack-age designed for the calculations of X-ray properties and the fast generation of completeness-optimized basis sets consisting of Gaussians.50This program takes the the number of primitives and the boundaries of the exponent range as input for each angular momentum shell separately, requiring the output of these
calculations to be concatenated manually. Completing these steps returns the completeness profiles for each shell and the list of exponents of the completeness-optimized basis.
A suitable basis set was found for these calculations by completeness optimizing the s and p shell of the dyall.ae2z basis set, while keeping the higher shells unaltered.51 Additionally, the s shell was augmented with either tight or diffuse functions, ranging from zero to three additional functions. For the s shell, the values of αmaxand αminwere found by increasing or decreasing the existing maximum or minimum of the dyall.ae2z basis set in an even-tempered fashion. Subsequently, these boundary values were being fed to Erkale, yielding the augmented s shell and the amount of exponents was increased accordingly.
Further basis set testing was performed by increasing the total size of the basis sets in a calcula-tions of the isotropic oscillator strength. This series of calculacalcula-tions was performed relativistically at the Hartree-Fock level of theory, using the dyall.ae2z, dyall.ae3z and dyall.ae4z basis sets.51For these calculations, the excitation 5p1
2 ← 1s is studied, which is the relativistic analogue of the 5px ← 1s transition. By performing this series of calculations, it could be assessed what the influence is of the basis set size. Furthermore, this series also serves as a comparison with the generalized velocity repre-sentation, for which equivalent calculations were performed. For these reasons, this series ought to be calculated relativistically, as there is no straightforward non-relativistic limit of (2.51).
As indicated in section 2.4.2, the high energies of X-ray transitions, require the use of improved virtual orbitals. Therefore, these orbitals were used to make sure no artefacts could have been caused by improper virtual orbitals. Furthermore, it was found that convergence of the solution vectors could only be reached if the criterion is based on the energy, with a threshold of 10−8, but conventional calculations base their criterion on the norm of the error vector used in the DIIS procedure.
In order to test whether the completeness-optimized basis sets are balanced, radial distributions were made from the transition moment densities for the contribution at 8thorder. Based on its behaviour upon augmentation, these distributions can serve as an indication of basis set requirements.
Radial distributions are functions of the radial distance and they are constructed by integrating out the angular dependency of some three dimensional function, thus for spherically symmetric systems, these distributions have a clear interpretation. In equation (3.3), radial distributions are defined,
P (r) = Z π 0 Z 2π 0 dθdφr2sin θf (r), (3.3)
where f (r) denotes some generic three dimensional scalar function. In our application, the anisotropic transition moment densities were used for these distributions, thus a frame of reference is required. In this frame, the polarization and wave vector were chosen such that former is aligned with the x-axis and the latter with the z-axis. With this configuration and at this order (2.50) becomes:
ˆ
Tmpol[8] = k8X8Z. (3.4)
3.1.2 Generalized Velocity Representation In the generalized velocity respresentation, the 5p1
2 ← 1s was studied by calculating this excitation relativistically at the Hartree-Fock level of theory for various basis set sizes, as was indicated in section 3.1.1. The used basis sets were dyall.ae2z, dyall.ae3z and dyall.ae4z and the results of this series are compared to equivalent calculations in the generalized length representation.51
Additionally, the truncated interaction operator in the generalized velocity representation was ap-plied in the calculation of oscillator strengths of a radium atom. The aim of these calculations was to find for which energy scale the divergences are introduced, thus it was desirable to perform these cal-culations with radium, as the divergences of this system are more pronounced. All studied excitations were of the type, 7p1
2 ← ns, where n denotes all occupied s orbitals of radium, i.e. 1, 2, · · · , 7 and the interaction operator was truncated at 12thorder. Furthermore, only the excitations of B1usymmetry
were studied to keep the analysis of the results straightforward and a frame was chosen where the polar-ization vector is aligned with the z-axis and the wave vector with the y-axis. In this frame and with B1u
symmetry, the interaction operator, (2.51), at order n becomes:
ˆ Tvel2n(ω) = e ω i2n (2n)!cαz ω c 2n y2n (3.5)
The purpose of these calculations was to test the convergence properties of the truncated interaction operator. Because the generalized velocity representation, (2.51), is derived from a Taylor expansion of the full interaction operator, (2.31), the former should converge to the latter in the limit n → ∞. The convergence behaviour was studied by tabulating the transition moments and the integrals of the truncated interaction operator. In order to evaluate the property integrals, basis functions are required, which were taken from the most dominant contribution of the gradient that describes the 7p1
2 ← 1s transition. By analysing the integrals and transition moments, it could be deduced at which step of the calculation the divergences are introduced.
Additionally, the convergence behaviour of the transition moments, (3.6), was studied by artificially varying ω in the vicinity of ω = c and plotting the result as function of ω for each order separately.
hf | ˆTvel[2m](ω) |ii = m X n=0 e ω i2n (2n)!hf | cαz ω c 2n y2n|ii (3.6)
In expression (3.6), square brackets denote accumulated quantities.
It was found that the largest elements of the gradient consist of py functions. For our frame of
reference and specific basis functions, the integrals can be expressed as:
ec ω i2n (2n)!k 2nhGα1 010| y 2n|Gα1 010i = ec ω i2n (2n)!k 2nNα1 1 N α2 1 (α1+ α2)− 2n+3 2 Γ 2n + 3 2 (3.7) and hGα1 010| ec ωe iky|Gα1 010i = ec 2ωN α1 1 N α2 1 √ π(α1+ α2)− 3 2e− k2 4(α1+α2) 1 − k 2 2(α1+ α2) , (3.8)
