Real-time solutions of coupled Ehrenfest-Maxwell-Pauli-Kohn-Sham equations: fundamentals, implementation, and nano-optical applications

TL;DR

This paper develops a density-functional approach for simulating coupled photon-electron-nuclei systems in non-relativistic quantum electrodynamics, enabling realistic modeling of light-matter interactions in nano-optics and related fields.

Contribution

It introduces a Kohn-Sham formalism for coupled light-matter systems and implements it in the Octopus code using a Schrödinger form of Maxwell's equations.

Findings

Formalism establishes a one-to-one correspondence between external fields and internal variables.

Implementation demonstrates feasibility of real-time simulations of coupled systems.

Applicable to various nano-optical and light-matter interaction scenarios.

Abstract

We present the theoretical foundations and the implementation details of a density-functional approach for coupled photons, electrons, and effective nuclei in non-relativistic quantum electrodynamics. Starting point of the formalism is a generalization of the Pauli-Fierz field theory for which we establish a one-to-one correspondence between external fields and internal variables. Based on this correspondence, we introduce a Kohn-Sham construction which provides a computationally feasible approach for ab-initio light-matter interactions. In the mean-field limit for the effective nuclei the formalism reduces to coupled Ehrenfest-Maxwell-Pauli-Kohn-Sham equations. We present an implementation of the approach in the real-space real-time code Octopus. For the implementation we use the Riemann-Silberstein formulation of classical electrodynamics and rewrite Maxwell's equations in Schr\"odinger form. This allows us to use existing time-evolution algorithms developed for quantum-mechanical systems also for Maxwell's equations. Overall, our approach is ideally suited for applications in nano-optics, nano-plasmonics, (photo) electrocatalysis, light-matter coupling in 2D materials, cases where laser pulses carry orbital angular momentum, or light-tailored chemical reactions in optical cavities to name but a few.

Tables

Table 1. Table 1: Table of acronyms that are used in the present work to indicate the level of EMPKS theory.

Acronym	Description
F@ED	Forward coupling with Electric Dipole term
FB@ED	Forward and Backward coupling with Electric Dipole term
F@(ED+MD+EQ)	Forward coupling with Electric Dipole, Mag- netic Dipole and Electric Quadrupole term
FB@(ED+MD+EQ)	Forward and Backward coupling with Electric Dipole, Magnetic Dipole and Electric Quadrupole term

Table 2. Table 2: Simulation parameters for the sodium dimer for distances d 1 = 0.1 subscript 𝑑 1 0.1 d_{1}=0.1 nm and d 2 = 0.5 subscript 𝑑 2 0.5 d_{2}=0.5 nm.

	distance $d_1$		distance $d_2$
variable	conv. units	[a.u.]	conv. units	[a.u.]
$\omega$	3.05 eV	0.112	2.83 eV	0.104
$k_x$	1.55 e-11 m-1	8.17 e-4	1.43 e-11 m-1	7.59 e-4
$\lambda$	406.5 nm	7681.84	438.1 nm	8279.02
$E_{0,z}$	5.142 e7 V/m	1.0 e-4	5.142 e7 V/m	1.0 e-4
Intensity	3.51 e12 W/m2	5.45 e-4	3.51 e12 W/m2	5.45 e-4
$\xi$	2034.08 nm	38438.5	2034.08 nm	41395.1
$x_0$	4068.16 nm	76877.0	4381.07 nm	82790.2
$L_{\mathrm{KS},x}$	1.993 nm	37.658	1.993 nm	37.658
$L_{\mathrm{KS},y}$	1.993 nm	37.658	1.993 nm	37.658
$L_{\mathrm{KS},z}$	3.347 nm	63.258	3.547 nm	67.037
$L_{\mathrm{Mx},x}$	2.646 nm	50.000	2.646 nm	50.000
$L_{\mathrm{Mx},y}$	2.646 nm	50.000	2.646 nm	50.000
$L_{\mathrm{Mx},z}$	4.498 nm	85.000	4.498 nm	85.000
$a_{\mathrm{Mx},x}$	2.170 nm	41.000	2.170 nm	41.000
$a_{\mathrm{Mx},y}$	2.170 nm	41.000	2.170 nm	41.000
$a_{\mathrm{Mx},z}$	4.022 nm	76.000	4.022 nm	76.000
$\Delta x_{\mathrm{KS}}$	0.053 nm	1.000	0.053 nm	1.000
$\Delta x_{\mathrm{Mx}}$	0.053 nm	1.000	0.053 nm	1.000
$\Delta t_{\mathrm{KS}}$	5.096 e-3 fs	0.211	5.096 e-3 fs	0.211
$\Delta t_{\mathrm{Mx}}$	1.019 e-4 fs	4.21 e-3	1.019 e-4 fs	4.21 e-3

Equations

$$g_{\mu \nu }=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& -1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& -1\end{array}\right).$$

$$\vec{A}$$

$$\vec{A}\cdot \vec{B}$$

$$\vec{\nabla }\cdot \vec{A}$$

$$A_k=g_{kl}{A}^l$$

$$\vec{\nabla }\times \vec{A}$$

$${\gamma }^0=\left(\begin{array}{cc}{1}_2& 0\\ 0& -1_2\end{array}\right),{\gamma }^k=\left(\begin{array}{cc}0& {\sigma }^k\\ -{\sigma }^k& 0\end{array}\right).$$

$${\sigma }^k{\sigma }^l$$

$$={\delta }^{kl}{1}_2-\mathrm{i}{ϵ}^{klj}{\sigma }_j.$$

$$S^1$$

$$S^2$$

$$S^3$$

$$\left(\frac{E^2}{c_0^2}-p^2\right)\vec{F}$$

$$-\left(\begin{array}{ccc}(p_1{)}^2& p_1{p}_2& p_1{p}_3\\ p_2{p}_1& (p_2{)}^2& p_2{p}_3\\ p_3{p}_1& p_3{p}_2& (p_3{)}^2\end{array}\right)\vec{F}=0,$$

$$\mathrm{i}\hslash {\partial }_0\vec{F}=-\mathrm{i}\hslash {\partial }_k{S}^k\vec{F}$$

$$-\mathrm{i}\hslash \vec{\nabla }\cdot \vec{F}$$

$$F_{\pm }^k=\sqrt{\frac{{ϵ}_0}{2}}\left(E^k\pm \mathrm{i}{c}_0{B}^k\right),$$

$$\vec{\nabla }\times \vec{E}$$

$$\vec{\nabla }\times \vec{B}$$

$$\vec{\nabla }\cdot \vec{E}$$

$$\vec{\nabla }\cdot \vec{B}$$

$$\left({\partial }_0^2+{\partial }_l{\partial }^l\right)A^k=0.$$

$$\hat{h}=\frac{1}{2M}{\left(-\mathrm{i}\hslash \vec{\nabla }-\frac{q}{c_0}\vec{a}(\vec{r},t)\right)}^2+q{a}^0(\vec{r},t)-\frac{q\hslash }{2M}\vec{S}\cdot \vec{b}(\vec{r},t).$$

$${\left[\hat{\Phi }(\vec{r},s),{\hat{\Phi }}^{†}(\vec{r}{}^{\prime },s{}^{\prime })\right]}_{\pm }={\delta }_{ss{}^{\prime }}{\delta }^3(\vec{r}-\vec{r}{}^{\prime }),$$

$${\hat{\Phi }}^{†}\hat{\Phi }$$

$${\hat{\Phi }}^{†}{S}_k\hat{\Phi }$$

$$\hat{H}$$

$$=-\int {\mathrm{d}}^{3}r\frac{1}{2M}{\hat{\Phi }}^{†}(-\mathrm{i}\hslash {\partial }_k+\frac{q}{c_0}{a}_k)(-\mathrm{i}\hslash {\partial }^k+\frac{q}{c_0}{a}^k)\hat{\Phi }$$

$$+\int {\mathrm{d}}^{3}rq{a}^0{\hat{\Phi }}^{†}\hat{\Phi }-\int {\mathrm{d}}^{3}r\frac{q\hslash }{2M}{\hat{\Phi }}^{†}{S}_k\hat{\Phi }\left(\frac{1}{c_0}{ϵ}^{klm}{\partial }_l{a}_m\right).$$

$$\left\{{\hat{\Phi }}_{(n)};{\hat{\Phi }}_{(n)}^{†};M_{(n)};q_{(n)};S_{(n)}\right\}$$

$${\hat{H}}^{(0)}$$

$$=\sum _{n=1}^N-\frac{1}{2M_{(n)}}\int {\mathrm{d}}^{3}r{\hat{\Phi }}_{(n)}^{†}\left(-\mathrm{i}\hslash {\partial }_k+\frac{q_{(n)}}{c_0}{a}_k\right)\left(-\mathrm{i}\hslash {\partial }^k+\frac{q_{(n)}}{c_0}{a}^k\right){\hat{\Phi }}_{(n)}$$

$$+\sum _{n=1}^N\int {\mathrm{d}}^{3}r{q}_{(n)}{a}^0{\hat{\Phi }}_{(n)}^{†}{\hat{\Phi }}_{(n)}$$

$$+\sum _{n=1}^N-\frac{1}{2M_{(n)}}\int {\mathrm{d}}^{3}r{\hat{\Phi }}_{(n)}^{†}{S}_k^{(n)}\hat{\Phi }\left(\frac{1}{c_0}{ϵ}^{klm}{\partial }_l{a}_m\right).$$

$$\left[{\hat{A}}_k(\vec{r});{ϵ}_0{\hat{E}}_l^{\perp }(\vec{r}{}^{\prime })\right]=-\mathrm{i}\hslash c_0{\delta }_{kl}^{\perp }(\vec{r}-\vec{r}{}^{\prime }),$$

$${\delta }_{ij}^{\perp }(\vec{r}-\vec{r}{}^{\prime })=\left({\delta }_{ij}{\partial }_i\frac{1}{\Delta }{\partial }_j\right){\delta }^3(\vec{r}-\vec{r}{}^{\prime }).$$

$${\hat{A}}^k(\vec{r})$$

$$\sqrt{\frac{\hslash c_0^2}{{ϵ}_0(2\pi {)}^3}}\int \frac{{\mathrm{d}}^{3}k}{\sqrt{2{\omega }_k}}\sum _{\lambda =1}^2\overrightarrow{ϵ}(\vec{k},\lambda )\left[\hat{a}(\vec{k},\lambda ){\mathrm{e}}^{\mathrm{i}\vec{k}\cdot \vec{r}}+{\hat{a}}^{†}(\vec{k},\lambda ){\mathrm{e}}^{-\mathrm{i}\vec{k}\cdot \vec{r}}\right],$$

$${\hat{E}}_{\perp }^k(\vec{r})$$

$$\sqrt{\frac{\hslash c_0^2}{{ϵ}_0(2\pi {)}^3}}\int \frac{{\mathrm{d}}^{3}k\mathrm{i}{\omega }_k}{\sqrt{2{\omega }_k}}\sum _{\lambda =1}^2\overrightarrow{ϵ}(\vec{k},\lambda )\left[\hat{a}(\vec{k},\lambda ){\mathrm{e}}^{\mathrm{i}\vec{k}\cdot \vec{r}}-{\hat{a}}^{†}(\vec{k},\lambda ){\mathrm{e}}^{-\mathrm{i}\vec{k}\cdot \vec{r}}\right],$$

Full text

We present an implementation of the approach in the real-space real-time code Octopus. For the implementation we use the Riemann-Silberstein formulation of classical electrodynamics and rewrite Maxwell’s equations in Schrödinger form. This allows us to use existing time-evolution algorithms developed for quantum-mechanical systems also for Maxwell’s equations. We introduce a predictor-corrector scheme and show how to couple the Riemann-Silberstein time-evolution of the electromagnetic fields self-consistently to the time-evolution of the electrons and nuclei. Furthermore, the Riemann-Silberstein approach allows to seamlessly combine macroscopic dielectric media with a microscopic coupling to matter currents. For an efficient absorption of outgoing electromagnetic waves, we present a perfectly matched layer for the Riemann-Silberstein vector. We introduce the concept of electromagnetic detectors, which allow to measure outgoing radiation in the far field and provide a direct way to record various spectroscopies. We present a multi-scale approach in space and time which allows to deal with the different length-scales of light and matter for a multitude of applications. We apply the approach to laser-induced plasmon excitation in a nanoplasmonic dimer system. We find that the self-consistent coupling of light and matter leads to significant local field effects which can not be captured with the conventional light-matter forward coupling. For our nanoplasmonic example we show that the self-consistent foward-backward coupling leads to changes in observables which are larger than the difference between local density and gradient corrected approximations for the exchange correlation functional. In addition, in our example we observe harmonic generation which appears only beyond the dipole approximation and can be directly observed in the outgoing electromagnetic waves on the simulation grid. The self-consistent coupling of the electromagnetic fields to the ion motion reveals significant energy transfer from the electromagnetic fields to matter on the scale of a few tens of femtoseconds.

Overall, our approach is ideally suited for applications in nano-optics, nano-plasmonics, (photo) electrocatalysis, light-matter coupling in 2D materials, cases where laser pulses carry orbital angular momentum, or light-tailored chemical reactions in optical cavities to name but a few.

1. I Introduction

2. II Fundamentals

3. II.1 Relativistic wave equations

4. II.2 Free matter Hamiltonians

5. II.3 Free photon Hamiltonian

6. II.4 Interaction Hamiltonians

7. II.4.1 Longitudinal interactions

8. II.4.2 Transverse interactions

9. III QEDFT for multi-species

10. IV Maxwell-Pauli-Kohn-Sham coupling

11. IV.1 Classical limit for Nuclei

12. IV.2 MPKS-Outlook

13. V Real-space implementation

14. V.1 Maxwell’s equations in Riemann-Silberstein representation

15. V.1.1 Maxwell’s equations in Schrödinger form

16. V.1.2 Current densities and integral kernels

17. V.1.3 Helicity coupling in a linear medium

18. V.2 Time-evolution operator for the Riemann-Silberstein vector

19. V.2.1 Forward-backward coupling and time-reversal symmetry

20. V.3 Maxwell boundary conditions

21. V.3.1 Boundary regions

22. V.3.2 Absorbing boundaries by mask function

23. V.3.3 Absorbing boundaries by perfectly matched layer

24. V.3.4 Incident waves boundaries

25. V.3.5 Incident waves boundaries plus absorbing boundaries

26. V.4 Time-evolution operator for Kohn-Sham orbitals

27. V.5 Full minimal coupling and multipole expansion

28. V.6 Multi-scale implementation

29. V.6.1 Multi-grid types

30. V.6.2 Multi-scale in time

31. V.6.3 Parallelization strategies

32. V.7 Predictor-corrector scheme

33. V.7.1 Forward coupling

34. V.7.2 Forward and backward coupling

35. V.8 Validation and comparison to FDTD

36. VI Applications

37. VI.1 System and simulation parameters

38. VI.1.1 Na297-dimer geometry

39. VI.1.2 Simulation boxes and grid alignment

40. VI.1.3 Measurement regions

41. VI.1.4 Laser pulse shape

42. VI.1.5 Propagators

43. VI.2 Results from EMPKS simulations

44. VI.2.1 Electric field enhancement

45. VI.2.2 Next order in multipole coupling and energies

46. VI.2.3 Electromagnetic detectors and harmonic generation

47. VI.2.4 Ion motion

48. VI.2.5 Comparison of different density functionals

49. VII Summary and conclusions

50. VIII Outlook

I Introduction

Low-energy quantum physics has been divided traditionally into different subfields, e.g., quantum chemistry, quantum optics or solid-state physics. Each of this subfield focuses on a specific part of coupled light-matter systems. Roughly speaking (see Ref. Ruggenthaler et al. (2018) for details), one either prescribes how the electromagnetic field looks like and then determines properties of the matter subsystem, e.g., in quantum chemistry or solid-state physics, or one prescribes the properties of matter and then determines how the photon subsystem behaves, as done in, e.g., quantum optics or photonics. This division is reflected also in the available theoretical methodologies, which either focus on the matter degrees of freedom (see, e.g., Refs Szabo and Ostlund (2012); Martin et al. (2016); Schollwöck (2011); Engel and Dreizler (2011)) or on the electromagnetic field (see, e.g., Refs Loudon (1988); Grynberg et al. (2010); Born and Wolf (2013)). This rough division is further differentiated depending on which part of the matter or photon degrees of freedom are investigated, e.g., nuclear dynamics that drive chemical reactions Tavernelli (2015).

However, there have been a lot of recent experimental results that question these traditional distinctions. For instance, when matter and light couple strongly and neither can be considered a perturbation of the other, then the properties of matter and light can be strongly modified and novel states of matter emerge, such as polaritons (light-matter hybrid states). This happens for single molecules in nanocavities Chikkaraddy et al. (2016) or microcavities Wang et al. (2017), where the confined light modes interact strongly with the matter degrees of freedom. But also in other situations, such as at interfaces or nanostructures Sukharev and Nitzan (2017), strong coupling can change well-established results such as the usual selection rules of quantum chemistry Yamamoto et al. (2014). Even for rather bad cavities and ambient conditions strong coupling can be achieved by, for example, increasing the number of molecules or atoms Ebbesen (2016), which in turn provides a novel and very robust tool to influence and control chemical properties. It has been observed that by merely coupling to the changed vacuum of the electromagnetic field, chemical reactions can be modified Hutchison et al. (2012), well-established limits for energy transfer can be broken Zhong et al. , or Raman processes can be enhanced Shalabney et al. (2015). Strong coupling has been achieved for many different physical systems, e.g., even for living bacteria Coles et al. (2014), and it can be used to engineer novel states of matter such as polariton condensates Byrnes et al. (2014). Besides these strong coupling situations many more cases are known where light and matter become equally important, such as in the case of screening, polarization and retardation effects as observed, e.g., in the energy transfer induced by attosecond laser pulses Sommer et al. (2016) or more traditionally in optical responses Maki et al. (1991). Furthermore there are many situations where usually neglected properties of the light field lead to substantial changes in the matter system, such as in strong-field physics Ludwig et al. (2014), when photons carry a large angular momentum Schmiegelow et al. (2016); Bialynicki-Birula and Bialynicka-Birula (2018); Yue et al. (2018), or when the emission spectrum is investigated in detail Taminiau et al. (2012). Indeed, considering the photon and matter degrees of freedom at the same time can lead to impressive novel practical applications such as daytime radiative cooling Fan (2017).

In the above examples the complex interplay between the basic constituents of coupled light-matter systems – electrons, nuclei and photons – are essential. In most theoretical treatments, however, a strong reduction to only a few important degrees of freedom is performed a priori Ruggenthaler et al. (2018), such as in electronic-structure theory Szabo and Ostlund (2012), where the nuclei and the photons are treated only as external perturbations. While many advanced methods to solve the resulting many-electron equation exist Booth et al. (2009); Orus (2014); Schollwöck (2005); Kotliar et al. (2006); Gull et al. (2011); Metzner et al. (2012); Aoki et al. (2014), they miss most effects of the correlation with nuclei and photons. Furthermore, many approaches have been developed that try to tackle the full electron-nucleus problem, such as exact-factorization approaches Abedi et al. (2010, 2012) or trajectory formulations Kapral (2015); Tully (2012); Miller (2012). However, what the correlation with the light-field is concerned, there are only a few such approaches available to date. On the one hand we have coupling with classical light fields, such as presented in Lorin et al. (2007); Lopata and Neuhauser (2009); Yabana et al. (2012); Lucchini et al. (2016); Li et al. (2016); Yamada and Yabana (2018a, b); Uemoto et al. (2018); Yamada et al. (2018), on the other hand we have also coupling to the full quantized field as discussed in Ruggenthaler et al. (2011); Tokatly (2013); Ruggenthaler et al. (2014a); Ruggenthaler (2015); Galego et al. (2015); Flick et al. (2018a); Feist et al. (2017); Ribeiro et al. (2018). Only very recently also practical formulations that include all three constituents explicitly have emerged de Melo and Marini (2016); Flick et al. (2017a); Flick and Narang (2018); Schäfer et al. (2018). These coupled matter-photon approaches allow to investigate in detail the change of chemical structures due to changes in the electromagnetic vacuum Flick et al. (2018b), coupled light-matter observables such as polariton states and novel potential-energy surfaces Flick et al. (2017b), or changes in Maxwell’s equations due to the interaction with matter Flick et al. (2018c). A further advantage of a coupled light-matter description is that observables that are measured via the light field, such as absorption or emission spectra, do no longer need to be approximately determined by matter degrees but are accessible directly by the calculated photon field Ruggenthaler et al. (2018).

However, the above presented approaches and techniques are themselves either simplifications of the full problem, e.g., by only assuming dipole coupling or neglecting the nuclear degrees of freedom, or have not been made practical for the general case. In this work we will close this gap and provide the first full ab-initio treatment of electrons, nuclei and photons on equal footing via a density-functional reformulation of a generalization of the Pauli-Fierz Hamiltonian of non-relativistic quantum electrodynamics (QED) Spohn (2004); Ruggenthaler et al. (2018) together with a numerical implementation of the resulting Maxwell-Pauli-Kohn-Sham (MPKS) equations. By applying the resulting multi-scale and multi-species formulation to a nanoplasmonic case study, we highlight how discarding electromagnetic and/or nuclear degrees of freedom can alter certain observables as well as how observables differ if they are computed directly from the Maxwell field as opposed to the usual approximate treatment. These results provide a completely new perspective on fundamental low-energy physics, where a disagreement between theory and experiment is often attributed to missing correlations among only one species of particles, i.e., electrons, nuclei or photons, but the discarded degrees of freedom are completely neglected. Furthermore, this unbiased approach allows to tackle many of the above mentioned experimental results and introduces a novel tool to employ the complex interplay between light and matter for the design and control of novel materials.

The paper is structured as follows. In section II, we introduce the fundamentals that lead to the generalized many-body Pauli-Fierz Hamiltonian that we consider as starting point of our approach. Next, we provide in section III a one-to-one correspondence between external and internal variables for the Pauli-Fierz field theory. Based on this we can establish a density-functional theory (DFT) of non-relativistic QED for photons, electrons and effective nuclei. In section IV, we introduce the Kohn-Sham construction for our generalization of quantum-electrodyamical density-functional theory (QEDFT), which leads in the mean-field approximation for the nuclei to coupled Ehrenfest-Maxwell-Pauli-Kohn-Sham (EMPKS) equations. All practical details for the implementation of a solution of these coupled equations are provided in section V. In particular, using the Riemann-Silberstein vector of classical electrodynamics, we rewrite Maxwell’s equations in Schrödinger form and introduce the corresponding time-evolution operators for the homogeneous and for the inhomogeneous cases. For an efficient absorption of outgoing waves, we introduce absorbing boundary conditions and a perfectly matched layer for the Riemann-Silberstein vector. We discuss full-minimal coupling and a multipole expansion and introduce a predictor-corrector scheme for self-consistent forward-backward coupling of light and matter. We place an emphasis on the spatial and temporal multiscale aspects of light-matter coupling and conclude the section with a validation of the method and a comparison to the finite-difference-time-domain (FDTD) approach for solving Maxwell’s equations. Finally, in section VI we illustrate the coupled EMPKS approach for a nanoplasmonic system. We investigate electric field enhancements, harmonic generation and analyze the role of nuclear motion. Since in our approach we also propagate the electromagnetic fields on a grid, we can define electromagnetic detectors in the far field which record all outgoing electromagnetic waves. We conclude the paper in section VII and provide an outlook in section VIII.

II Fundamentals

Let us start by defining some notation. We will use the usual vector notation alongside the relativistic covariant notation. From a fundamental point of view this is convenient since we can easily connect to QED and its relativistic equations, i.e., the Dirac and Maxwell equations. We therefore make a difference between upper and lower indices, which are connected via the Minkowski metric with signature $g\equiv(+,-,-,-)$,

[TABLE]

We denote by $x^{0}=c_{0}t$ the (temporal) zero component of the four-component vector $x^{\mu}$, where Greek indices go over all space-time dimensions, $\mu=\{0,1,2,3\}$, and roman letters in the covariant notation go over only the three spatial dimensions, $k=\{1,2,3\}$. Further we denote the speed of light as $c_{0}$, which is related to the vacuum permeability $\mu_{0}$ and the vacuum permittivity $\epsilon_{0}$ via $c_{0}=1/\sqrt{\mu_{0}\epsilon_{0}}$. To make switching between notations easier, we provide the following table, where the Einstein summation convention over repeated upper and lower indices is implied:

[TABLE]

Here $g_{kl}$ is the three-dimensional (spatial) submatrix of $g_{\mu\nu}$ and $\epsilon^{klm}$ corresponds to the totally anti-symmetric Levi-Civita symbol that arises due the Pauli-matrix algebra (see Eq. (4) below). Accordingly we define $\epsilon_{klm}=\epsilon^{abc}g_{ak}g_{bl}g_{cm}$.

II.1 Relativistic wave equations

Next, we consider the relativistic wave equations that form the basis of coupled light-matter systems. Here we follow the seminal work of Dirac and introduce specific matrix algebras that allow to rewrite a second-order partial-differential equation, which usually corresponds to the relativistic energy-momentum relation $E^{2}=m^{2}c_{0}^{4}+p^{2}c_{0}^{2}$, into a first-order partial-differential equation. The type of matrix algebra to use depends on the statistics of the particle one wants to describe. The most well-known case is of course the Dirac equation, where the spin-$1/2$ nature of the electrons dictates the use of

[TABLE]

Here $\mathbb{1}_{2}$ are two-dimensional identity matrices and the Pauli matrices $\sigma^{k}$ obey

[TABLE]

We see here the Levi-Civita symbol appearing and also that the indices of the matrices are connected via the Minkowski metric, e.g., $\sigma_{k}=g_{kl}\sigma^{l}$. Using these definitions we find that the Dirac operator $\mathrm{i}\hbar c_{0}\gamma^{\mu}\partial_{\mu}$ applied twice to a solution of the Dirac equation $\mathrm{i}\hbar c_{0}\gamma^{\mu}\partial_{\mu}\psi=mc_{0}^{2}\psi$ for a four component wave function $\psi$ implies the relativistic energy-momentum relation, i.e., the Klein-Gordon equation $\hbar^{2}c_{0}^{2}(\nabla^{2}-\tfrac{1}{c_{0}^{2}}\partial_{t}^{2})\psi=m^{2}c_{0}^{4}\psi$.

A similar procedure can also be applied (to some extent at least Gersten (1999)) to particles with other spin. Of particular importance among those are photons, massless spin-$1$ particles. In this case, instead of the spin-$1/2$ Pauli matrices we use the corresponding spin-$1$ matrices in a somewhat non-standard form Gersten (1999)

[TABLE]

While the matrices also obey the spin-algebra $[S^{k},S^{l}]=-\mathrm{i}\epsilon^{klj}S_{j}$ and $\vec{S}^{2}=\mathbb{1}_{3}$, they do not fulfill the same algebra as the Pauli matrices in Eq. (4). Therefore we cannot find a similar simple form of the relativistic energy-momentum relation but instead have side conditions Gersten (1999). This leads to the famous Riemann-Silberstein formulation of electrodynamics Silberstein (1907); Oppenheimer (1931); Keller (2011). To be more specific, we find with the above spin-$1$ matrices Gersten (1999)

[TABLE]

where $E/c_{0}\equiv\mathrm{i}\hbar\partial_{0}$, $p_{k}\equiv-\mathrm{i}\hbar\partial_{k}$ and $\vec{F}$ is the Riemann-Silberstein vector. This vector is a three-component wave function such that we have an entry for each spin state. The above equation holds if we equivalently satisfy

[TABLE]

where we have expressed the momentum mix-term as a side condition on $\vec{F}$. Taking into account that also the complex conjugate of the above equations leads to the right energy-momentum relation and by defining the positive $(+)$ and negative $(-)$ Riemann-Silberstein helicity states

[TABLE]

where $E^{k}$ is the electric field and $B^{k}$ is the magnetic field, the above equations can be recast as the usual homogeneous Maxwell equations in vacuum Gersten (1999)

[TABLE]

We emphasize here that the above side condition translates into the Gauss laws and that the $\hbar$ has canceled out. Only upon quantizing the classical fields (see Sec. II.3) will Planck’s constant reappear in the electromagnetic field equations. In analogy to the Dirac equation we can consider $\vec{F}_{\pm}$ as the single-photon wave function in real space Keller (2011). This analogy will become of practical importance when we actually want to solve the equations of motion for the electromagnetic field numerically (see Sec. V.1). Furthermore, in the case of the photons we refer to the spin as helicity. In the following, in order to make explicit that we only have two physically allowed independent (circular) polarizations, we will always employ the Coulomb gauge when using the vector-potential formulation to couple to matter. In this case the above homogeneous Maxwell equations can be compactly reformulated by introducing vector potentials $A_{\mu}$ that obey the Coulomb gauge condition $\partial_{k}A^{k}=0$ and the second-order relativistic wave equation

[TABLE]

To connect to the previous versions of the homogeneous Maxwell equations we only need the relation between the vector potential and the physical fields, i.e., $E^{k}=-\partial_{0}A^{k}$ and $B^{k}=-\frac{1}{c_{0}}\epsilon^{klm}\partial_{l}A_{m}$.

One can extend the above considerations to arbitrary spins (also for particles with mass), which leads to the Bargmann-Wigner equations. However, finding simple and physical side conditions, as Gauss’ law in the above photon case, is not always possible. This is one reason why in the following we will use the non-relativistic limit of the Dirac equation (also for the electrons) and then replace the spin-$1/2$ matrices by different spin matrices for each species of nuclei. In this way we only keep the most important degrees of freedom of the nuclei in our considerations, i.e., quantized translations, rotations and vibrations, and do not describe the protons and neutrons explicitly, which themselves are effective spin-$1/2$ particles. The second, more relevant reason why we will not use relativistic equations to describe the matter degrees of freedom (while the photons are treated fully relativistically) is that we would like to have a stable ground state and a mathematically well-defined non-perturbative theory Spohn (2004). Without further restrictions the Dirac equation does not have a ground state, as can be seen from the fact that the spectrum is in general unbounded from below Thaller (2013). Taking the non-relativistic limit of the Dirac equation minimally coupled to a classical external vector potential $a_{\mu}$, i.e., $\partial_{\mu}\rightarrow\partial_{\mu}+\mathrm{i}qa_{\mu}/\hbar c_{0}$ with $q$ being the charge of the particle species, leads in first order of $1/Mc_{0}^{2}$ to the Pauli-equation

[TABLE]

Here we have already replaced the Pauli matrices with general spin matrices $S^{k}$, and $M$ is the mass of the particle and we used $\vec{b}(\vec{r},t)=\frac{1}{c_{0}}\vec{\nabla}\times\vec{a}(\vec{r},t)$. This generalized Pauli equation for arbitrary spin will be, together with the homogeneous Maxwell equations for the photons, our mathematical representation of the basic building blocks of our theory for electrons, nuclei and photons. Furthermore, we use the notation convention that an external field, i.e., a field that is not part of the modelled light-matter system such as a classical pump or probe pulse, is denoted with a lower-case letter. We will next combine all these ingredients into one general framework, which will be a generalized form the Pauli-Fierz Hamiltonian Spohn (2004). For completeness, we note that semi-relativistic extensions of the Pauli-Fierz Hamiltonian exist Miyao and Spohn (2009); Könenberg et al. (2011); Hiroshima (2014), but as starting point we stay within the non-relativistic limit for the matter subsystem. This limit is already enough for a vast set of applications.

II.2 Free matter Hamiltonians

Let us start by merging the different mathematical representations of our fundamental building blocks. We first consider the matter subsystems, i.e., electrons and effective nuclei. Since at this point we are not yet coupling the matter subsystems to the quantized electromagnetic field, i.e., the photons, these particles are not interacting. In the end the photons are the gauge bosons that make charged particles interact. We therefore follow the usual construction of quantum field theory to establish interacting theories Greiner and Reinhardt (2013).

In order to consider many non-interacting particles we will lift our single-particle description of Eq. (15) to arbitrarily but finitely many particles. Formally this is done most easily by working in Fock space. We therefore introduce Fock-space creation and annihilation field operators that obey

[TABLE]

where $s$ corresponds to the different possible spins of the particles, and $\pm$ refers to anti-commutation (fermions) or commutations (bosons) relations, respectively. Mathematically these objects are somewhat inconvenient Thirring (2013) but they allow for very efficient formal manipulations. Thus we introduce the “second-quantized” notation for computational convenience. However, since we work with number-conserving Hamiltonians, we can always switch back to the usual “first quantized” form in which every object can be made well-defined. Only for the photons the Fock space is necessary. In this case, however, the mathematical problems can be kept in check Spohn (2004). If we then introduce the conventions

[TABLE]

we can lift the single-particle Pauli equation to the Fock space via

[TABLE]

We can do this now for every species of particles. If we have $N$ different species of particles, i.e., electrons and effective nuclei, we then have a direct sum of Fock spaces. The resulting Hamiltonian on this sum of Fock spaces is given with the definition of the respective field operators, masses, charges and spin matrices

[TABLE]

as

[TABLE]

We note here that while all of the different particles do not interact with each other, we still assume that they all see the same external field $a_{\mu}$ in Coulomb gauge. By lifting this external classical field to a quantum field we will make the particles interact.

II.3 Free photon Hamiltonian

Before we do so, we will first quantize the electromagnetic field. We will follow the standard procedure Greiner and Reinhardt (2013), but will also connect to the Riemann-Silberstein formulation of QED Keller (2011). Since we have chosen the Coulomb gauge, the canonical quantization procedure only affects the transversal fields Greiner and Reinhardt (2013) and the canonical commutation relations read

[TABLE]

where we employed the transversal delta distribution

[TABLE]

Here $1/\Delta$ is the inverse of the Laplacian $\Delta\equiv\vec{\nabla}^{2}$. Due to this quantization procedure in Coulomb gauge, the longitudinal part of the electromagnetic field stays classical and does not influence the quantized degrees of freedom. This can be seen most easily if we construct the Hamiltonian of the free Maxwell field. To do so we introduce the vector-potential operator in terms of creation and annihilation field operators in momentum space

[TABLE]

where $\omega_{k}=c_{0}|\vec{k}|$ and $\vec{\epsilon}\,(\vec{k},\lambda)$ is the transversal polarization vector that obeys $\vec{k}\cdot\vec{\epsilon}\,(\vec{k},\lambda)=\vec{\epsilon}\,(\vec{k},1)\cdot\vec{\epsilon}\,(\vec{k},2)=0$ Greiner and Reinhardt (2013). The momentum-space annihilation and creation field operators obey the usual commutation relations. The purely transversal electric field is then given in accordance to the classical case by $\partial_{0}\hat{A}^{k}=-\hat{E}^{k}_{\perp}$ as

[TABLE]

and the magnetic field via $\hat{B}^{k}=-\tfrac{1}{c_{0}}\epsilon^{klm}\partial_{l}\hat{A}_{m}$ as

[TABLE]

Again, following the classical definition of the energy of the electromagnetic field, we find

[TABLE]

where we used normal ordering, denoted as $::$, to discard the constant energy shift Greiner and Reinhardt (2013) and included the coupling to a classical external charge current $j^{\mu}$. Using the Green’s function of the Laplacian in real-space representation, i.e., $G(\vec{r},\vec{r}^{\prime})=\left\langle\,\vec{r}\,\left|\,(\Delta^{-1})\vec{r}\hskip 1.42262pt^{\prime}\,\right.\right\rangle$, we can further express the zero component of the field as

[TABLE]

This, together with $\vec{E}_{\parallel}(\vec{r},t)=-\vec{\nabla}A_{0}(\vec{r},t)$, allows us, after partial integration, to rewrite Eq. (26) as

[TABLE]

Here we have defined

[TABLE]

However, the last term in Eq. (28), which corresponds to the longitudinal degrees of freedom of the field, is merely a constant (not an operator) and can thus be discarded in the current case. It commutes with all observables. When we couple to the quantum particles, the external current will correspond instead to an operator-valued field and thus will no longer vanish. Indeed, it is this part that will lead to the longitudinal Coulomb interaction among the charged particles (see Sec. II.4.1).

Before we move on and derive the inhomogeneous Maxwell equation from the above Hamiltonian, let us rewrite the previously introduced quantum field in a form such that we can easily connect it to the Riemann-Silberstein formulation. In accordance to the classical case we can introduce the Riemann-Silberstein operators

[TABLE]

Thus, we see that we can use the expectation value of the Riemann-Silberstein vectors $\vec{F}_{\pm}$ to re-express the transversal electric and magnetic fields. Some further analysis Keller (2011) shows that one can furthermore decompose the operators and the expectation values in positive and negative frequency parts, which then give rise to helicity creation and annihilation field operators and helicity single-photon wave functions.

Irrespective thereof, from the Hamiltonian of Eq. (28) describing a photon field coupled to a classical external current, we can derive the operator form of the inhomogeneous Maxwell equation in Coulomb gauge by applying the Heisenberg equation of motion twice, i.e.,

[TABLE]

Here we have assumed in the last step that the external classical charge current $j^{\mu}$ obeys the continuity equation $\partial_{\mu}j^{\mu}=0$, such that by the Helmholtz decomposition only the transversal part of the charge current $j^{i}_{\perp}$ couples to the purely transversal photon field. As pointed out before, the longitudinal component of the photon Hamiltonian does not influence the quantized degrees of freedom since it commutes with the field operators. The classical component is determined by Eq. (27).

II.4 Interaction Hamiltonians

After providing the basic uncoupled Hamiltonians of non-interacting particles of different species in Eq. (21) and of uncoupled photons in Eq. (28), we now join them. Following the minimal-coupling prescription of QED Greiner and Reinhardt (2013), we know that we merely have to use the conserved total charge current of all species and, in accordance to classical electrodynamics, couple it linearly to the vector-potential operator Ruggenthaler et al. (2014b). This can be done by promoting the above external classical fields to operator-valued fields. Since we work in Coulomb gauge we can do the coupling conveniently in two consecutive steps: first only due to the longitudinal and then due to the transversal electromagnetic field.

II.4.1 Longitudinal interactions

The zero component of the charge-current operator for the multi-species case is

[TABLE]

In accordance to the units of a charge current we have multiplied the usual densities $\hat{\Phi}^{\dagger}_{(n)}\hat{\Phi}_{(n)}$ not only by the respective charge of the species $q_{(n)}$ but also by the velocity of light $c_{0}$. The longitudinal part of the photon field now gets an operator-valued contribution and thus can no longer be discarded. The corresponding term then reads

[TABLE]

Here we used normal ordering to bring the interaction in the usual Coulomb form. We note that in the third line, since we consider different species, the normal ordering does not affect the expression as the operators commute. This term we refer to as the inter-species Coulomb interaction, i.e., how the different effective nuclei and electrons act on each other via longitudinal photons. The last line describes how the particles of each individual species interact longitudinally with each other. We call this term the intra-species Coulomb interaction. For later reference and also to highlight the nature of the somewhat unfamiliar term of the inter-species Coulomb interaction, we note that $w(\vec{r},\vec{r}\hskip 1.42262pt^{\prime})$ is the real-space Green’s function of the Poisson equation. This allows us to define the operator-valued scalar field that one specific species feels due to all the other species, i.e.,

[TABLE]

Though now operator-valued, this term still only affects the matter subsystem and commutes with the photon field observables. Since we have expressed the longitudinal interaction purely in matter degrees of freedom we get an interaction among the particles but not with the transversal photon field. The resulting equations would be the interacting Pauli equation for a multi-species and multi-particle problem and uncoupled (transversal) photons. However, to make the resulting Pauli equation physically reasonable we still need to change the masses of the particles. Indeed, to agree with the observed spectrum, we need to use the renormalized physical masses that take into account the effect of the transversal photon degrees instead of the bare (unobservable) masses with which we built the coupled problem Ryder (1996); Spohn (2004). The physical mass is always a sum of the bare plus the electromagnetic mass. The electromagnetic mass just subsumes the energy that is stored in the transversal photon field , which is always non-zero when coupled to charges. In this way even when we discard the coupling to the transversal photon field, it implicitly shows up in our physical masses of electrons and effective nuclei.

II.4.2 Transverse interactions

Let us next consider the coupling between the transversal degrees of freedom of matter and light. To be as general as possible at this point we will not only consider the internal degrees of freedom but also couple to classical external fields. To do so we introduce total fields that have an operator-valued and a classical contribution, i.e.,

[TABLE]

Here, due to the fact that the Coulomb gauge only quantizes the physical (transversal) photon degrees of freedom, an external charge density $j^{0}$ couples via the Coulomb kernel directly to the charged particles. Therefore, physically an external classical density becomes equivalent to an external scalar potential $a^{0}$. Since the external charge density effectively only changes the external classical scalar potential, we need to avoid this sort of “double counting” when we want to establish a Runge-Gross-type mapping. We then couple these total fields to the three spatial components of the total conserved charge current of the matter system

[TABLE]

The first term is the total paramagnetic current density given by

[TABLE]

the second is the magnetization current due to the Stern-Gerlach (Pauli) term

[TABLE]

and the last term is the diamagnetic term. That the photon field becomes part of the charge current is due to the quadratic part in the Pauli minimal coupling, which in turn is due to expressing the anti-particle (positronic) degrees of freedom by the particle degrees of freedom in the non-relativistic limit of the Dirac equation Ruggenthaler et al. (2014b). Since we only couple to the transversal part of the photon field it is also only the transversal part of the total conserved charge current $\hat{J}_{\perp}^{k}$ that is needed. Therefore, the zero component of the charge current, i.e., $\hat{J}^{0}$, does not couple to the transversal photons but only to the previously studied longitudinal part of the electromagnetic field. The resulting fully coupled generalized non-relativistic QED Hamiltonian is given in the next section. For completeness we note that in order to have a well-defined self-adjoint operator, a square-integrable mode mask function should be used Spohn (2004). In its simplest form this is just a cut-off such that the integrals over the photon modes stop at some highest allowed frequency. Physically, since we treat non-relativistic particles, a sensible choice is the rest-mass energy of the electrons at about $\hbar\omega_{k}\approx 0.5110$ MeV. For such high energies the Pauli equation becomes unreliable. Note that this cut-off implies that also the choice of the bare mass depends on this cut-off. For instance, if the cut-off is chosen to be zero, i.e., we only have longitudinal coupling to the Maxwell field, the masses are the physical masses.

III QEDFT for multi-species

Let us now collect all the previous results of the last section and give the generalized Pauli-Fierz Hamiltonian of interacting multi-species matter systems coupled to photons and to external classical electromagnetic fields and charge currents

[TABLE]

with the following definition of the canonical momentum

[TABLE]

Here and in the following we have assumed that also the external electromagnetic field is given in Coulomb gauge

[TABLE]

We also assumed that we have only an external classical scalar potential $a^{0}$ and no external classical charge density $j^{0}$, which due to Eq. (35) would only modify the scalar potential. Thus different configurations of $a^{0}$ and $j^{0}$ can lead to the same physics. Consequently, since we assumed that the external classical current obeys the continuity equation (see Subsec. II.3), we also find

[TABLE]

We note how enforcing that we do not “double count” physically equivalent situations leads to an equivalent “gauge” condition on the classical charge current. This reduction to inequivalent external fields is a prerequisite for a Runge-Gross-type mapping. Furthermore, one can check for the consistency of our construction by calculating the operator-valued continuity equation. By determining the Heisenberg equation of motion of the zero component of the charge current one readily finds

[TABLE]

i.e., the total conserved current is the afore defined operator of Eq. (36). Here it is important to note a common pitfall. Since the magnetization current $\hat{J}_{\rm{mc}}^{k}$ is given in terms of a curl, by construction it does not appear in the continuity equation. Only upon careful consideration of how to express the energy in terms of the linear coupling between current and vector potential or via the Gordon decomposition of the relativistic charge current Engel and Dreizler (2011); Ruggenthaler et al. (2014b) does the full charge current appear. This is the reason why the magnetization current is often overlooked.

In a next step we want to establish a bijective mapping between the external fields $(a_{\mu},j_{k})$, also called the external pair, and the internal pair $(J^{\mu},A^{k})$, which are given by the expectation value of the corresponding operators for the wave function $\Psi(t)$ that is determined by propagating a fixed initial state $\Psi_{0}$ with the Hamiltonian of Eq. (39). Hereby we indicate the dependence of the Hamiltonian on the chosen external pair by $\hat{H}[a_{\mu},j_{k}]$ and, correspondingly, the propagated wave function also becomes dependent on the external pair, i.e., $\Psi([a_{\mu},j_{k}],t)$. Establishing a one-to-one correspondence between $(a_{\mu},j_{k})$ and $(J^{\mu},A^{k})$ would then allow to re-express the dependence of the wave function of the fully coupled system in terms of the internal pair, i.e., $\Psi([J^{\mu},A^{k}],t)$. This makes all observables expressible in terms of the internal pair only and allows to re-express the full Hamiltonian equation as an exact quantum-fluid equation Ruggenthaler et al. (2014b). This in turn leads to the possibility of a Kohn-Sham-type construction for the coupled matter-photon problem at hand (see Subsec. IV). In order to establish such a one-to-one correspondence we will follow the approach of van Leeuwen for the purely electronic case van Leeuwen (1999) and combine it with the derivations of Ref. Ruggenthaler et al. (2014b), where electrons are coupled to photons in a similar setting as the one considered here. As a first step we therefore establish the fundamental equations of motion for the (still operator-valued) internal pair. Since the first time-derivative of $\hat{A}^{k}$ merely leads to the conjugate field $\hat{E}_{\perp}^{k}$, we need to go to the second time-derivative, which leads to the inhomogeneous Maxwell equation in Coulomb gauge

[TABLE]

Here the last term on the right-hand side merely takes care to only count the transversal part of the current operator, i.e., it subtracts the longitudinal component with an operator-equivalent of the Helmholtz decomposition. Alternatively, we could also write $\hat{J}^{k}_{\perp}$, which implies we only take the transversal component of $\hat{J}^{k}$. Again we note that the external current is chosen to be purely transversal, i.e., it obeys Eq. (41). For the current we can use directly the first time derivative. Instead of summing over all contributions directly let us give the equation of motion for the species current, which after lengthy derivations reads as

[TABLE]

In accordance to the usual electronic case, we have defined here the momentum stress tensor

[TABLE]

and the interaction stress tensor

[TABLE]

Furthermore, in correspondence with the current density, we have also defined a magnetization density of the $n$-th particle species

[TABLE]

as well as the corresponding expressions for the $n$-th paramagnetic current. In this somewhat complicated expression, which is consistent with previous electron-photon Ruggenthaler et al. (2014b) and electron-only results Stefanucci and Van Leeuwen (2013), the term of main interest in the following will be

[TABLE]

Here the first term on the right-hand side will be the reason why we are able to establish a one-to-one correspondence, i.e., a Runge-Gross-type result. Due to the fact that we have a $q^{2}_{(n)}$ factor, all different particles species add up positively and the expectation value for a physical wave function will be usually strictly positive in all of space, i.e., non-zero everywhere. Specifically, in the following we will choose an initial state $\Psi_{0}$ that obeys $\sum_{n}\left\langle\,\Psi_{0}\,\left|\,\hat{\Phi}_{(n)}^{\dagger}\hat{\Phi}_{(n)}\Psi_{0}\,\right.\right\rangle>0$ in all of space. While we already need to assume such a strict positivity for electrons only, we now have a positive contribution from each particle species leaving this assumption rather trivial to be fulfilled.

Before proving a bijective mapping (under certain assumptions) we introduce some further convenient notation. We rewrite the equation of motion for the charge current in the following compact form

[TABLE]

where we have subsumed all the complex expressions of how particles and photons couple with each other in

[TABLE]

Here it is important to note that no time-derivatives of the external potentials appear in $\hat{Q}^{k}$. In a next step we define higher-order time derivatives of expectation values of operators at $t=0$ by

[TABLE]

and accordingly for the classical external fields by

[TABLE]

where we put the time-derivative in brackets, i.e., $\{\alpha\}$, to distinguish them from other components. Here $\alpha\in\mathbb{N}_{0}$ and we assume from now on that all external fields and expectation values are real analytic in time, i.e., that their Taylor series in time has a finite convergence radius. In general one can get away with fewer conditions Ruggenthaler et al. (2015) but for the sake of simplicity we stick with these rather stringent ones. Having assumed that all time-derivatives at zero exist (which implies a rather well-behaved initial state) we can, based on the Eq. (49), deduce how all the different time-derivatives at $t=0$ are connected. Following van Leeuwen van Leeuwen (1999) we find

[TABLE]

Accordingly we can define a similar equation based on Eq. (43) for the vector-potential expectation value and the external classical charge current

[TABLE]

We can now show constructively that there is a bijective mapping $(a_{\mu},j_{k})\leftrightarrow(J^{\mu},A^{k})$ by prescribing the Taylor series of an internal pair and then constructing a unique external pair that generates $(J^{\mu},A^{k})$ by propagating the given initial state $\Psi_{0}$ with the corresponding Hamiltonian $\hat{H}[a_{\mu},j_{k}]$.

To do so we first of all need to guarantee that the prescribed internal current and vector potential at $t=0$ agrees with the fixed initial state $\Psi_{0}$. Therefore it needs to hold that

[TABLE]

The last condition guarantees that also $J^{0}$ is exactly reproduced via the continuity equation as we go along in time. We note that due to the fact that the current at the initial time also contains $a^{k,\{0\}}$ we get a first condition on the chosen current as well as the $\alpha=0$ component of the external field by

[TABLE]

where all other quantities are given via the initial state. Due to the fact that we have chosen the Coulomb gauge, the left hand side should obey $\partial_{k}a^{k,\{0\}}=0$. This in turn implies that also the right-hand side should be purely transversal. If this is not the case, and since we know that the charge current is gauge-independent, we would need to perform a gauge transformation (which we now know)

[TABLE]

on the initial state such that it complies with the Coulomb gauge or we use a different gauge choice. For the other $\alpha=0$ components we merely employ the equations of motions from Eqs. (36) and (52), which leads to

[TABLE]

Here $Q^{k,\{0\}}$ contains only the previously determined components of the external vector potential. For the higher orders of the external pair we now only use the equations of motion, i.e.,

[TABLE]

The $a^{k,\{1\}}$-term we get again in two steps. From the first equation above (in the case $\alpha=0$ the sum over $\beta$ is to be discarded) we determine an $\tilde{a}^{k,\{1\}}$, and from the condition that the resulting field is purely transversal we determine the corresponding gauge function in next order, i.e.,

[TABLE]

which then leads to $a^{k,\{1\}}$. We indicate this procedure with a $\perp$ at the right-hand side of the equation. This procedure is possible since only the previously determined orders appear in $Q^{k,\{1\}}$. The first-order for the external current is simply determined from plugging in the prescribed internal pair. We can now repeat this order by order and with this construct the corresponding external pair

[TABLE]

In this procedure we see the necessity of the afore-mentioned (after Eq. (48)) strict positivity of $\rho^{\{0\}}$. If $\rho^{\{0\}}$ would be zero somewhere, the external field would be undetermined by this construction. Let us finally point out that within a given gauge choice only one such pair can exist. For the $j^{k}$ this is trivial due to the linearity of the Maxwell equation and the “gauge” condition $\partial_{k}j^{k}=0$. For the vector potential in zeroth order, i.e., Eq. (III), we immediately see that the difference can only be due to a gauge transformation. Fixing the gauge leaves us with only one choice. This in turn fixes the zero component of the scalar potential, i.e., Eq. (57). In the next order, i.e., Eq. (52) for $\alpha=1$, we are then again left with only a gauge choice and by fixing $\partial_{k}a^{k,\{1\}}=0$ we find the unique representative. In this way we can go through all orders and find one and only one $a^{\mu}$.

IV Maxwell-Pauli-Kohn-Sham coupling

Although we could have shown the one-to-one correspondence more easily by following the original Runge-Gross approach Runge and Gross (1984), the constructive approach by van Leeuwen van Leeuwen (1999) is beneficial for showing the existence of a Kohn-Sham-type system. But before, let us briefly explain how the above one-to-one correspondence becomes relevant for practical applications.

In the previous section we established that in principle the full wave function of the generalized Pauli-Fierz Hamiltonian of Eq. (39) is determined uniquely by the given initial state $\Psi_{0}$ and the internal pair $(J^{\mu},A^{k})$. From this we realize that instead of solving the wave-function-based Hamiltonian formulation of the problem (that due to the infinitely many degrees of freedom of the photon field is impossible even for only one particle) we can alternatively solve two coupled non-linear fluid equations

[TABLE]

Here $\rho([J^{\nu},A^{l}];\vec{r},t)$ is a functional of the internal pair and initial state respectively (we do not indicate the dependence on the initial state for simplicity) and $Q^{k}([J^{\nu},A^{l};a^{\tau}];\vec{r},t)$ does also depend on the chosen external vector potential (see the partial definition in Eq. (III) and Eq.(44), respectively). It is easy to rewrite the term $Q^{k}$ in such a way as to make the dependence on $a^{\mu}$ explicit, but for the sake of brevity we refrain from this here. The above coupled equation are the exact quantum Navier-Stokes equations Tokatly (2005); Ruggenthaler et al. (2014b, 2015) and make calculations of particles coupled to the photon field practical. In principle, for a given external pair the solution of these two coupled equations would lead to the exact internal pair, from which we could construct all physical observables. The only severe drawback is that we do not know the expressions of $\rho$ and $Q^{k}$ in terms of the currents and fields, but only in terms of the intractable wave function, so we have to apply some approximations in practice. Instead of trying to find approximations in terms of $J^{\mu}$ and $A^{k}$ directly, we will follow the well-established approach to employ a numerically simpler quantum system that shares as much similarities with the full generalized Pauli-Fierz problem and then build approximations on top of this ersatz system. To put it differently, we will use a still numerically solvable system to build approximations to $\rho$ and $Q^{k}$ in terms of the resulting auxiliary wave functions. This procedure is called a Kohn-Sham construction Kohn and Sham (1965); Tokatly (2005); Ruggenthaler et al. (2014b, 2015).

While we have many possibilities, the simplest choice for an auxiliary Kohn-Sham-type system is a non-interacting (and consequently also uncoupled) system as already introduced for the electron-photon case Ruggenthaler et al. (2011, 2014b). Similarly to the fully interacting case, we can find a one-to-one correspondence between external and internal pairs

[TABLE]

Here we follow the usual convention to denote the external fields for a non-interacting system with a subindex “$\rm{s}$”. By following the construction of Sec. III, we can also provide a one-to-one correspondence for the non-interacting case and thus generate two maps

[TABLE]

of a given internal to different external pairs for given initial states $\Psi_{0}$ (of the interacting) and $\Phi_{0}$ (of the non-interacting), respectively. We therefore see that the above extended Runge-Gross-type approach helps us in showing that a fully interacting problem is representable by a non-interacting problem. Since the Kohn-Sham photon field is uncoupled, instead of solving for the Kohn-Sham photon wave function for infinitely many degrees, we can just solve the corresponding classical inhomogeneous Maxwell equation Ruggenthaler et al. (2014b). Both the classical as well as the fully quantized yet uncoupled Kohn-Sham photon equation lead, by construction, to the same field $A^{k}$. This allows us to use a matter wave function $\Phi_{0}$ of non-interacting particles only, which we assume to obey

[TABLE]

The initial state of the fully coupled wave function $\Psi_{0}$ then further provides the initial condition for the classical photon field, i.e., Eqs. (III). For the external current, due to the linearity of the Maxwell equation, the mappings are the same and hence lead to

[TABLE]

Following the usual Kohn-Sham construction Ruggenthaler et al. (2015) we merely need to introduce a mean-field exchange-correlation (Mxc) potential

[TABLE]

to take care of the difference between the physical and the auxiliary system. Alternatively we can directly construct the Kohn-Sham vector potential $a_{\rm{KS}}^{\mu}=a^{\mu}+a^{\mu}_{\rm{Mxc}}$ by employing the Eqs. (58) for the non-interacting case together with Eq. (52) for the interacting case to express the unknown $J^{k,\{\alpha+1\}}$ in terms of the basic variables. With these definitions in place, assuming that the initial state is a tensor product of Slater determinants (fermions) and permanents (bosons), we find the auxiliary Maxwell-Pauli-Kohn-Sham (MPKS) equations

[TABLE]

[TABLE]

[TABLE]

where the internal charge current is determined by the auxiliary wave function $\Phi(t)$. This wave function consists of auxiliary single-particle spin-orbitals $\phi_{(n,i)}(t)$, where we have a corresponding spin-degree of freedom $s_{(n)}$ for each species $n$. To distinguish the indices $n$ and $i$ from covariant indices we put them in parenthesis. We stress that, provided we have the exact Mxc potential, the coupled auxiliary MPKS problem predicts the exact internal pair $(J^{\mu},A^{k})$ for a given generalized Pauli-Fierz Hamiltonian $\hat{H}[a^{\mu},j^{k}]$. The infinite number of degrees of freedom of the quantized photon field has been exactly subsumed in the classical Maxwell equation and in the non-linear coupling to the matter subsystem. In this way the MPKS equation takes into account the full photon and phonon bath in an exact manner.

IV.1 Classical limit for Nuclei

At this point we note that we have quite some freedom in establishing the mappings as well as the MPKS systems. For instance, we can use instead of the uncoupled particles with their respective bare masses just the uncoupled particles with their physical masses. That is, we already subsume as usually done in quantum mechanics, the fluctuations of the bare electromagnetic vacuum in the physical (dressed) masses Ruggenthaler et al. (2018). This we will do in the following, i.e., the MPKS systems is build in practice by using the physical masses of the particles. Furthermore, we could look at each individual particle-species’ internal current and find that also each species’ internal current can be used to establish a mapping individually. We could then (unphysically) assume that each species sees a different external field and then establish purpose built $a^{\mu}_{\rm{Mxc},(n)}$. Such an approach might be easier in establishing more accurate approximations to the ultimately unknown Mxc potentials. Here we will not follow this route but instead try to find a first approximation that simplifies the MPKS construction slightly. Even though we have rewritten the coupled generalized Pauli-Fierz problem in terms of single-particle quantum equations, the solution of a large amount of these non-linear equations is still numerically demanding. Furthermore, for the initial states, which can be determined from a ground-state reformulation of the generalized Pauli-Fierz problem following Ref. Ruggenthaler (2015), it is often beneficial to make the Born-Oppenheimer approximation and treat the nuclei semi-classically. We point out, that due to the coupling of the photons to the electron-nucleus system the problem without external fields is now invariant with respect to the total momentum and total angular momentum of the coupled matter-photon system Spohn (2004), not the coupled matter system only. Thus when translating or rotating only the matter-system in real space, the photon field is changed, which breaks the usual real-space translational and rotational symmetry. This can be beneficial to overcome the usual drawback of electron-nucleus quantum mechanics (and the corresponding DFT formulation) that the densities of such matter systems are homogeneous and no direct molecular structure is apparent. The physical rationale to then perform the classical limit for the nuclei subsystem is that the nuclei are much heavier than the electrons and hence they behave more classically. In the following we will therefore simplify the MPKS construction slightly and describe the nuclei classically. We note that more advanced alternatives exists that, e.g., are based on the exact factorization of electron-nuclei wave functions Abedi et al. (2010, 2012).

In order to see whether this is possible in a simple manner let us consider the MPKS equations for the nuclear orbitals. A first approximation will be to discard the Stern-Gerlach term. Next we rewrite the spatial orbitals in polar representation

[TABLE]

If we then plug this into the MPKS equation without the Stern-Gerlach term (which makes the solution spin-independent) we find a Hamilton-Jacobi-type equation for the phase Min et al. (2017)

[TABLE]

where we denote $a^{\mu}_{\rm{KS}}=a^{\mu}+a^{\mu}_{\rm{Mxc}}$. Using the above equation for the phase we will turn the quantum evolution equation in a classical equation in the following. The physical argument is based on the fact that the $\hbar^{2}/2M_{(n)}$ factor becomes very small for heavy particles such as nuclei, and can thus be discarded in a first approximation. But before we do so, let us connect the relevant quantum observables to their classical counter parts. To do so we first consider the conserved current for Eq. (68) without the Stern-Gerlach term, which is $\vec{J}_{(n,i)}(\vec{r},t)=\vec{J}_{\rm{pmc}}^{(n,i)}(\vec{r},t)+\vec{J}_{\rm{dmc}}^{(n,i)}(\vec{r},t)=\tfrac{q_{(n)}}{M_{(n)}}|\phi_{(n,i)}(\vec{r},t)|\vec{\nabla}S_{(n,i)}(\vec{r},t)-\tfrac{q_{(n)}^{2}}{M_{(n)c_{0}}}\vec{a}_{\mathrm{KS}}(\vec{r},t)|\phi_{(n,i)}(\vec{r},t)|$. With this the total velocity field becomes

[TABLE]

And accordingly we can define the total momentum field $\vec{p}_{(n,i)}=\vec{v}_{(n,i)}/M_{(n)}$.

If we discard the quantum-potential part in Eq. (71) assuming the classical limit $\hbar\rightarrow 0$ for the slowly moving nuclei, then the resulting equation becomes

[TABLE]

where $\vec{b}_{\mathrm{KS}}(\vec{r},t)=\tfrac{1}{c_{0}}\vec{\nabla}\times\vec{a}_{\mathrm{KS}}(\vec{r},t)$. If we then add to both sides the partial time-derivative of the Kohn-Sham vector potential and define the total derivative for a co-moving reference frame that moves with the velocity field $\vec{v}_{(n,i)}$, i.e.,

[TABLE]

then the above equation becomes

[TABLE]

Here we used that $-\partial_{0}\vec{a}_{\mathrm{KS}}=\vec{E}_{\perp}^{\mathrm{KS}}$ gives the transversal electric field and $-\vec{\nabla}a^{0}_{\mathrm{KS}}=\vec{E}_{\parallel}^{\mathrm{KS}}$ gives the longitudinal electric field such that $\vec{E}_{\mathrm{KS}}=\vec{E}_{\perp}^{\mathrm{KS}}+\vec{E}_{\parallel}^{\mathrm{KS}}$. This is just the classical Lorentz equation which can be solved by the method of characteristics, i.e., we can follow a specific classical trajectory that starts at $\vec{r}_{(n,i)}$ and $\vec{p}_{(n,i)}$. The initial wave function then gives us the initial distribution of these trajectories. Using this classical approximation we can determine the charge current of the nuclei that contribute to the total current and thus the effective Kohn-Sham field $a^{\mu}_{\rm{KS}}$. In the case that we use classical trajectories for the nuclei in our MPKS approach, we call the resulting simplification in analogy to matter-only quantum dynamics the Ehrenfest-Maxwell-Pauli-Kohn-Sham (EMPKS) approach.

IV.2 MPKS-Outlook

As is obvious from the extent of the physics included in our considerations – from matter-only quantum mechanics to quantum optics and beyond – the unknown exact effective fields $a^{\mu}_{\rm{Mxc}}$ need to contain all knowledge of the correlated quantum dynamics. Already for electron-only quantum mechanics the quest for such exact expressions is a herculean effort. And this will not become simpler for the full matter-photon problem of electrons, nuclei and photons. Specifically challenging is that the description of the matter subsystem is based on the current density and not on the density as is the usual case in density-functional theory Gross and Dreizler (2013); Ullrich (2011); Engel and Dreizler (2011). Hence most approximations developed for the effective fields only cover the zero component and little is known besides linear-response kernels Vignale and Kohn (1996); Ullrich and Vignale (2002) for the more advanced current-density-functional theory Vignale and Rasolt (1988); Vignale (2004). Therefore a necessary step to develop more accurate functionals beyond treating the zero component via density-functional approximations and the spatial components on the mean-field level includes also more advanced current-density functionals for electron-only theory. Of course we also need to account for new terms that are due to the transversal photon interaction and that are not covered by matter-only theories. It is helpful in this regard that for processes in free space where many photons are involved, e.g., when the system is perturbed by a weak laser pulse, the mean-field limit is very accurate Ruggenthaler et al. (2018); Flick et al. (2015). This changes if we consider situations where it is the fluctuations of the photon field that are important, such as changes in the ground-state of the combined light-matter system Flick et al. (2018b), and few-photon strong coupling situations Flick et al. (2018c). But for such specific cases a slight simplification of the presented theory to dipole coupling Tokatly (2013); Ruggenthaler et al. (2014b); Ruggenthaler (2015) becomes a sensible alternative. In this case it is only the density of the matter-subsystem that couples to the photon field Schäfer et al. (2018) and standard density-functional approximation schemes become applicable. The currently most sophisticated is an extension of the optimized-effective potential approach to matter-photon systems Pellegrini et al. (2015), which has been already applied to real systems Flick et al. (2018b). In this work an exact-exchange-type of approximation for the transversal and longitudinal matter-photon interactions is employed. This means for the novel transversal part that, e.g., leads to the Lamb shift of the matter states, a single-photon approximation is used. Thus multi-photon processes of the fluctuating photon field are not included. An interesting alternative to include all orders of field-fluctuations are trajectory-based methods. In this case instead of one Maxwell field many realization are propagated at the same time, which is akin to a field-quantization procedure in the Riemann-Silberstein formulation of electrodynamics Bialynicki-Birula (1996); Bialynicki-Birula and Bialynicka-Birula (2013). Furthermore, the present formulation allows to inculude many different species of particles. This suggests to model quantum dynamics in a complex environment by, e.g., treating certain particle types with a crude yet efficient orbital-free quantum Navier-Stokes description based directly on Eqs. (61) and (62), while the system of interest is treated with a more accurate Kohn-Sham description.

V Real-space implementation

After presenting the theoretical fundamentals, we introduce in the following sections our first implementation of the coupled EMPKS equations in real-time. Our implementation is based on finite-difference discretizations and real-space grid representations for both the matter wave functions and the electromagnetic fields. While not the only possible choice, this representation has the advantage to allow for a uniform and unbiased description of the combined light-matter system. More importantly, this choice also simplifies the description of coupling between matter and radiation, since also the QED couplings are prescribed in real space. Moreover, the real-space representation is suited for the multi-scale aspects of the coupling, and allows to reuse simulation techniques and algorithms for the matter and the radiation subsystems. We have integrated our EMPKS implementation in the real-space real-time code Octopus Andrade et al. (2015); Castro et al. (2006), an open source simulation package for quantum-mechanical ab-initio calculations based on time-dependent density-functional theory.

We start by going into more detail of the Maxwell’s equations in Riemann-Silberstein representation in section V.1. In section V.2, we introduce a time-evolution operator for the Riemann-Silberstein vector, while section V.3 is focussing on different boundary conditions. We continue our discussion in section V.4 with details on the time evolution of the matter degrees of freedom before we turn our attention in sections V.5, V.6, and V.7 on the coupling of radiation and matter which involves a multipole expansion for the coupling Hamiltonian, a multi-scale implementation, and a predictor-corrector scheme for a self-consistent coupling of radiation and matter. We conclude the discussion of our implementation in section V.8 with a validation of the method and a comparison to FDTD methods.

V.1 Maxwell’s equations in Riemann-Silberstein representation

Most common electromagnetic simulations are based on the FDTD method, which was already introduced in 1966 by Yee Yee (1966); Taflove and Hagness (2005). FDTD uses two grids, shifted by half of the chosen grid spacing. One grid represents the electric field and the other one the magnetic field. Based on this description, the time evolution in FDTD is given by two update equations one for each electromagnetic field. In contrast, in the Octopus code several methods are implemented to approximate the quantum-mechanical time-evolution operator Castro et al. (2004) for the propagation of wave functions in real-time. By transforming Maxwell’s equations into a six-component Riemann-Silberstein representation Mohr (2010), the underlying equation of motion for the Maxwell fields can be expressed as a Schrödinger-like equation. This fact gives us the opportunity to propagate also the electromagnetic field with quantum-mechanical time-evolution methods modified for Maxwell fields. Besides the fact that we can immediately reuse with Octopus an existing time-evolution engine that is very efficient and highly parallelized, the reformulation of Maxwell’s equations in terms of the six-dimensional Riemann-Silberstein vector also has the advantage over FDTD that higher order discretizations for the spatial derivatives can be used which in turn allow for much larger grid spacings and time steps. Furthermore, the very same grid can be used for the electric and for the magnetic field which simplifies the coupling to microscopic matter charge currents. Also, from our experience it allows for more stable time-stepping of coupled radiation and matter, and improves the maintainability of the implementation.

In Eq. (10) we have introduced already the complex Riemann-Silberstein vector which can also be written in the form

[TABLE]

As already stated before, the sign of the imaginary part of the Riemann-Silberstein vector corresponds to different helicities. To convert Maxwell’s equations to Riemann-Silberstein form, we have to be able to keep track of superposition states of different helicities which correspond to a linear combination of $\vec{F}_{+}$ and $\vec{F}_{-}$. Hence, it is useful to combine both vectors in a six-dimensional vector

[TABLE]

The inverse transformation for given Riemann Silberstein vectors $\vec{F}_{+}(\vec{r},t)$ and $\vec{F}_{-}(\vec{r},t)$ then takes simply the form

[TABLE]

[TABLE]

which allows always to reconstruct the electric and magnetic fields from the Riemann-Silberstein vector.

V.1.1 Maxwell’s equations in Schrödinger form

In this section, we generalize the Riemann-Silberstein description of the homogeneous Maxwell’s equations Eqs. (11-14) to the inhomogeneous case. While here we consider classical external charge and current densities, later, for a fully microscopic description, we use the expectation values of the previously derived quantum-mechanical charge density Eq. (32) and current density Eq. (36) which will lead to a modification of the usual Maxwell’s equations.

We start by considering the two Gauss’ laws for the electric and the magnetic field

[TABLE]

where we introduced as the external charge density $\rho(\vec{r},t)=\tfrac{1}{c_{0}}j^{0}(\vec{r},t)$. Later we will then use also the internal charge density $n(\vec{r},t)=\tfrac{1}{c_{0}}J^{0}(\vec{r},t)$. Both Gauss’ laws can be combined into one equation

[TABLE]

The remaining two Maxwell’s equations, the Ampère’s and Faraday’s law

[TABLE]

can also be combined into one evolution equation for the Riemann-Silberstein vectors as

[TABLE]

The first term on the right-hand side of Eq. (85) describes the curl operation and can be represented by the already previously introduced spin-1 matrices for the photons in Eq. (5). Taking them into account, Ampère’s and Faraday’s laws in Riemann-Silberstein form can be written as

[TABLE]

This equation has the form of an inhomogeneous Schrödinger equation with single particle ”photon” Hamiltonian

[TABLE]

As consequence, Maxwell’s equations in Riemann-Silberstein form can be interpreted as the first quantized wave equation for a single photon in real-space. Nevertheless, we are dealing still with a classical equation, since $\hbar$ can be cancelled in all terms in Eq. (86).

V.1.2 Current densities and integral kernels

As already noted, to reach a microscopic description of the electromagnetic fields, the current and charge densities that appear in our Maxwell’s equations in Riemann-Silberstein form have to correspond to the multi-species currents that we introduced in Eq. (36). The current densities consist of three terms, the paramagnetic current term $J_{\mathrm{pmc}}^{k}$, the diamagnetic current term $J_{\mathrm{dmc}}^{k}$, and the magnetization current term $J_{\mathrm{mc}}^{k}$. These contributions to the total current can be expressed in terms of $I_{(n)}$ auxiliary one-body Kohn-Sham orbitals $\phi_{(n,i)}(\vec{r}s_{(n)},t)$, and the Kohn-Sham densities $n_{(n,i)}(\vec{r}s_{(n)},t)=|\phi_{(n,i)}(\vec{r}s_{(n)},t)|^{2}$ for the corresponding species $n$, which gives rise to the charge density $n(\vec{r},t)=\tfrac{1}{c_{0}}J^{0}(\vec{r},t)=\sum_{n=1}^{N}q_{(n)}\sum_{i,s_{(n)}}^{I_{(n)}}n_{(n,i)}(\vec{r}s_{(n)},t)$ in the form

[TABLE]

[TABLE]

[TABLE]

where the summations go over all spin-states $s_{(n)}$, Kohn-Sham orbitals $i$ and species types $n$. All three current terms as well as the total charge density depend on the particle charge $q_{(n)}$, the particle mass $M_{(n)}$, and the single-particle wave functions. While the paramagnetic current and magnetization current do not depend explicitly on the Maxwell fields, the diamagnetic current depends on the vector potential $a_{\mathrm{KS}}^{k}=a^{k}+a^{k}_{\rm{xc}}+A^{k}$, which is implicitly determined by the Riemann-Silberstein vectors $F_{\pm}^{k}$. The mean-field vector potential $A^{k}$ can be expressed in terms of the magnetic field via a Poisson equation for each magnetic field vector component as

[TABLE]

But we can also use the transversal part of the electric field such that the vector potential $A^{k}$ results from an integral over time starting at the initial time $t=t_{0}$ and ending at the current time $t$

[TABLE]

Here $\perp$ again indicates that only the transversal degrees are to be considered. This can be made explicit with the help of the Helmholtz decomposition, as also shown in Eq. (43). Also, the initial value $A^{k}(\vec{r},t_{0})$ is in principle determined by the initial wave function of the interacting system. Alternatively, provided the external field obeys itself the homogeneous Maxwell’s equation $(\partial_{0}^{2}+\partial_{l}\partial^{l})a^{k}=0$, we can also use instead the total field $A^{k}_{\rm{tot}}=a^{k}+A^{k}$ in the above considerations. This is the standard case and in the implementation we usually work directly with the total field and the corresponding Riemann-Silberstein vector. Since the diamagnetic current term is obtained by one of the two different equivalent vector potential expressions from above, the charge current $J_{\mathrm{dmc}}^{k}$ carries non-local contributions, which originate from the electromagnetic fields. Here we make the distinction between an internal part that depends directly on $A^{k}$ and an external part of the diamagnetic current that comes from the external field $a^{k}$ and from the missing exchange-correlation forces $a^{k}_{\rm{xc}}$. Furthermore, as the internal diamagnetic current depends explicitly on the Riemann-Silberstein vectors, we can combine it with the curl operation of Eq. (86) to define a single integral operator. The Hamiltonian-like integral operator $\bm{\hat{\mathcal{H}}}$ then acts on the Riemann-Silberstein vector, which in the following we denote by

[TABLE]

The operator $\bm{\hat{\mathcal{H}}}$ is determined by its corresponding integral kernel $\hat{\bar{\mathcal{H}}}$ and here given by

[TABLE]

The matrix $\hat{\bar{\mathcal{H}}}_{(0)}$ represents the curl operation with the spin-1 matrices in equation (86) without any diamagnetic current

[TABLE]

Here and in the following, we describe the 6x6 matrices which act on the six-dimensional Riemann-Silberstein state $\mathcal{F}$ as a Kronecker product of a 2x2 and a 3x3 matrix. The first 2x2 matrix which we call ”coupling” matrix shows whether the two Riemann-Silberstein vectors $\vec{F}_{\pm}$ couple to each other. If the off-diagonal of this 2x2 matrix is zero, the resulting 6x6 matrix does not couple the two vectors. Furthermore, the second 3x3 kernel matrix, contains all necessary physical variables and operations to satisfy the underlying Maxwell’s equations. Due to the delta functions in $\hat{\bar{\mathcal{H}}}_{(0)}$, the application of $\bm{\hat{\mathcal{H}}}_{(0)}$ then results in a local linear operator acting on the Riemann-Silberstein vector, i.e.,

[TABLE]

where

[TABLE]

In contrast to this, the integral kernel originating from the diamagnetic current contribution in equation (94) carries a non-locality in space due to

[TABLE]

The integrals over the kernel have to be kept explicitly. Note, that the 2x2 coupling matrix contains off-diagonal entries in this case. The prefactor $\kappa$ is given by

[TABLE]

This non-local term that arises due to the microscopic treatment of the matter system therefore leads already for a classical treatment of the photon field to effective ”photon-photon” interactions and changes the usual Maxwell’s equation. This is discussed in more detail in Flick et al. (2018c). The remaining current densities in Eqs. (88) and (90) as well as the external diamagnetic currents can be added to the external current such that we have a inhomogeneous current J_{\mathrm{inh}}^{k}=J^{k}_{\mathrm{pmc}}+J^{k}_{\mathrm{mc}}-\sum_{n}\tfrac{q^{2}_{(n)}}{M_{(n)c_{0}}}\bigl{[}\sum_{i,s_{(n)}}n_{(n,i)}\bigr{]}(a^{k}+a^{k}_{\rm{xc}})+j^{k}, which constitutes a six-dimensional vector for the Riemann-Silberstein Maxwell’s equation given by

[TABLE]

Combining the previous considerations, allows us finally to formulate Ampère’s and Faraday’s law in Riemann-Silberstein representation

[TABLE]

The Maxwell fields that follow this equation are due to microscopic quantum-mechanical currents from Eq. (100). We call these fields in the following microscopic Maxwell fields. We can, however, also include other types of contributions in our implementation. In the following Sec.V.1.3, we will discuss how we can in a similar manner include also macroscopic Maxwell fields due to, e.g., lenses or surfaces. Further note that in contrast to the quantum mechanical second-order derivative operator for the kinetic energy, the integral kernel $\hat{\mathcal{H}}$ contains only first-order derivatives. Without inhomogeneity this reflects the linear dispersion relation for photons. On the other hand, including the inhomogeneity introduces a non-linear dispersion for the photons.

In the discussion so far, we have seen from the matrix operator $\hat{\mathcal{H}}$ in Eq. (94) and its underlying matrix operators $\hat{\mathcal{H}}_{(0)}$, and $\hat{\mathcal{K}}$ in Eqs. (97) and (98) that the two different Riemann-Silberstein vectors $F_{+}^{k}$ and $F_{-}^{k}$ only couple in the presence of a diamagnetic current $J_{\mathrm{dmc}}^{k}$. Without this contribution, the combined Ampère’s and Faraday’s law in Riemann-Silberstein form in (101) decouple into two three-component equations

[TABLE]

for positive helicity fields $F_{+}^{k}$, and

[TABLE]

for negative helicity fields $F_{-}^{k}$.

V.1.3 Helicity coupling in a linear medium

So far, we have only focused on the microscopic Maxwell’s equations and have obtained an explicit coupling between $F_{+}^{k}$ and $F_{-}^{k}$ caused by the diamagnetic current. The formulation in terms of integral kernels, however, allows also to seamlessly combine a microscopic propagation with macroscopic electrodynamics for linear media. This has the merit that macroscopic optical elements like lenses or mirrors can easily be incorporated in the Riemann-Silberstein propagation. In the following, we briefly discuss this relation.

The coupling matrix $\mathcal{K}$ bears a direct similarity to the Riemann-Silberstein Maxwell’s equation in a linear medium Khan (2005); Bialynicki-Birula (1996), where the electromagnetic fields are described by the electric displacement field $\vec{D}(\vec{r},t)$ and the magnetic $\vec{H}(\vec{r},t)$ field. In linear optics, it is assumed that $\vec{D}$ and $\vec{B}$ depend approximately linearly on $\vec{E}$ and $\vec{H}$

[TABLE]

Here, $\epsilon(\vec{r},t)$ and $\mu(\vec{r},t)$ denote the electric permittivity and magnetic permeability in the medium which in general depend on space and time. Referring to the Riemann-Silberstein definition in vacuum in Eq. (30) with constant vacuum electric permittivity and magnetic permeability, the Riemann-Silberstein vector in a linear medium $\vec{F}_{\pm,\mathrm{lm}}$ has to be adapted according to

[TABLE]

Consequently, the speed of light in the medium also depends now on space and time with

[TABLE]

Based on this new definition for the Riemann-Silberstein vectors in Eq. (105), Ampère’s and Faraday’s law for a linear medium Jackson (1998) are then given by Khan (2005); Bialynicki-Birula (1996)

[TABLE]

where the Hamiltonian $\hat{\mathcal{H}}_{\mathrm{lm}}$ is given by

[TABLE]

The first term describes the free evolution in the linear medium

[TABLE]

Note that, due to $c(\vec{r},t)$ this term contains now an explicit spatial and temporal dependence. The remaining term $\hat{\mathcal{K}}_{\mathrm{lm}}$ describes the properties of the linear medium and takes the form Bialynicki-Birula (1996); Khan (2005)

[TABLE]

where $\dot{\epsilon}(\vec{r},t)$ and $\dot{\mu}(\vec{r},t)$ represent the total time derivative of $\epsilon(\vec{r},t)$ and $\mu(\vec{r},t)$. From the 2x2 coupling matrices contained in $\hat{\mathcal{K}}_{\mathrm{lm}}$ it becomes evident that the medium is capable to couple $F_{+}^{k}$ with $F_{-}^{k}$, as expected from linear optics.

Based on the present formulation, it is instructive to further analyze the relation of the coupling that is induced by the diamagnetic contribution of the current density and the coupling that appears in linear and non-linear media. This provides a route to directly connect microscopic electrodynamics with a macroscopic formulation and will be considered in an upcoming publication.

V.2 Time-evolution operator for the Riemann-Silberstein vector

In the previous section, we have illustrated how to express Ampère’s law and Faraday’s law as an inhomogeneous Schrödinger-like equation in Riemann-Silberstein form. In the following, we construct the time-evolution operators for this equation in the homogeneous and inhomogeneous case. We start with the homogeneous case, i.e.,

[TABLE]

The goal is to find an expression for a time-evolution operator that is acting on the Riemann-Silberstein vector $\mathcal{F}(\vec{r},t)$ in analogy to the time-evolution operator used in quantum mechanics for matter wave functions. Since the right hand-side of Eq. (111) contains in general an integral operator with integral kernel $\hat{\mathcal{H}}$, we use as ansatz for the time-evolution operator $\bm{\hat{\mathcal{U}}}(t,t_{0})$ also an integral operator. This operator should obey the evolution equation

[TABLE]

with the following boundary condition and group composition laws

[TABLE]

From Eqs. (112) we then find the implicit definition

[TABLE]

which fulfills the conditions in Eq. (113). Iterating this equation allows to construct a series expansion for the time evolution-operator

[TABLE]

where we used (93), $\vec{r}_{0}=\vec{r}$ and $\mathcal{\hat{T}}$ is the time-ordering operator such that earlier times go to the right. In the homogeneous microscopic case, when the diamagnetic current and $\mathcal{J}$ are zero, the time evolution equation in Eq. (115) simplifies further to

[TABLE]

In this limit the Hamiltonian $\hat{H}_{\mathrm{Mx},(0)}$ is time-independent which allows to write the time-evolution operator for the Riemann-Silberstein vector formally as exponential

[TABLE]

which has the familiar form of evolution operators for matter wave functions. In the inhomogeneous case, but with $j_{\mathrm{mdc}}$ equal to zero, we start with the inhomogeneous Schrödinger-like equation

[TABLE]

Such inhomogeneous time evolution equations are also found in quantum mechanics Serban et al. (2005); Ndong et al. (2009), which can be readily transferred to our Riemann-Silberstein case. By multiplying the time-evolution operator $\hat{\mathcal{U}}_{(0)}(t,t_{0})$ with an additional time-dependent operator $\hat{\mathcal{Z}}(t)$, we propose

[TABLE]

as an ansatz for the solution of the inhomogeneous equation. The time derivative of this expression leads to the time-evolution equation in terms of $\hat{\mathcal{U}}_{(0)}(t,t_{0})$ and $\bm{\hat{\mathcal{Z}}}(t)$. Taking Eq. (112) into account leads to

[TABLE]

This equation has to be equivalent to Eq. (118) such that the last term on the right-hand side of Eq. (120) is equal to $-\mathcal{J}$. This condition gives us that the $\hat{\mathcal{Z}}(t)\mathcal{F}(\vec{r},t_{0})$ operation has to be equal to

[TABLE]

Substituting this result into the ansatz in Eq. (119) yields the final time-evolution equation for the inhomogeneous Maxwell’s equation in Riemann-Silberstein form

[TABLE]

For a non-zero diamagnetic current $j_{\mathrm{mdc}}$, we have to replace the first term on the right-hand side of Eq. (123) by the corresponding Eq. (115). The final equation provides the correct propagation including the diamagnetic current and reads

[TABLE]

We note, that due to the continuity equation, the time-evolution Eq. (123) conserves the Gauß side condition of the electromagnetic field during the propagation.

In general, the time-evolution of the Riemann-Silberstein vector can not be expressed analytically in closed form. However, the form of Eq. (123) is already suitable for a numerical time-stepping scheme. Iterating the time-evolution $(m+1)$ times for a time-step $\Delta t$ yields

[TABLE]

For convergence, the length of $\Delta t$ is chosen such that the propagation stays stable and results in an accurate evolution of the Riemann-Silberstein vector with negligible error. With this, we have now reached a numerical time-stepping scheme for the Maxwell’s equations coupled to microscopic as well as macroscopic matter.

V.2.1 Forward-backward coupling and time-reversal symmetry

The time-evolution of the Pauli-Fierz Hamiltonian in Eq. (39) is given by the time-dependent Schrödinger equation of the fully coupled matter-photon system. If there are no time-dependent external fields acting on the coupled light-matter system it is by construction time reversal symmetric when $t\rightarrow-t$. This symmetry is only strictly given, if we consider the whole coupled light-matter system where both subsystems influence each other. Breaking the forward and backward coupling between matter and electromagnetic fields also destroys the time-reversal symmetry.

As consequence of the time-reversal symmetry for the full photon-matter coupling, the corresponding Maxwell time stepping has to obey the property that a reverse time step from $\mathcal{F}(\vec{r},m\Delta t)$ leads again to the previous result $\mathcal{F}(\vec{r},(m-1)\Delta t)$. We can construct a numerical time-evolution equation based on the enforced time-reversal symmetry (ETRS) Castro et al. (2004) propagator to also numerically ensure time-reversal symmetry for the fully coupled case. The underlying condition of the ETRS algorithm requires that a propagation forward by one step $\Delta t$ starting from $\mathcal{F}(\vec{r},m\Delta t)$ and then half a step backwards with $\Delta t/2$ is equivalent to propagating half a step forward. As result, a modified numerical time-evolution equation for ETRS based on Eq. (124) takes the form

[TABLE]

The integrals in Eq. (124) and (125) can be approximated by the trapezoidal rule so that the numerical time-evolution equations take the explicit form

[TABLE]

for the simple direct propagation, and

[TABLE]

for the ETRS propagation.

Note, that a strict time-reversal symmetry can only be obtained if the propagation of the matter wave function is considered simultaneously with an ETRS scheme and if both, matter and electromagnetic fields are coupled self-consistently in every time-step. We provide more details for such a self-consistent coupling in section V.7.

The time step parameter $\Delta t$ has to be chosen such that the propagation remains stable and accurate. Both criteria are not defined strictly, but it is possible to reduce the value $\Delta t$ until convergence is reached. A well-known criterion for the time step $\Delta t$ is given by the Courant-Friedrichs-Lewy (CFL) condition Courant et al. (1928, 1967)

[TABLE]

and depends basically on the grid spacings $\Delta x$, $\Delta y$, $\Delta z$ for the three-dimensional grid, which is taken to be equidistant in each dimension, and the Courant number $S_{\mathrm{CFL,max}}$. The Courant number $S_{\mathrm{CFL,max}}$ varies for different propagation methods. In case of FDTD in three dimensions, $S_{\mathrm{CFL,max}}$ is equal to one Courant et al. (1928, 1967). We have performed convergence tests of our Riemann-Silberstein time-evolution and have found that we can chose a Courant number slightly larger than one.

To summarize this section, we have introduced a formulation of Maxwell’s equations in Riemann-Silberstein form, which allows to use time-evolution techniques that have been developed for the time-stepping of the time-dependent Schrödinger equation. Since the same grid is used for the electric and the magnetic fields, this simplifies the implementation, in particular the coupling to matter, and allows for larger time steps compared to FDTD since higher order finite-difference discretizations can be employed. Choosing the Hamiltonian-like kernel $\hat{\mathcal{H}}$ and the microscopic paramagnetic current $J_{\mathrm{pmc}}^{k}$, diamagnetic current $J_{\mathrm{dmc}}^{k}$ and magnetization currents $J_{\mathrm{mc}}^{k}$ as inhomogeneity, allows to evolve the Maxwell fields coupled to matter. On the other hand, replacing $\hat{\mathcal{H}}$ with an integral kernel for a linear medium $\hat{\mathcal{H}}_{\mathrm{lm}}$, and using a free classical current density, allows to simulate with the same implementation macroscopic Maxwell fields. In principle it is also possible to combine both cases in the same real-time simulation. Such a combined microscopic-macroscopic propagation can be used for example to describe a molecule inside an optical cavity geometry or nanoplasmonic systems close to an interface.

V.3 Maxwell boundary conditions

For the Riemann-Silberstein implementation, it is very important to have flexible boundary conditions for the electromagnetic fields. Since we are dealing with a finite simulation box for the radiation fields it is of course necessary to have absorbing boundaries, such that outgoing electromagnetic fields propagate without any reflections at the box boundaries. While this is a standard procedure in FDTD simulations, we emphasize here that such absorbing boundaries effectively allow to turn our coupled light-matter system into an open quantum system from first principles. No artificial bath degrees of freedom have to be introduced as commonly done in the description of open quantum systems. Another class of boundary conditions, that arise from the different length scales of typical wavelengths for light and matter, are incident wave boundary conditions. Such boundaries allow to feed electromagnetic waves with much larger wavelength into the microscopic simulation box. Finally, the boundaries can also serve as electromagnetic detectors. Provided the charge densities and currents of the matter system have decayed sufficiently and are effectively zero in the boundary region, the propagation of the electromagnetic fields from there on then only corresponds to a propagation in vacuum. In other words, whatever arrives in the boundary region would propagate to the far field and contributes to what can be measured in the far field by a detector.

In the following, we discuss the different types of boundary conditions relevant for our coupled light-matter simulations.

V.3.1 Boundary regions

To properly define different boundary conditions, we separate the simulation box into two regions, the inner free Maxwell-propagation region and the outer boundary region with specified Maxwell’s equation to fulfill the appropriate simulation condition. Such a region splitting is shown for a two-dimensional cut of the three-dimensional simulation box in Fig. 1. The outer box limits are determined by $L_{u}$ for direction $u\in(x,y,z)$, whereas the boundary region is limited by $b_{u}$. We note that $L_{u}$ and $b_{u}$ are always positive and the box center is always located at the Cartesian origin. The total box dimension in each direction is $-L_{u}$ to $+L_{u}$, and the inner borders of the free propagation region are $-b_{u}$ and $+b_{u}$. Hence, the area between $b_{u}<|u|<L_{u}$ describes the boundary region which contains different conditions to fulfill the corresponding simulation system which we consider in the following.

V.3.2 Absorbing boundaries by mask function

A simple, robust and easily applicable method for absorbing boundaries can be reached by multiplying the Riemann-Silberstein vector after each time step by a mask function which decreases the electromagnetic fields at the boundaries. A proper mask function $f_{\rm{mask}}^{1\mathrm{D}}(u)$ in direction $u\in(x,y,z)$, which is already implemented in Octopus and satisfies the condition of damping the fields smoothly in the boundary region, is given by Andrade et al. (2015)

[TABLE]

For a three-dimensional simulation box, the corresponding mask function $f_{\rm{mask}}^{3\mathrm{D}}(x,y,z)$ is factorized into three one-dimensional mask functions for each direction and therefore takes the form

[TABLE]

The only parameter that improves the absorbing effect of the mask function is the boundary width $L_{u}-b_{u}$ in the corresponding direction. In case of relatively large absorbing boundary regions and small outgoing electromagnetic field amplitudes compared to the inner box fields, absorbing boundaries in terms of mask functions are a simple and effective method for simulating open Maxwell systems. Whenever larger electromagnetic field amplitudes have to be absorbed, or when using larger padding regions for the mask function become computationally too costly, it becomes advantageous to use a perfectly matched layer, which we discuss next.

V.3.3 Absorbing boundaries by perfectly matched layer

A more accurate method for open Maxwell systems is the perfectly matched layer (PML) absorbing boundary condition. We have implemented such a PML analogous to the Berènger layer for the Yee finite difference time domain algorithm Taflove and Hagness (2005); Bérenger (2007), but now modified for the Riemann-Silberstein Maxwell propagation. The basic idea of the Berènger layer is to complement Maxwell’s equations with an artificial lossy layer, which is described by a non-physical electric conductivity $\sigma_{\mathrm{el}}$ and non-physical magnetic conductivity $\sigma_{\mathrm{mag}}$. These conductivities are defined such as to yield minimal reflections at the boundaries, but have no physical meaning otherwise. The loss due to the conductivity parameters is linear in the corresponding $\vec{E}(\vec{r},t)$ and $\vec{B}(\vec{r},t)$, so that Faraday’s and Ampère’s law without current density take the form

[TABLE]

[TABLE]

We note here, that the PML is in principle not restricted to vacuum conditions but also valid for other homogeneous dielectric conditions. In other words, the Riemann-Silberstein PML that we have implemented also works in a linear medium, but with the constraint that $\epsilon$ and $\mu$, and consequently $c$ are constant at the border and inside the boundaries. Transforming Eq. (78-79) to frequency domain, where $\tilde{F}$ denotes the Riemann-Silberstein vector in frequency space, leads to the underlying Riemann-Silberstein Maxwell’s equation for the absorbing layer

[TABLE]

where the first term on the right-hand side describes the curl operation. The principle of a PML is to propagate the respective field components in the absorbing boundary region which are necessary for a correct propagation inside the free Maxwell region, and to damp the remaining components without causing strong reflections back into the free Maxwell region. For this purpose, Berènger’s method splits up Maxwell’s equations for each direction in two equations which form the basis for the so-called split PML Berenger (1996); Bérenger (2007); Taflove and Hagness (2005). The field component in $k$-direction is split into one component for $l$ and one for $m$ with $k\neq l\neq m$, so that the vector k component $\tilde{F}_{\pm}^{k}$ is given by

[TABLE]

The field component $\tilde{F}_{\pm}^{k}$ is split such that the $\tilde{F}_{\pm}^{k,(l)}$ part is responsible for the field propagation parallel to direction $l$ and accordingly $\tilde{F}_{\pm}^{k,(m)}$ parallel to direction $m$. In other words, there are two separate propagations which simulate only the free propagation along the corresponding direction. Thus, one propagation could be where the field enters the PML region while the other part is still in the free propagation box. The damping of the fields is applied by the electric and magnetic conductivities $\sigma_{\mathrm{el}}$, $\sigma_{\mathrm{mag}}$ which are artificially modified and depend now on the splitted direction, i.e., $l$ direction for $\tilde{F}_{\pm}^{k,(l)}$, not the field direction. In addition, each equation only contains one part of the two curl terms. Applying all these considerations to the six components of the Maxwell’s equations in Riemann-Silberstein form yields twelve relations for the PML. Explicitly, the two equations for the $x$ component in equation (135) are

[TABLE]

[TABLE]

for the $y$ component

[TABLE]

[TABLE]

and for the $z$ component

[TABLE]

[TABLE]

Analogous to Berènger’s split field PML derivation for the Yee-Algorithm, we want to include also in our case the frequency $\omega$ and the electric and magnetic conductivity $\sigma_{\mathrm{el}}$ and $\sigma_{\mathrm{mag}}$ in a factor multiplied by the corresponding split field Berenger (1996); Bérenger (2007); Taflove and Hagness (2005) before we recombine the two split field equations. Using the two factors

[TABLE]

[TABLE]

the system of the split equations in Eqs. (137 - 142) can be rearranged equivalently to

[TABLE]

[TABLE]

for the $x$ component, and

[TABLE]

[TABLE]

for the $y$ component, and

[TABLE]

[TABLE]

for the $z$ component. Finally, adding each of the two Eqs. (145-146), (148-147), and (149-150) using Eq. (136) yields the PML equations in frequency domain for the Riemann-Silberstein representation

[TABLE]

In our PML implementation for simplicity, we do not introduce new correction terms in $\tilde{\eta}^{(l)}$ or $\tilde{\xi}^{(l)}$ to improve the PML and to reduce low-frequency reflections like it is commonly applied for the Yee algorithm Berenger (1996); Bérenger (2007); Taflove and Hagness (2005). While such extensions are possible in future refinements of our implementation, we found the simple form without correction terms already to provide good absorbance at the boundaries. By back transforming Eq. (151) from frequency domain into time domain, we arrive at

[TABLE]

[TABLE]

The electric conductivity $\sigma_{\mathrm{el}}$ and the magnetic conductivity $\sigma_{\mathrm{mag}}$ have to be chosen such that the reflection becomes minimal. As is well-known in FDTD Berenger (1996); Bérenger (2007); Taflove and Hagness (2005), the relation between the electric conductivity $\sigma_{\mathrm{el}}$ and the magnetic conductivity $\sigma_{\mathrm{mag}}$ to minimize the reflection coefficient has to obey

[TABLE]

at the border between the free Maxwell simulation box and the absorbing boundaries. Using this relation between the two conductivities, it is convenient to use only one conductivity with $\sigma=\sigma_{\mathrm{el}}$, and the updated forms of the expressions $\zeta^{(l)}(t)$, and $\xi^{(l)}(t)$ are

[TABLE]

[TABLE]

As a result, the back transformation of Eq. (151) becomes

[TABLE]

which contains several convolutions in time. Whereas the first convolution on the right-hand side in (157) is simply

[TABLE]

the remaining convolutions are explicitly given by

[TABLE]

[TABLE]

This completes the construction of the PML for our Riemann-Silberstein formulation. Since we have already illustrated how to include a linear medium in the Riemann-Silberstein time-evolution in the previous sections, it becomes now straight forward to combine the PML with our existing implementation. Adding the adequate PML expressions to the integral kernel for the Riemann-Silberstein Hamiltonian $\hat{\bar{\mathcal{H}}}$ in Eq. (94), we arrive at a propagation scheme with perfectly matched layer boundaries

[TABLE]

with

[TABLE]

where the 6x6 PML matrix $\hat{\bar{\mathcal{G}}}(\vec{r},t,\tau)$ is given by

[TABLE]

The left factor of the second Kronecker product in Eq. (163) has entries in the off-diagonal, and therefore the two Riemann-Silberstein vectors $\vec{F}_{\pm}$ always couple in the PML region.

In principle, the PML terms in Eq. (161) have to be calculated for each time step, which massively increases computational cost. However, taking a closer look at Eqs. (159) and (160), we notice that the two functions $\zeta^{(k)}(t-\tau)$ and $\xi^{(k)}(t-\tau)$ contain exponential factors. Therefore it is possible to obtain a rather accurate approximation of the terms by using a recursive-convolution method Luebbers and Hunsberger (1992) with finite time steps $\Delta t$. The recursive-convolution method allows to express integrals of the form

[TABLE]

in terms of

[TABLE]

where we have taken $t=m\Delta t$. For finite yet sufficiently small time steps $\Delta t$, it can be assumed that the function $f_{2}(\tau)$ in the last integral term on the right-hand side of equation (165) is constant. This allows to take $f_{2}((m-1)\Delta t)$ outside of the integral. In the next step, we substitute $\alpha=\sigma$ and the functions $f_{1}(t)$ and $f_{2}(t)$ with the corresponding ones in Eq. (159) and Eq. (160). Therefore, we obtain $f_{1,+,l}(t),f_{2,+,l}(t)$ for Eq. (159) and $f_{1,-,l}(t),f_{2,-,l}(t)$ for Eq. (160) with

[TABLE]

[TABLE]

However, the above recursive convolution applied to Eqs. (161) and (163) does not allow to express the term $\hat{\bm{\mathcal{G}}}\bm{\mathcal{F}}$ in Eq. (161) as a matrix vector multiplication with an approximated matrix $\hat{\mathcal{G}}$ and vector $\mathcal{F}$. Nevertheless, it is possible to replace the whole $\hat{\bm{\mathcal{G}}}\bm{\mathcal{F}}$ term by a 6x6 dimensional matrix, denoted as $\tilde{\mathcal{G}}$

[TABLE]

The matrix $\tilde{\mathcal{G}}$ contains four 3x3 dimensional matrices, which are defined recursively and depend on the current $t=j\Delta t$ and the previous time $t^{\prime}=(j-1)\Delta t$. These recursive matrices are given by

[TABLE]

The auxiliary variables $a_{k}$ and $b_{\pm,k}$ result from the last line of Eq. (165). Adapting them to the corresponding PML equation yields

[TABLE]

[TABLE]

Collecting all steps, we can express the Maxwell Riemann-Silberstein Hamiltonian $\hat{\mathcal{H}}_{\mathrm{PML}}$ from Eq. (161) in discretized form

[TABLE]

With this definition it is then straightforward to insert the above PML expression $\hat{\bar{\mathcal{H}}}_{\mathrm{PML}}$ in the Maxwell propagator $\bm{\hat{\mathcal{U}}}$ of the numerical propagation equations in Eq. (126) or Eq. (127) to enable the simulation of open quantum systems via PML absorption.

In the last step, we have to determine the conductivity $\sigma_{(u)}$ adequately to get an optimal PML. In FDTD, several useful profiles for the conductivity $\sigma_{(u)}$ were found and we have chosen for our Riemann-Silberstein PML the FDTD polynomial grading profile which has the form Taflove and Hagness (2005)

[TABLE]

with direction coordinate $u\in(x,y,z)$ where $b_{u}$ and $L_{u}$ denote the corresponding boundary dimensions in Fig. 1. The last variable $\sigma_{\mathrm{max}}$ for the grading profile is determined by

[TABLE]

where the tolerated reflection error for normal angle incidence equal to zero can be set manually. The parameter $q$ for the grading profile is usually set in a range between $2.0\leq q\leq 4.0$ to get numerically sufficient absorbance.

V.3.4 Incident waves boundaries

A very frequently found experimental situation corresponds to cases where the matter system is a molecule or nanoparticle on the scale of a few tens or hundreds of Ångtröms, while the incident laser fields are typically in the infrared, optical, or ultraviolet range. This corresponds to a spatial extension of the optical waves which is a few orders of magnitude larger than the matter system. Simulating both, radiation and matter, in the same simulation box, would therefore require boxes with an unfeasible size. In addition, the used light pulse can frequently be approximated with a mathematically closed description such that the time-evolution is known analytically. For such situations, it is not necessary to chose a large Maxwell simulation box such that the external light signal is completely inside the box. A convenient method to simulate such electromagnetic waves in our Riemann-Silberstein simulation box can be obtained by using the boundary region to match with the outside free evolution of, e.g., a laser pulse, by updating the values of the grid points according to the Dirichlet boundary conditions of the free evolution. In Fig. 2, we illustrate such a scenario. We show in a 2D cut the analytically calculated outer frame around the simulation box. In the interior region a Gaussian shaped pulse envelope propagates parallel to its wavevector $\vec{k}$, which is performed numerically. The transition between analytical boundary and numerical electromagnetic wave in the interior is seamless.

In general, an arbitrary shaped analytical wave can be obtained by a superposition of different linearly independent waves with

[TABLE]

where the $i^{th}$ wave is represented by its wavevector $\vec{k}_{(i)}$ and frequency $\omega_{(i)}$ and by a Riemann-Silberstein vector $\mathcal{F}_{\mathrm{pw},(i)}$ as initial vector. We select the width of the boundary region where this incident analytical wave is prescribed as Dirichlet boundary condition according to the number of the grid points that are used for the stencil of the finite-difference discretization. In this way artefacts at the interface of boundary region and interior region can be avoided.

We emphasize here that the incident plane-wave boundary condition amounts to simulating an open quantum system since energy enters the system through the analytically prescribed boundaries. Time-reversal symmetry does not hold for open systems, especially in the presence of magnetic fields, and consequently the ETRS propagator in Eq. (127) does not hold. However, we assume in the present work that the full coupled Hamiltonian stays approximately time-reversal since the main breaking of the symmetry arises if we consider the magnetic field propagation without any back reaction of the matter.

The incident waves boundaries give us also the opportunity to simulate pump-probe experiments if we propagate two different signals with arbitrary angles and time delay that hit the molecule and calculate the resulting electromagnetic field.

V.3.5 Incident waves boundaries plus absorbing boundaries

The incident waves boundaries that we just introduced allow to simulate signals that enter the simulation box. On the other hand, the PML dampens all outgoing electromagnetic fields. However, for most purposes it is desirable to combine both boundary conditions. Due to the analytical behavior of the incident waves, we know the incoming fields for all times. In addition, the outgoing electromagnetic signals should not cause any reflection at the boundaries. In the previous section, we have shown how the PML looks like only for the absorbing boundary condition. Since the Maxwell’s equations for the incoming fields is just linear we can add it to the matter-coupled (internal) field that is absorbed at the boundaries. At the same time, if we have given the total electromagnetic field (internal and incoming fields) we can easily subtract the incoming field in the whole simulation box to determine the internal fields. This allows us to apply the PML only to the internal fields of the coupled light-matter system. To combine incident waves with absorbing boundaries, we therefore split the boundary region into two regions: an outer region for the incident waves, and an inner one for the PML as it is shown in Fig. (3).

A time step is then performed such that first the freely propagating wave that arises from the incident wave is subtracted from the numerically propagated Riemann-Silberstein vector. Next, the PML is applied to the remaining field. In the last step, the freely propagating wave is added again to the field. We note, that there are two options to compute the values of the incident field. As first option, one uses simply the analytical wave values for the current time step. The second option requires a second auxiliary propagation, where the incident wave is propagated on the numerical grid, however, without any coupling to the matter. In other words, in this case, the free wave propagation is calculated by the discretized Maxwell time-evolution operator. This method avoids numerical reflection artefacts at the boundary due to the fact that there are always small numerical discrepancies between the analytical wave and a numerically calculated one.

V.4 Time-evolution operator for Kohn-Sham orbitals

Since various time-stepping schemes for the time-dependent Kohn-Sham equations are already well established and have been extensively discussed in the literature Castro et al. (2004); Gómez Pueyo et al. (2018), we here only briefly summarize the basic equations and introduce the corresponding notation.

According to Eq. (68), the one-particle MPKS Hamiltonian takes the form

[TABLE]

where the canonical momentum is given by

[TABLE]

Summing these one-particle non-interacting Hamiltonians $\hat{h}_{\mathrm{MPKS}}^{(n)}(\vec{r}_{i})$, gives the non-interacting many-particle MPKS Hamiltonian for the species $n$, i.e.,

[TABLE]

The time evolution of the Kohn-Sham orbitals $\phi_{(n,i)}$ from starting time $t_{0}$ to time $t$ is given by

[TABLE]

where $\hat{u}_{\mathrm{MPKS}}^{(n)}(t,t_{0})$ denotes the corresponding time-ordered MPKS time-evolution operator for the species $n$

[TABLE]

Based on this, the Kohn-Sham wave function evolves in time according to

[TABLE]

where

[TABLE]

and $N^{(n)}$ denotes the number of occupied Kohn-Sham orbitals for species ${(n)}$. Note that the evolution operators do not need to be (anti-)symmetrized, since the symmetry of the initial state is kept.

For our coupled Maxwell-Kohn-Sham system, we want to use for both, matter and radiation, consistent time-evolution techniques. We therefore select, as in the case of the Riemann-Silberstein time-stepping, the ETRS time-evolution operator to propagate the Kohn-Sham orbitals. The numerical ETRS time-evolution equation for the $(m+1)\Delta t$ time step of the Kohn-Sham orbital $\phi_{(n,i)}(\vec{r},(m+1)\Delta t)$ is given by

[TABLE]

where we have introduced the ETRS time-evolution operatorCastro et al. (2004)

[TABLE]

Similar to the Maxwell time-evolution, the time step parameter $\Delta t$ has to be selected to yield a stable and accurate propagation. However, in contrast to the Maxwell system, there is no CFL criterion for the Kohn-Sham evolution since the speed of matter waves in our non-relativistic setting is not capped by the speed of light. Nevertheless, since we restrict our considerations to low-energy problems where the non-relativistic approximation is justified, by construction the Kohn-Sham orbitals evolve much slower than the speed of light. Hence, the maximum $\Delta t$ is in most cases much larger than the one for the Maxwell fields (for more details see the different case studies presented below).

V.5 Full minimal coupling and multipole expansion

In this section we discuss the different levels of theory that can be used to couple light and matter.

Full minimal coupling extends the Dirac equation to include the coupling to an electromagnetic field, taking both, Lorentz- and gauge-invariance into account. Although we partially violate the Lorentz-invariance property in our non-relativistic limit that leads to the MPKS approach based on the Hamiltonian in Eq. (176), it still represents an accurate approximation provided the kinetic-energies of the particles are well below relativistic values. In that case higher-order relativistic corrections to the Pauli equation become important. Such semi-relativistic version of the Pauli-Fierz Hamiltonian exist and can then be used instead. But for most applications we envision, the MPKS approach should be sufficiently accurate. Let us in the following see how we get from this full minimal-coupling description to simplified couplings that are often employed and that we want to use later to investigate the influence of these approximations on physical processes and observables.

As first step, we separate the Hamiltonian of Eq. (184) into a kinetic Hamiltonian $\hat{H}_{\mathrm{kin}}$ and an interacting Hamiltonian $\hat{H}_{\mathrm{int}}$ which includes Maxwell and matter variables

[TABLE]

where the kinetic piece is given by

[TABLE]

and the light-matter coupling is contained in

[TABLE]

The total vector potential $a_{\mathrm{KS}}^{\mu}=a^{\mu}+A^{\mu}+a^{\mu}_{\rm{xc}}$ in Eq. (186) is determined by the Riemann-Silberstein vector via Eq. (91) or (92). We note here, that the total vector potential, especially the scalar potential component $a^{0}_{\mathrm{KS}}$, in principle includes all electronic and nuclear potentials.

In many applications, the full minimal coupling is expensive to calculate or not needed since the length scales of matter and radiation differ vastly. Therefore, the minimal coupling is often approximated by a multipole expansion using the electric and magnetic fields variables. As is well-known, the ubiquitous electric dipole approximation is equivalent to the lowest order term of this multipole expansion. Since this type of multipole approximations are ubiquitous in quantum physics it is very interesting to investigate how observables change order-by-order. In the following, we briefly summarize the derivation of the multipole expansion based on the Power-Zienau-Woolley transformation (cf. chapter 5.2 of Ref. Loudon (1988)) and adapt it to the present MPKS case. As a first step, we introduce the polarization

[TABLE]

In Coulomb gauge with $\vec{\nabla}\cdot\vec{a}_{\mathrm{KS}}=0$, the vector potential is always transversal and hence the unitary Power-Zienau-Woolley transformation $\hat{U}_{\mathrm{PZW}}$ is defined by

[TABLE]

By then using the MPKS equation, however, transformed with the unitary operator $\hat{U}_{\rm{PZW}}$ for the wave function $\Phi^{\prime(n)}=\hat{U}_{\mathrm{PZW}}^{-1}\Phi^{(n)}$, leads together with

[TABLE]

as well as an extra term due to $\partial_{t}\hat{U}_{\mathrm{PZW}}^{-1}(t)$ to

[TABLE]

In this form it becomes easy to perform a multipole expansion for the electric dipole term

[TABLE]

the magnetic dipole term

[TABLE]

and the electric quadrupole term

[TABLE]

all expanded around a point $\vec{r}_{0}$, which can be chosen in good approximation either as center of mass or center of charge of the matter system. We note that we here used a multipole expansion based on classical (Kohn-Sham) fields. Without time-dependence in the classical-field case no electric transversal field can appear. If we instead performed the Power-Zienau-Woolley transformation on the quantized level (so the field in Eq. (188) is a time-independent operator) two major differences appear: First, instead of with electric fields we work with displacement fields and a novel term, the dipole self-energy term appears Loudon (1988); Rokaj et al. (2018). This term is physically and mathematically necessary. Secondly, since the fields are now (unbounded) operators, the Taylor expansion in its usual form is not applicable. Thus in the case of a Power-Zienau-Woolley transformation with quantized fields extra care needs to be taken Rokaj et al. (2018); Spohn (2004).

For our present MPKS implementation, we have neglected higher order multipole terms that contain non-linear orders of Maxwell-field variables and the Stern-Gerlach term. Whereas the magnetic dipole term depends on the total magnetic field, the electric dipole and quadrupole terms depend only on the transverse component of the electric field. Consequently, we have to decompose the Riemann-Silberstein vector into its transverse and longitudinal components. In general, the Helmholtz-decomposition for the Riemann-Silberstein vector is

[TABLE]

The first term on the right-hand side of Eq. (195) can be computed efficiently by a Poisson solver since it is the solution of the Poisson equation. The second integral in Eq. (195) is a surface integral which is necessary, if the Riemann-Silberstein vector does not vanish at the simulation box boundaries, such as for, e.g., periodic systems. Since the Riemann-Silberstein vector in the multipole Hamiltonian does only depend on the expansion point $\vec{r}_{0}$ of the multipole expansion and its corresponding Riemann-Silberstein vector inside the box, it is sufficient to calculate the values of the surface integral only for the stencil points that correspond to the expansion center $\vec{r}_{0}$ of the multipole expansion. This reduces the computational cost for the boundary term significantly.

V.6 Multi-scale implementation

In the previous sections we have already seen that the time-evolution for matter and electromagnetic fields can be expressed mathematically in a very similar way. While both share a similar Schrödinger-type form, the physical parameters for many applications of the theory imply rather different length scales for light and matter. For example, when molecular systems interact with infrared, optical, or ultraviolet laser pulses, the matter wave functions are mainly localized in a relatively small volume with possibly rather strong fluctuations. In contrast, the Maxwell fields that result from typical experimental resonators are localized on the scale of the pulse envelope and are often smoother than the matter wave functions. For optical lasers, the length scales of radiation and matter waves differ by about two or three orders of magnitude. To handle such situations we use a finer grid for the matter wave functions and coarser grids for the Riemann-Silberstein vector. Our MPKS implementation in Octopus is however not restricted to such laser-molecule interaction applications. In the following, we introduce the possible combinations of matter and Maxwell grids that can be selected for MPKS simulations with Octopus.

V.6.1 Multi-grid types

In Fig. 5, we illustrate the different possibilities to combine matter and Maxwell grids. We limit the graphical illustration here to the case when the matter grid spacing is smaller or equal to the Maxwell grid spacing. Our implementation is however not limited to these cases and can also be used for the reverse situation which is suitable to simulate, e.g., hard x-rays.

The grids a) and b) in Fig. 5 correspond to simulation boxes with identical dimensions for both systems. In general, like in a), the Maxwell grid points do not necessarily have to lie on top of a matter grid point. For both grid types a) or b), it is not possible to use the field values at given grid points directly in the respective coupling terms of the propagation equations for light and matter. Instead, as is used routinely for multigrid methods in literature, prolongation and interpolation schemes have to be used. For example, the finer matter grid points have to be grouped in clusters which map to the next nearest Maxwell grid point and different weighted mean schemes can be applied to yield the coupling value for that Maxwell point. Vice versa, a Taylor series expansion of the Maxwell values can be used to determine the coupling values for the matter points. These grid types are well suited to simulate periodic systems.

The two schemes in Figure Fig. 5 c) and d) show matter grid types, in the first case without common grid points, and in the second case with common grid points, but compared to the previous grids in a) and b) with smaller dimensions for the matter grid than for the Maxwell one. As before, the values for the corresponding coupling terms have to be calculated by weighted mean and interpolation. Such grid combinations can describe efficiently bound molecules and nanoparticles, especially when focussing on electromagnetic far-fields.

If near field effects are of interest, where the electromagnetic field fluctuations correlate strongly with the matter wave functions, it is a good choice to select the same grid spacings for both grids and to place matter and Maxwell grid points on top of each other as shown in Fig. 5 e) and f). In this case, the values for both respective coupling terms can be obtained directly from the field values at the respective grid point. Besides near-field simulations, the grid type f) with larger Maxwell grid dimensions is suited to study the onset of the electromagnetic far-field and allows to define electromagnetic detectors at the box boundaries.

Finally, a further grid combination is illustrated in Fig. 5 g), where the matter grid is chosen much finer than the Maxwell grid. Only one Maxwell grid point lies in the middle of the matter grid. Here, it is assumed that the Maxwell field is approximately constant for all matter grid points. Vice versa, the coupling value for the Maxwell grid is obtained by the mean value of all matter points.

V.6.2 Multi-scale in time

So far we have focussed on different length scales and the related choices for spatial grid spacings and box dimensions. In this section we turn our attention to the different time scales and the related choices for numerical time steps for the matter and Maxwell propagation.

While it is in principle possible to chose identical time steps for the matter and Maxwell propagation, given the constraint that both propagators have to run stable, this is for most physical cases not the most efficient choice. The reason for this can be seen directly in the Maxwell Hamiltonian-like form in Eq. (87). The gradient in this expression is multiplied with the speed of light $c_{0}$ (roughly 137 in atomic units). This factor imposes the speed for the electromagnetic waves on the grid. On the other hand, the matter Hamiltonian in our non-relativistic Pauli limit is lacking the factor of $c_{0}$, yielding a much smaller spectral range of the maximum and minimum eigenvalue of the Hamiltonian for a given grid and stencil, and hence a much slower wave motion. As consequence of the ”fast” photon motion compared to the ”slow” motion of matter, in the Maxwell case a much smaller maximum time step $\Delta t_{\mathrm{Mx,max}}$ has to be selected compared to the maximum time step of the matter $\Delta t_{\mathrm{KS,max}}$.

Such a situation of different physical time scales arises already in electron-nuclear dynamics, where the large nuclear mass leads to a rather slow motion of the nuclei compared to the faster motion of the lighter electrons. Coupling the propagation of the electromagnetic fields also with the electron-nuclear dynamics is adding a third time-scale to the problem. In our numerical time-stepping scheme, we exploit the different time-scales explicitly to increase the computational efficiency. Our simulations have shown, that the coupled propagation of nuclei, electrons, and electromagnetic fields keeps relatively accurate, stable and converged, if we perform several Maxwell propagation steps $\Delta t_{\mathrm{Mx}}$ in-between the Kohn-Sham propagation steps for the electrons and Ehrenfest steps for the nuclei Xavier Andrade (2009). For convenience, we select the Kohn-Sham time step $\Delta t_{\mathrm{KS}}<\Delta t_{\mathrm{KS,max}}$ as the basic time step parameter for the coupled MPKS system, and the number $N_{\mathrm{Mx-steps}}$ of intermediate Maxwell steps is automatically chosen such that

[TABLE]

where $\Delta t_{\mathrm{Mx,CFL}}$ denotes the Courant time step for the Maxwell grid given in equation (128). Performing these intermediate Maxwell steps, assumes that the intermediate Maxwell propagation between $m\Delta t_{\mathrm{KS}}$ and $(m+1)\Delta t_{\mathrm{KS}}$ does not affect the matter propagation significantly, which means that the current density for the corresponding $i^{th}$ intermediate time step is well approximated by the linear expansion

[TABLE]

The ETRS time-evolution equation for the $i^{th}$ intermediate step then takes according to Eq. (123) the form

[TABLE]

with

[TABLE]

To reduce the computational cost even further, we can assume in most cases that the inhomogeneity terms in equation (123) are approximately constant for all intermediate time steps during the time interval $\Delta t_{\mathrm{KS}}$. Thus, we use in this case for the current density $\mathcal{J}$ in equation (198) the arithmetic mean of $\mathcal{J}(\vec{r},m\Delta t)$, and $\mathcal{J}(\vec{r},(m+1)\Delta t)$ which reduces the amount of necessary computational expensive $\bm{\hat{\mathcal{U}}}$ and $\hat{\mathcal{U}}_{(0)}$ operations. Clearly, these simplifications introduce approximations to the full time-evolution and convergence has always to be checked for a given application.

V.6.3 Parallelization strategies

For a speedup of the time-evolution of the Riemann-Silberstein vector, we employ a domain parallelization for the real-space grid. Since we have implemented our propagation scheme in Octopus by using the existing gradient operations that have been optimized for matter degrees of freedom, we immediately inherit the existing domain parallelization of the code Castro et al. (2006) also for the Maxwell case. A parallelization in ”states” or ”orbitals/k-points” as is used for matter wavefunctions, is not so effective for the Maxwell case, as there are always only six ”orbitals” to propagate, namely the six vector components of the Riemann-Silberstein vector. We have therefore restricted ourselves for the calculations in the present work solely to the domain parallelization of the Maxwell grid.

V.7 Predictor-corrector scheme

So far we have described separately the propagation schemes for matter and electromagnetic fields as well as the coupling Hamiltonian. In the present section, we discuss the predictor-corrector method which we employ to enforce a self-consistent propagation of the subsystems.

V.7.1 Forward coupling

In most studies in the literature light-matter coupling is restricted to forward Maxwell-matter coupling: the electromagnetic fields influence the matter degrees of freedom, but the matter currents are not allowed to influence the propagation of the electromagnetic fields. This situation is illustrated on the left-hand side of Fig. 6. In the forward coupling case, the external Maxwell field propagates without any perturbation by the matter and is calculated separately, either analytically or numerically. Looking at the time-evolution operator in equation (183), we note that the operator $\hat{U}_{\mathrm{MPKS}}^{\mathrm{ETRS},(n)}((m+1)\Delta t_{\mathrm{KS}},m\Delta t_{\mathrm{KS}})$ depends on the Hamiltonian operator $\hat{H}_{\mathrm{MPKS}}^{(n)}((m+1)\Delta t_{\mathrm{KS}})$ at the future time $t=(m+1)\Delta t_{\mathrm{KS}}$. This future Hamiltonian is not only determined by the external Maxwell fields but also by the motions of the ions and electrons and their interactions. Therefore, it is necessary to apply a predictor corrector cycle for the matter propagation. In a first step, the future Hamiltonian is estimated by an extrapolation Marques et al. (2003) and the calculated time propagation returns an estimated Kohn-Sham potential which is again used for an updated extrapolation of the Hamiltonian. These steps are repeated until the absolute value of the variance of two subsequently Kohn-Sham potentials falls below a small threshold value. For our calculations, we set a threshold value of $1e^{-6}$ in atomic units for the potential variance and adjust the time step $\Delta t_{\mathrm{KS}}$ for the propagation so that the matter system is converged in at least two iterations if the system is only disturbed very weakly by the external field. During the full run and stronger perturbations, we notice that the number of iterations is barely larger than five.

V.7.2 Forward and backward coupling

The back-reaction of matter on the Maxwell fields appears due to the current density in equation (36) (in the MPKS reformulation due to its expectation value), which is caused by the motion of matter. The three electronic current types, paramagnetic, diamagnetic and magnetization current, as well as the ionic current influence the Maxwell propagation equation (123). The influence of the paramagnetic current, the magnetization current and the external diamagnetic current as well as optional external currents result directly in the inhomogeneity $\mathcal{J}$ term in equation (100). The internal diamagnetic current implicitly influences the time-evolution due to the modified Maxwell time-evolution operator $\bm{\hat{\mathcal{U}}}$ for this case given in (115). The full forward and backward coupling scheme is shown on the right-hand side of Fig. 6.

In a fully self-consistent scheme, both systems and accordingly their time-evolution propagation equations in (123) and (178) couple to each other. First we apply the extrapolation of the future matter Hamiltonian to get a prediction for the Kohn-Sham orbitals. These orbitals and the initial ones give us the necessary current density which couples to the Maxwell fields. Using the first predicted updated current density at time $(m+1)\Delta t_{\mathrm{KS}}$ leads to an updated Riemann-Silberstein vector. At this point, the predictor-corrector loop restarts by updating the Kohn-Sham orbitals but now with a corrected matter Hamiltonian which includes the updated Riemann-Silberstein vector. As a consequence, the previously predicted variables get a correction closer to the values which make the coupled system self-consistent. We additionally check the consistency of the Maxwell fields by comparing the Maxwell energy inside the simulation box for two successive updated Riemann-Silberstein vectors. For consistency we use the same threshold value $1e^{-6}$ as for the matter case. Additionally, we chose the system propagation time $\Delta t$ such that the predictor-corrector step iterates at least two times until the self-consistency thresholds are fulfilled for weak perturbations. Again, from our experience the number of iterations for strong perturbation periods should not be larger than five steps. Otherwise the time steps have to be reduced.

V.8 Validation and comparison to FDTD

To validate our Maxwell implementation, we compare in this section our results directly with the MIT Electromagnetic Equation Propagation (MEEP) Oskooi et al. (2010) package, which is a common simulation software for electromagnetic field propagation. The implemented Maxwell field propagation in MEEP is based on the FDTD method and the Yee-algorithm Yee (1966). The underlying electromagnetic simulation grid is split into two grids shifted by half of the grid spacing of the corresponding direction. As a consequence, the requested spatial and time derivative points for the propagation equation are in the middle of two adjacent grid points. Therefore, the center finite-difference method leads to the first-order derivative equation (here for simplicity discussed in one dimension)

[TABLE]

In case of the Yee grid, there are no grid points at the derivative points $x_{0}$ but next to it at $x_{0}-\Delta x/2$ and $x_{0}+\Delta x/2$. In our finite-difference discretization in Octopus, we also use the center finite-difference formula to get first-order derivatives, but the derivative points lie always on top of a grid point. Thus, the derivative equation takes the form

[TABLE]

which means that the error of $f^{\prime}(x_{0})$ is smaller for the Yee-algorithm if we consider the same order terms of the finite difference method. However, the MEEP finite difference stencil operation is always of order two whereas the Octopus stencil order can be set to higher orders to obtain better accuracy for the derivative operators.

As a test scenario, we simulate a spatial and temporal shaped external current density inside the simulation box, which is prescribed analytically, and plot several relevant physical variables. For this reference calculation, all physical units in Octopus are set equal to MEEP internal units. The spatial and time dependent current density that we impose externally takes the form

[TABLE]

We select a cubic simulation box of size $20.0$ and a PML boundary region width of $4.0$ in each direction. Referring to our sketch of the simulation box in Fig. 4, the corresponding parameters are $L_{x}=L_{y}=L_{z}=44.0$ and $a_{x}=a_{y}=a_{z}=40.0$. The grid spacing in each dimension is $\Delta x=\Delta y=\Delta z=0.2$, which leads to a mutual time step in MEEP and Octopus of $\Delta t=0.1$. We chose a time delay $t_{0}=10.0$ and a frequency $\omega=2.0$.

To compare our Octopus implementation to MEEP, we evaluate the electric field, the magnetic field, the electromagnetic energy density at the box point $(5,0,0)$, and the energy inside the simulation box.

Figure 7 shows the comparison of MEEP and our Maxwell propagation implementation in Octopus. In the first panel $a)$, we plot the initial current density $\vec{j}(\vec{r},t)$ in z-direction at point $(5,0,0)$. It can be seen that the maximum is according to the time delay at time $t=t_{0}$. Due to the spatial shape of the current density, the maximum value of $\vec{j}(\vec{r},t)$ is damped by the factor $e^{-5^{2}/2}$. The next two panels $b)$ and $c)$ in Fig. 7 show the electric field in z-direction respectively the magnetic field in y-direction also both at the point $(5,0,0)$. Both field signals are shifted by $\Delta t=5.0$ to the future compared to the maximum current density signal. This follows from the fact that the current maximum which causes the electromagnetic field reaction arises along the z-axis for $x=y=0$ at time $t=t_{0}$. Hence, the center of the electromagnetic reaction signal arrives at $t=t_{0}+\Delta t=15.0$. We see that both simulations lead to the same electric and magnetic field behavior in time at the selected grid point. Furthermore, the Maxwell energy densities at this point calculated by Octopus and MEEP also match, and are plotted in Fig. 7 $d)$. The last two graphs in panel $e)$ show the total Maxwell energy inside the box with $-40.0\leq x\leq 40.0$, $-40.0\leq y\leq 40.0$, and $-40.0\leq z\leq 40.0$. Also here we observe very good agreement with MEEP, which together with the previous cases shows the correctness of our Riemann-Silberstein implementation of Maxwell’s equations.

VI Applications

In this section, we demonstrate the significance of taking the fully self-consistent coupling of the time-dependent Kohn-Sham equations for the electrons, Ehrenfest equations for the nuclei, and Maxwell’s equations for the electromagnetic fields into account. We use our EMPKS implementation in the Octopus code that we introduced in the previous sections to study several different coupling scenarios. These range from conventional forward light-matter coupling in dipole approximation with clamped nuclei to a theory level with forward-backward self-consistent light-matter coupling including electric quadrupole and magnetic dipole terms where we also include the motion of the ions and classical Lorentz forces on the ions. An overview of the various EMPKS theory levels that we use in the present work can be found in Tab. 1. The possibility to switch on and off different degrees of freedom and their mutual coupling has the benefit that we can better isolate the impact and significance of physical mechanisms in various applications.

VI.1 System and simulation parameters

VI.1.1 Na297-dimer geometry

As a first test system for our EMPKS implementation, we have selected a nanoplasmonic system which was already investigated previously in a work by Varas et. al. Varas et al. (2015). This system consists of two almost spherical sodium nanoparticles with 297 sodium atoms each, which are arranged in a dimer configuration. In Fig. 8, we illustrate the geometry of the sodium dimer for two different distances between the sodium nanoparticles. We use the default Troullier-Martins pseudo-potentials in Octopus, which treat one valence electron per sodium atom. Consequently, the total system comprises 594 sodium atoms and 594 valence electrons Noya et al. (2007); et. al. (2018).

We have performed standard geometry optimizations for the ground state of the system. The icosahedral polyhedron is the most stable geometry for a single sodium nanoparticle. The quite large polyhedron is approximately a sphere with an effective diameter of $2R$. From the optimization we find an effective radius of $R\approx 2.61$nm. The key parameters for the composition of the dimer are the distance between the two nanoparticles and their relative orientation. We use the distance $b$ of the two centered sodium atoms of each icosahedron so that $d$ is defined as $d=b-2R$ and does not depend on the relative orientation of the two nanoparticles to each other. The two icosahedrons can be oriented in several constellations. We have used a relative orientation such that the 3-atom edges of the hexagons on the surface of the nanoparticles are lying face to face. This so called E2E configuration Varas et al. (2015) is illustrated in Fig. 8. The dimer axis is oriented parallel to the z-axis, therefore the dimer is symmetric with respect to x- and y-axis.

To investigate the effect of the total internal dipole of the system on the coupled time-evolution, we consider in the following two different distances of the nanoparticles with $d_{1}=0.1$ nm and $d_{2}=0.5$ nm which result in different dipoles. By computing the optical absorption cross section of the dimer, it was shown in Ref. Varas et al. (2015) that the system with distance $d_{1}$, shows a prominent quadrupole (Q mode) localized surface-plasmon resonance, whereas the second system with distance $d_{2}$ has a strong dipole (D mode) resonance. In the following examples, we excite the nanoplasmonic dimer with an incoming laser pulse, where we select the corresponding resonance frequencies of the Q and D mode as carrier frequency of the laser.

VI.1.2 Simulation boxes and grid alignment

For the real-time simulations of the nanoplasmonic dimer, we use in the following a Maxwell-Kohn-Sham simulation box which corresponds to grid type f) in Fig. 5. The matter Kohn-Sham grid is smaller than the Maxwell grid, but both grids have the same grid spacings in each spatial direction and all Kohn-Sham grid points lie on top of Maxwell grid points. The Kohn-Sham grid geometry is based on the so called minimum box construction Castro et al. (2006). The minimum box of a molecule consists of the union of all Cartesian grid points which lie inside a fixed radius around each ion of the system. For all simulations, we select a radius of $R_{\mathrm{min}}=0.794$ nm ($15$ a.u.). Taking the corresponding geometries for the dimer distances $d_{1}$ and $d_{2}$ into account, we obtain maximal extensions $L_{\mathrm{KS},x}$, $L_{\mathrm{KS},y}$, $L_{\mathrm{KS},z}$ in each direction given in Tab. 2. The matter grid is surrounded by a significantly larger parallelepiped shaped box for the Maxwell grid points with the extensions $L_{\mathrm{Mx},x}$, $L_{\mathrm{Mx},y}$, $L_{\mathrm{Mx},z}$ in negative and positive direction which is illustrated in Fig. 4. As grid spacing for both grids we select $0.053$ nm ($1.0$ a.u.) and a stencil order of 4 for the finite-difference discretization.

For the Kohn-Sham orbitals we use a zero Dirichlet boundary condition, whereas for the Maxwell grid we employ the combined incident plane-wave plus absorbing boundaries via PML as introduced in section V.3.5. Hence, the Maxwell simulation grid is separated into two areas, one outer for the incident plane-wave boundaries and one inner for the PML. The incident plane-wave boundary width depends on the derivative order for the finite-difference stencil times the grid spacing. In the present case, the width of the plane-wave boundary region is $0.212$ nm ($4.0$ a.u.), and we use in addition $0.265$ nm ($5.0$ a.u.) as PML region. The total inner simulation box for the Maxwell propagation is therefore limited by $a_{\mathrm{Mx},x}$, $a_{\mathrm{Mx},y}$, and $a_{\mathrm{Mx},z}$ also given in Tab. 2.

VI.1.3 Measurement regions

To distinguish near-field and far-field effects, we define different measurement regions in the simulation box where we record the electromagnetic fields, and the current density as function of time during the laser pulse propagation. Besides the mid point (mp) $\vec{r}_{\mathrm{mp}}=(0,0,0)$ that was already considered in the work of Varas et. al., we also consider two far-field points. The point (ffpx) $\vec{r}_{\mathrm{ffpx}}=(1.957~{}$nm$,0,0)=(37.0~{}$a.u.$,0,0)$ is shifted along the laser propagation axis, and (ffpy) $\vec{r}_{\mathrm{ffpy}}=(0,1.957~{}$nm$,0)=(0,37.0~{}$a.u.$,0)$ along the y-axis. Furthermore, we define a detector surface given by the parametrization

[TABLE]

The detector surface includes the far-field point $\vec{r}_{\mathrm{ffpx}}$ and the extension is determined by the box limits. We chose the limits such that all points have sufficient distance to the absorbing PML region. For ease of comparison, we compute the average of the electric field over the detector surface.

VI.1.4 Laser pulse shape

For the incident laser pulse that approaches our nanoplasmonic dimer, we use the following analytical expression for the external electric field

[TABLE]

where $\theta$ denotes the usual Heavyside step function. The laser pulse propagates with wavevector $\vec{k}=(k_{x},0,0)$ along the x-axis. The electric-field polarization is oriented along the z-axis and consequently the magnetic field oscillates parallel to the y-axis. The corresponding magnetic field $\vec{B}_{\mathrm{pw}}$ and Riemann-Silberstein vectors $\vec{F}_{\pm,\mathrm{pw}}$ are with $\vec{k}\cdot\vec{E}_{\mathrm{pw}}=\vec{k}\cdot\vec{B}_{\mathrm{pw}}=0$ and Eq. (76) given by

[TABLE]

[TABLE]

As discussed in section V.3.4, this analytical form of the Riemann-Silberstein vector $\vec{F}_{\pm,\mathrm{pw}}$ is used to update the incident plane wave boundaries for each propagation step. While the incoming laser pulse is prescribed analytically in the boundary region, we propagate the electromagnetic fields fully numerically in the interior of the simulation box.

VI.1.5 Propagators

Inside the Maxwell propagation region with ($-l_{\mathrm{Mx},x}\leq x\leq l_{\mathrm{Mx},x}$), ($-l_{\mathrm{Mx},y}\leq y\leq l_{\mathrm{Mx},y}$), ($-l_{\mathrm{Mx},z}\leq z\leq l_{\mathrm{Mx},z}$), we propagate the Kohn-Sham system with the matter ETRS propagator from Eq. (183) using the Power-Zienau-Woolley transformed MPKS Hamiltonian from Eq. (190) with multipole expansion. The Maxwell system is evolved in time by the Maxwell ETRS propagator from Eq. (127) with corresponding Hamiltonian kernel $\hat{\mathcal{H}}$ from Eq. (94), where we use only $\hat{\mathcal{H}}$ and switch off the internal diamagnetic current kernel $\hat{\mathcal{K}}$, and the current density term $\mathcal{J}$ given in Eq. (100) only contains the paramagnetic current contribution. We propagate the Riemann-Silberstein vector corresponding to the total vector potential $A^{k}_{\rm{tot}}=a^{k}+A^{k}$ of external and internal fields. The nuclei are treated classically and are propagated with Ehrenfest equations of motion Xavier Andrade (2009). Since we have the electromagnetic fields available in the simulation box, we also include the classical Lorentz force that acts on the ions. For the transversal Kohn-Sham field we use the mean-field approximation $a^{k}_{\rm{KS}}\approx a^{k}+A^{k}$ as well as the physical mass of the particles to take into account the bare vacuum fluctuations of the photon field. For the longitudinal Kohn-Sham field we use $a^{0}_{\rm{KS}}\approx a^{0}+A^{0}+a^{\rm{LDA}}_{\rm{xc}}$, where $a^{\rm{LDA}}_{\rm{xc}}$ is the adiabatic LDA exchange-correlation approximation.

In addition to the fully coupled EMPKS simulation, we propagate in addition the unperturbed Maxwell system inside the inner simulation box to get the required values for the incident plane wave plus PML boundaries according to section V.3.5. Hence, the Maxwell Hamiltonian has to be updated inside the PML boundaries by the additional PML kernel $\hat{\mathcal{G}}$ in Eq. (163).

VI.2 Results from EMPKS simulations

In the previous section, we have introduced the nanoplasmonic dimer, the incident laser pulse, and the parameters for the EMPKS time-stepping. In the following, we now focus on the actual EMPKS simulations for the dimer. In all cases, we compare our results from the different theory levels of self-consistently coupled light-matter propagations with the conventional forward-coupling only case with dipole approximation and highlight the differences.

To give an overview of the dynamics of the nanoplasmonic dimer, we provide in Ref. Jestädt et al. (2018a) and Ref. Jestädt et al. (2018b) movies which show the interaction of the dimer with the incoming laser pulse. In Fig. 9 and 10, we show two representative snapshots at different time steps for the dimer with distance $d_{1}=0.1$ nm. The movie frames are divided into four panels. Each panel contains a contour plot which shows the x-y plane of the 3D data. The upper left panel shows the absolute value of the induced paramagentic current. In the upper right panel, we show the electron localization function (ELF) to visualize the electron motion. In the lower left panel, we plot the electric field enhancement. We define this enhancement by subtracting the longitudinal electric field that corresponds to the ground state Kohn-Sham potential from the full induced electric field at a given time step. The same difference is used in the lower right panel, which shows the Maxwell energy difference of the induced field and the ground state field energy. Above these four panels, we also show a contour plot of the laser pulse propagation which includes at the center of the image also the geometry of the nanoplasmonic dimer.

VI.2.1 Electric field enhancement

In the work of Varas et. al. a large field enhancement right at the mid point between the two sodium nanoparticles was found. This enhancement originates solely from longitudinal photon degrees of freedom since in their work only a conventional time-dependent Kohn-Sham calculation without coupling to Maxwell’s equations was performed. The longitudinal component of the electric field can be obtained in this case from the scalar Kohn-Sham potential $a_{\mathrm{KS}}^{0}=a^{0}+a^{0}_{\rm{xc}}+A^{0}$ according to

[TABLE]

which is just the external and the Hartree potential. On the other hand, we can expect that also field enhancements occur for the transverse components of the electromagnetic field. To investigate this, we use our EMPKS approach in the following to study the field enhancements in the fully coupled case.

In Fig. 11, we show in panel a) the electric field amplitude of the incident laser, in panel b) the current density at the mid point $\vec{r}_{\mathrm{mp}}$, and in panels c)-e) the electric field at the mid point $\vec{r}_{\mathrm{mp}}$, the electric field at the far field point $\vec{r}_{\mathrm{ffpx}}$, and the average of the electric field over the detector surface, respectively.

All blue curves in Fig. 11 refer to the F@ED theory level. These results are fully identical with the previous results of Varas et. al., but have here been obtained from our EMPKS implementation by switching off the matter to Maxwell back-reaction. The exact overlap with the data of Varas et. al. provides a further consistency-check of our implementation.

The black curve in panel a) and the light gray curve in panel c) display the same initial cosinusoidal shaped free laser pulse that passes through the simulation box without any matter interaction. We include the pulse from panel a) again in panel c) to facilitate the comparison of the incident field amplitude (gray) in comparison to the actual electric field values when light-matter coupling is taking place (blue and green). As was noted by Varas et al. for the F@ED theory level, a field enhancement of approximately a factor of three compared to the incident laser amplitude can be found at the mid point as can be seen from the blue curve in panel c). Varas et. al. also noticed a delay of the total induced electric field maximum compared to the maximum field amplitude of the driving laser. This shift of maxima can also be seen in panel c) by comparing the maxima of the gray and blue curves.

We now switch on the back-reaction of the matter on the electromagnetic fields. We denote this with theory level FB@ED (cf. Tab. 1) and use green curves in the figure. Including the back-reaction, we find in panel c) a similar enhancement for very short times, a slightly increased enhancement for intermediate times and smaller beating for longer times. Overall, adding the back-reaction is increasing the field enhancement at the mid point while the external laser is active. Similar to the F@ED case of Varas et. al., a similar delay of the field maximum compared to the driving laser also appears in the forward-backward coupling case. The situation is quite different in the far field which is illustrated in panel d). Here the fully self-consistent light matter coupling (green curve) shows a field enhancement which is almost twice as large as in the forward-only coupling case. While the selected far field point is lying on the laser propagation axis and could be considered special, this picture for the field enhancement in the far field is confirmed if we compare the field enhancements averaged over the detector surface as shown in panel e). It is also worth noting, that the electric field is mostly in phase with the incident laser in the far field, for both F@ED and FB@ED. On the other hand a small phase and frequency shift occurs in the fully coupled case in the near-field in panel c).

By construction, we know that the total electric field for F@ED coupling and consequently all related field enhancements in this case are only longitudinal. To analyze further the nature of the field enhancement when the coupling to Maxwell’s equations is enabled, we have performed a Helmholtz-decomposition of the electric field in the fully-coupled FB@ED case. We find that also for FB@ED coupling the main contribution to the field enhancement arises from the longitudinal component. This implies that the additional back-reaction is causing the longitudinal components to leak much farther into the far-field. There is also a small enhancement of the transverse field, but but this enhancement is an order of magnitude smaller compared to the longitudinal fields. Since longitudinal fields are always bound to charges in space and since we do not consider electronic loss, the longitudinal field enhancement that we measure at the boundary surface is not leaving the system. However, since we also have a transversal field enhancement there is now also a contribution that is radiating and that is reaching our absorbing boundary PML region.

Next, we turn our attention to the case of a larger nanoparticle separation of $d_{2}=0.5$ nm with a larger internal dipole. In this case, we drive the plasmonic dimer with the frequency $\omega_{2}=2.83$ eV $(0.104$ a.u.$)$ that corresponds to the dipole localized surface-plasmon resonance (D mode). We show in Fig. 12 the results for this case and use the same ordering and labels for the panels as before in Fig. 11. Due to the larger distance between the two sodium nanoparticles, the absolute value of the current density is significantly smaller compared to the previous case with the smaller distance. In panel b) of Fig. 12 it can be seen that the induced current densities for the forward coupling F@ED case are larger than for the forward and backward coupled FB@ED case. Looking at the electric field enhancements at the mid point shown in panel c), we find a significantly smaller field enhancement in the FB@ED coupled case. This is in contrast to the smaller dimer separation in Fig. 11, where at least for short times the FB@ED coupling has induced larger field enhancements compared to the only forward coupled case. As before in the case of the $d_{1}$ separation, we find here that the average of the electric field over the detector surface in the far field as shown in panel e) is mostly locked to the phase of the incident laser. In the near field at the mid point between the two nanoparticles the phase and frequency shifts are again larger for F@ED coupling at short times and the phase turns even to the opposite sign compared to FB@ED coupling for intermediate times. As before we have also performed here a Helmholtz-decomposition and also for the larger nanoparticle separation the longitudinal component of the electric field is leaking much more into the far field in the FB@ED coupled case than in the F@ED coupled case.

To illustrate the relative magnitude of longitudinal and transverse components of the electric field quantitatively, we show in Figures Fig. 13 and Fig. 14 the temporal evolution of both field components in z-direction at point $\vec{r}_{\mathrm{mp}}$ and $\vec{r}_{\mathrm{ffpx}}$ respectively. As reference, we also plot the total field with solid lines. Besides the fact that the longitudinal enhancement is about one order of magnitude larger for the $d_{1}$ distance and even two orders of magnitude larger for the $d_{2}$ distance, we see some phase and frequency shifts between the different fields. The phase shift between the longitudinal and the total field in Fig. 13 a) and Fig. 14 a) is very small for both distances. Although the phase shift to the transverse field is rather large, its small amplitude leads only to a minor contribution for the total field. This behavior differs at the detector surface point $\vec{r}_{\mathrm{ffpx}}$ illustrated in Fig. 13 c) and Fig. 14 c). Here, both the longitudinal and the transverse field have almost the same magnitude and show a clear phase shift. The behaviour in the far field in Fig. 13 d) and Fig. 14 d) exhibits besides a phase shift also a slight frequency shift. Consequently, the incident laser pulse interferes with the induced transverse field wich results in a frequency modification of the outgoing laser. Since the transverse field, which reaches the far-field detector region, propagates freely, our detector point measures this frequency shift.

Up to now, we have looked at the electric field enhancement as function of time. In Fig. 15, we show contour plots of the transversal electric field enhancement as function of space in the x-z plane of the nanoplasmonic dimer. The two plots in the top row correspond to the point in time where the incoming laser pulse reaches its maximum, whereas the two plots in the bottom row correspond to the point in time where the electric field enhancement reaches its maximum. The two plots in the left column have been computed with light-matter forward coupling only and in the two plots in the right column we have used self-consistent forward-backward coupling. As can be seen, the forward coupled cases show a rather uniform electric field in the plane which is due to the dipole approximation and the fact that the incident wavelength is rather large on the scale of the dimer. On the other hand for the fully coupled case on the right hand side local field effects are clearly visible. In particular in the plot on the bottom right it can be seen that at the maximum of the field enhancement the transversal field contribution in fact counter acts the longitudinal contribution since it has turned to a negative sign in most regions of space.

VI.2.2 Next order in multipole coupling and energies

So far we have focussed only on the electric dipole approximation. In the following, we include also the next order in the multipole expansion in the Kohn-Sham Hamiltonian, namely the electric quadrupole (EQ) term and the magnetic dipole (MD) term from Eqs. (LABEL:eq_electric_quadrupole_hamiltonian) and (192). As before, we consider light-matter forward coupling as well as self-consistent forward-backward coupling. This leads to the theory levels F@(ED+MD+EQ) and FB@(ED+MD+EQ) respectively. In Fig. 16 we show for the $d_{1}=0.1$ nm configuration as before in panel a) the electric field amplitude of the incoming laser. In panel b) we plot the electronic energy of the Kohn-Sham system. The energies of the four different coupling theory levels F@ED, F@(ED+MD+EQ), FB@ED, and FB@(ED+MD+EQ) split up significantly in time and stay almost constant when the laser has passed the box. In all cases the matter absorbs energy during the time when the external laser propagates through the box. Similar to the electric fields enhancements, which exhibit a clear delay of reaction to the initial laser, also the energy gain for the matter is slightly delayed with respect to the external laser pulse. The blue curve shows the conventional F@ED in dipole approximation without back reaction of the matter to the field. If we add the second order multipole coupling terms, the corresponding F@(ED+MD+EQ) run (yellow curve) gains more energy than in dipole approximation. This can be understood by the fact that we drive with the frequency of the incoming laser in this case the Q mode of the nanoplasmonic dimer which has a quadrupole nature. The EQ and MD coupling terms in the Hamiltonian can efficiently couple to this mode which results in a larger transfer of energy to the electrons.

Switching on the backward coupling reduces the energy absorption of the matter. Both, the FB@ED and FB@(ED+MD+EQ) cases remain energetically below the reference F@ED run. Again, the additional multipole terms in the FB@(ED+MD+EQ) run increase the energy curve compared to the FB@ED run. As before, this shift can be attributed to the more efficient energy transfer mediated by the additional MD+EQ coupling term.

In panel c) we show the total Maxwell energy, which corresponds to an integration of the Maxwell energy density over the whole simulation box. First, we note that all Maxwell energies oscillate with twice the frequency of the initial laser. This simply arises due to the squared electric and squared magnetic fields in the expression for the energy density. On the other hand, the peak positions depend on the respective phase shift of the electric and magnetic fields. As we already noticed for the electric field enhancement, the dominant part of the total electric field is given by the longitudinal component, so that also for the Maxwell energy the largest contribution originates from the longitudinal part of the electric field. The transverse electric and the magnetic field have only a small contribution to the Maxwell energy. This can be seen by comparing the scale of the black curve, which corresponds to the energy of the purely transversal incoming laser pulse, with the blue curve which shows the Maxwell energy for the forward coupled case in dipole approximation. By also adding here the electric quadrupole and magnetic dipole terms to the Hamiltonian, we find a substantial gain in electromagnetic energy which exceeds the scale of the energy of the incoming laser by more than a factor of five. Again most of this energy is stored in the longitudinal component of the electric field.

When we look at the two forward-backward coupled cases, FB@ED and FB@(ED+MD+EQ), we find similar to the energy transfer to the electrons also a smaller transfer of energy to the Maxwell fields. Overall, the energy oscillations are about six times weaker and clearly phase shifted compared to the forward-coupled case. Additionally, we note that in the FB@ED and FB@(ED+MD+EQ) cases extra energy is stored inside the electromagnetic fields once the incoming laser has passed the system. This is evident from the energy plateau in both cases after about 20 fs, which is above the energy value at the initial time.

Finally, we turn our attention to the second dimer configuration with $d_{2}=0.5$ nm. The matter and Maxwell energies for this case are shown in Fig. (17). At this point is is useful to recall that we drive the nanoplasmonic dimer for the $d_{2}$ distance with a frequency that is resonant to the D mode which has dipole character. While overall a similar situation emerges as for the smaller $d_{1}$ distance, a notable difference is here that the inclusion of the higher order multipole terms, MD and EQ, has little effect on the energy gains for the electrons and the electromagnetic fields. The F@ED and F@(ED+MD+EQ) energies, as well as the FB@ED and FB@(ED+MD+EQ) energies are almost on top of each other. In other words, since we couple to a dipole mode of the system, the MD and ED coupling terms have almost no effect. In a perturbative analysis this can be understood from the selection rules of the involved states for the MD and EQ coupling Hamiltonians. From this we can conclude that it depends on the symmetry of the excited modes if higher order multipole terms become important.

The common fact, that for both distances $d_{1}$ and $d_{2}$ the forward- and backward coupling matter energies remain always below the forward coupling runs demonstrates that the matter absorbs less energy if the back-reaction is taken into account. In addition to the larger absorption of energy, the forward coupling causes larger Maxwell energy amplitudes inside the simulation box. This is remarkable since we observed in Fig. 11 and Fig. 12 that the self-consistent forward-backward coupling yields a larger enhancement of the field. As consequence, in some regions of the dimer large field enhancements occur, but the mean amplification is clearly weaker than for the only forward coupled cases. Furthermore, the forward coupling runs break energy conservation, since the laser pumps the matter system without any loss. In the forward-backward coupled simulations this is not as severe anymore, and explains that the energy absorption and the mean Maxwell field enhancement is always smaller compared to the forward coupling runs. The situation would be entirely different if we could enclose the laser pulse completely in the Maxwell box. Then the pulse would not be external anymore and in the forward-backward case full energy conservation holds. For optical wavelengths this requires enormous Maxwell simulation boxes if atomic scale grid spacings are used. But it becomes feasible for hard x-rays, where much smaller Maxwell grids are needed due to the shorter wavelength.

VI.2.3 Electromagnetic detectors and harmonic generation

It is common practice in most quantum simulations to use matter expectation values to approximate optical spectra. A prime example is the dipole expectation value. The Fourier transform of the dipole is routinely used to compute absorption spectra in the linear case or high-harmonic spectra in the non-linear case. With a Maxwell propagator at hand it becomes now feasible to directly look at the emitted radiation. In other words, it is not necessary anymore to take a detour and to approximate the emitted radiation from matter observables. Rather, we can look directly at the emitted radiation. This provides a paradigm shift, since we now can essentially define ”electromagnetic detectors” in the far-field close to the box boundaries of the Maxwell simulation box which allow to perform numerical simulations that very closely resemble the experimental situation.

To showcase this new paradigm, we plot in Fig. 18 in panels b)-e) dipole expectation values for our nanoplasmonic dimer. In panel f) we show the x-component of the electric field in the far field, which corresponds to our electromagnetic detector in this case. Furthermore, in Fig. 19, we show the Fourier transforms of the last two panels from Fig. 18. Panel a) from Fig. 19 is the Fourier transform of panel e) in Fig. 18 and panel b) is the Fourier transform of panel f).

Several observations can be deduced from these results. First, comparing the forward coupling cases, F@ED and F@(ED+MD+EQ), with the fully coupled dipole signals, FB@ED and FB@(ED+MD+EQ), shows a damping in the amplitude which is consistent with the observations we already made for the electric field enhancements and the energies. Second, adding MD and EQ coupling terms to the Hamiltonian produces in both, the only forward coupled case as well as in the forward-backward coupled case dipole signals which oscillate with twice the frequency of the incoming laser. In other words, second harmonic generation is only found in this case if we go beyond the dipole approximation. And finally, looking at the electric field in the far-field we also see the second harmonic signal. This signal only emerges in the fully forward-backward coupled case, since in the forward coupled case the matter can not influence the electromagnetic fields. This is a prime example that a self-consistent forward-backward coupling to Maxwell’s equations is needed to achieve a comprehensive physical picture.

VI.2.4 Ion motion

In all the simulations that we have considered so far, we have used the Born-Oppenheimer approximation and have clamped the classical nuclei at the optimized ground state geometry. In this section we now release this constraint and we also allow the nuclei to move according to the Ehrenfest equations of motion Xavier Andrade (2009) and the classical Lorentz forces that we introduced in Sec. IV.1. For the Lorentz forces on the nuclei we can take directly the electromagnetic fields that we propagate on our Riemann-Silberstein grid. This allows to capture nuclear forces due to local field effects. For all the following cases, we take as initial condition for the Ehrenfest equations the atomic positions of the optimized ground state and set the initial velocities to zero. This effectively corresponds to a rather ”cold” nuclear subsystem. More sophisticated velocity distributions could be used, e.g. thermalized velocity distributions from molecular dynamics runs coupled to a thermostat, but we leave such temperature studies for the future.

In Fig. 20 we show again matter and Maxwell energies for the nanoparticle distance $d_{1}$, but now for moving ions. As before, panel a) shows the incoming laser, panel b) the matter energies and panel c) the Maxwell energies. In addition we now also add panel d) which shows the sum of the kinetic energy of all nuclei as function of time.

The additional ionic motion causes on this rather short time scale of about 30 fs some additional fluctuations in the matter energy evolution, but the main behavior is very similar to the case with fixed ions. However, looking at the Maxwell energies in panel c) reveals a small overall decrease of the Maxwell energy in the forward coupled and a rather strong decrease in the self-consistent forward-backward case. Since the electronic energy remains almost identical to the case of clamped ions, the losses in the Maxwell energy are directly transferred to the nuclei. This can be seen in panel d) where the kinetic energy of the nuclei grows. It is quite remarkable that this increase is rather fast (30 fs) as nuclear motion is typically attributed to take place on a pico-second time scale.

VI.2.5 Comparison of different density functionals

All of the result shown so far have been computed exclusively with TDLDA as choice for the approximate exchange-correlation functional for the longitudinal part of the light-matter interaction. To assess the relative importance of exchange-correlation effects versus self-consistent light-matter interaction, we have performed also simulations with the PBE functional. In Fig. 21 we show in panel a) the electric field at the origin, in panel b) the electromagnetic energy, and in panel c) the electric dipole in z-direction. In all panels, we compare the TDPBE results (dashed lines) with the TDLDA results (solid lines). From the results it is apparent that the difference between TDPBE and TDLDA (dashed vs. solid lines) is much smaller than the difference between only forward coupling and self-consistent forward-backward coupling (blue vs. green lines). While the change in amplitude is minor for the electric field and the dipole, a clear frequency shift for forward-backward coupling is visible already after short time. This shift has also an effect on the electromagnetic energy which shows a larger deviation also in amplitude. In other words, for the present example it is more important to switch on self-consistent forward-backward light-matter interactions than to include further exchange-correlation contributions to the effective Kohn-Sham potentials. This supports the need for a self-consistent coupling to Maxwell’s equations to achieve a comprehensive description of light-matter interactions.

VII Summary and conclusions

In the present work we provide a comprehensive derivation of a formally exact density-functional approach for fully self-consistent light-matter interactions of photons, electrons, and effective nuclei. This approach corresponds to an exact reformulation of the full quantum problem of non-relativistic QED in terms of non-linear quantum Navier-Stokes equations. We also introduce a Kohn-Sham construction, which allows to map the dynamics of interacting particles to the effective dynamics of non-interacting particles with the same densities. This turns the approach into a computationally feasible method. In the mean-field approximation for the effective nuclei, our approach corresponds to coupled Ehrenfest-Maxwell-Pauli-Kohn-Sham equations.

Using the Riemann-Silberstein formulation of classical electrodynamics, we then show how to couple the Ehrenfest equations for the nuclei and the Pauli-Kohn-Sham equations for the electrons self-consistently to the time-evolution of the electromagnetic fields. The reformulation of Maxwell’s equations in a Schrödinger-like form allows us to use time-evolution techniques that have originally been developed for the dynamics of wavefunctions in quantum mechanics. With help of the Riemann-Silberstein vector we can seamlessly combine fully-microscopic electrodynamics with macroscopic electrodynamics for linear media. This allows us to easily incorporate macroscopic optical elements like lenses or mirrors in a microscopic description and opens the path to many applications, such as, e.g., molecules in optical cavities or close to boundaries. For the coupling between light and matter, we employ in practice a multipole expansion based on the Power-Zienau-Woolley transformation. We introduce a predictor-corrector scheme for a self-consistent coupling of light and matter. To address the different time and length-scales which arise for the coupling of molecular or nano-scale systems to light with optical wavelengths, we introduce a multi-scale approach in space and time. For an efficient absorption of outgoing electromagnetic waves, we present a perfectly matched layer for the Riemann-Silberstein vector. We have implemented our approach in the real-time real-space code Octopus and illustrate the approach with an example of light-matter interactions for a nanoplasmonic sodium dimer. The results from the example show that a self-consistent forward-backward coupling of light and matter is necessary for a comprehensive physical description of the system. Including the propagation of the microscopic electrodynamical fields in the real-time simulation allows us to define electromagnetic detectors in the far-field of the simulation box. In this way, it is not necessary anymore to take a detour and to approximate the emitted radiation from matter observables. Rather, we can look directly at the emitted radiation. This is a distinguished feature of our implementation and has direct application for many spectroscopies.

For the nanoplasmonic system that we have selected here as an example, we find that the difference for observables computed with the local density approximation and gradient-corrected functionals is minor compared to the difference of only forward coupling and self-consistent forward-backward light-matter interactions. Very often in literature the discrepancies between theory and experiment are attributed to missing exchange-correlation effects. From the present results, we can at least say that self-consistent forward-backward light-matter couplings should always be included and that their contribution can easily match the magnitude of longitudinal exchange-correlation contributions or even exceed them.

While the Ehrenfest-Maxwell-Pauli-Kohn-Sham limit is the level of the theory for which we present all the results in the present work, the exact density-functional formulation of the Pauli-Fierz field theory provides a sound framework to go beyond this mean-field limit. Further steps along these lines will be the subject of future investigations.

VIII Outlook

Although we have touched in this work many aspects that enter fully self-consistent light-matter coupling, there are still several tasks that remain to be addressed in the future. As most important task, functionals for the exchange-correlation terms that we introduced in the density-functional framework for non-relativistic QED have to be constructed. Including such exchange-correlation terms will allow to go beyond the classical mean-field description for photons (and likewise for the effective nuclei). In particular, this allows to include vacuum fluctuations of the electromagnetic field in the description beyond the physical-mass approximation.

Further, we have considered in the implemenetation so far only the lowest orders of a multipole expansion based on the Power-Zienau-Woolley transformation. Clearly, full-minimal coupling is required for a complete description of non-relativistic QED. In our implementation, we already compute the full vector potential with spatial resolution on the grid. The next step is to plug this potential back into the Hamiltonian and to treat full minimal coupling with full spatial resolution. It is then valuable to have both, the full minimal-coupling, as well as the multipole expansion available in the same implementation. When certain effects appear for a given application in full minimal-coupling, it is then possible to go back to the multipole expansion and to check in which order the effects start to emerge. One such example is the second-harmonic generation that we discussed which occurred in our nanoplasmonic system only by including magnetic dipole and electric quadrupole coupling.

For the back-coupling from matter to the electromagnetic fields, we have used in the present work the standard paramagnetic current density. We still need to include the diamagnetic current and magnetization current contributions of the electrons, and also the ionic current of the nuclei for a complete description of the back-coupling.

Another line of research that is required in the future is the dependence of the pseudo-potentials on the self-consistent light matter coupling. In the present work we have used the conventional field free pseudo-potentials. However, when constructing the pseudo-potentials the proper starting point should be a light-matter coupled all-electron Hamiltonian which introduces a dependence on the vector potential in the pseudo-potential. All the computations that we have shown start from the field-free ground-state of the matter. It is clear that also the ground state should be found when light-matter coupling is already switched on. This requires to solve stationary Maxwell equations self-consistently coupled to the stationary Kohn-Sham equations. Temperature effects can be included in our framework by starting from an equilibrated thermal initial state for the ions and initiating the non-equilibrium dynamics afterwards. Finally, we have focussed here on finite quantum systems. Extending the present formulation to periodic systems opens a whole new class of possibilities, ranging from 1D and 2D to 3D periodic structures.

Having all these additional aspects taken into account, allows to perform full-featured ab-initio simulations for a vast set of physical applications. Examples include systems studied in nano-optics and nano-plasmonics, (photo) electrocatalysis, light-matter coupling in 2D materials, cases where laser pulses carry orbital angular momentum, or light-tailored chemical reactions in optical cavities to name only a few.

Acknowledgements.

We acknowledge useful discussions with Nicolas Tancogne-Dejean, Aaron Kelly, Johannes Flick, Christian Schäfer, Norah Hoffman, Alberto Castro, and Hardy Gross. We thank Michele Compostella for developing the visualization software that has been used for the present work. Angel Rubio acknowledges financial support by the European Research Council (ERC-2015-AdG-694097) and Grupos Consolidados (IT578-13). The Flatiron Institute is a division of the Simons Foundation.