Quantum Electrodynamical Density-Functional Theory: Bridging Quantum Optics and Electronic-Structure Theory

In this work we give a comprehensive derivation of an exact and numerically feasible method to perform ab-initio calculations of quantum particles interacting with a quantized electromagnetic field. We present a hierachy of density-functional-type theories that describe the interaction of charged particles with photons and introduce the appropriate Kohn-Sham schemes. We show how the evolution of a system described by quantum electrodynamics in Coulomb gauge is uniquely determined by its initial state and two reduced quantities. These two fundamental observables, the polarization of the Dirac field and the vector potential of the photon field, can be calculated by solving two coupled, non-linear evolution equations without the need to explicitly determine the (numerically infeasible) many-body wave function of the coupled quantum system. To find reliable approximations to the implicit functionals we present the according Kohn-Sham construction. In the non-relativistic limit this density-functional-type theory of quantum electrodynamics reduces to the density-functional reformulation of the Pauli-Fierz Hamiltonian, which is based on the current density of the electrons and the vector potential of the photon field. By making further approximations, e.g. restricting the allowed modes of the photon field, we derive further density functional-type theories of coupled matter-photon systems for the corresponding approximate Hamiltonians. In the limit of only two sites and one mode we deduce the according effective theory for the two-site Hubbard model coupled to one photonic mode. This model system is used to illustrate the basic ideas of a density-functional reformulation in great detail and we present the exact Kohn-Sham potentials for our coupled matter-photon model system.


I. INTRODUCTION
The behavior of elementary charged particles, like electrons and positrons is governed by quantum electrodynamics (QED). In this theory the quantum particles interact via the exchange of the quanta of light, i.e. the photons [1][2][3]. Thus in principle we have to consider the quantum nature of the charged particles as well as of the light field. However, in several important cases we can focus almost exclusively on either the charged particles or the photons, while employing crude approximations for the other degrees of freedom.
In condensed matter physics and quantum chemistry the quantum nature of light can usually be ignored and the interaction between the charged quantum particles is approximated by the instantaneous Coulomb interaction. However, even then the resulting quantum mechanical equations (usually the many-body Schrödinger equation), where the electromagnetic fields are treated classically through the solution of the Maxwell equations, are solvable only for very simple systems. This lies ulti-mately in our incapability of handling the huge number of degrees of freedom of many-particle systems and consequently in our inability to determine the many-body states. This so-called many-body problem spawned a lot of interest into the question whether one can devise a closed set of equations for reduced quantities which do not involve the explicit solution of the full quantum mechanical equations and in which the many-body correlations can be approximated efficiently. Pursuits in this direction have led to various approaches such as, among others, many-body Green's function theories [4,5], (reduced) density-matrix theories [6] and density-functional theories [7][8][9][10]. These approaches differ in the complexity of the reduced quantity, which is used to calculate the various observables of interest. Especially density-functional theories, which are based on the simplest of those (functional) variables, the one-particle density (current), have proven to be exceptionally successful [11]. Their success can be attributed to the unprecedented balance between accuracy and numerical feasibility [12], which allows at present to treat several thousands of atoms [13]. Although the different flavors of densityfunctional theories cover most of the traditional problems of physics and chemistry (including approaches that combine classical Maxwell dynamics with the quantum particles [14][15][16][17][18]), by construction these theories cannot treat problems involving the quantum nature of light.
In quantum optics, on the other hand, the focus is on the photons, while usually simple approximations for the charged particles are employed, e.g. a few-level approximation. However, even in this situation the solution of the resulting equations [19,20] is only possible in simple cases (again due to the large number of degrees of freedom) and usually simplified model Hamiltonians, e.g. the Dicke model realized in a cavity [21][22][23], are employed to describe these physical situations. Already the validity of these effective Hamiltonians and their properties can be a matter of debate [24][25][26] and often further simplifications are adopted such as the Jaynes-Cummings model in the rotating-wave approximation. The rapid progress in quantum-optical experiments on the other hand, especially in the field of cavity QED [27][28][29][30] and circuit QED [31,32] , allows to study and control multi-particle systems ultra-strongly coupled to photons [33][34][35][36], where such a simple approximative treatment is no longer valid [37]. This new regime of light-matter interaction is widely unexplored for, e.g., molecular physics and material sciences [38]. Possibilities like altering and strongly influencing the chemical reactions of a molecule in the presence of a cavity mode or setting the matter into new non-equilibrium states with novel properties, e.g. lightinduced superconductivity [39], arise. Specifically in such situations an oversimplified treatment of the charged particles may no longer be allowed and an approach that considers both, the quantum nature of the light field as well as of the charged particles is needed.
In this work we give a comprehensive derivation of an exact and numerically feasible method that generalizes ideas of time-dependent density-functional theory (TDDFT). This method bridges the gap between the above two extreme cases and provides a scheme to perform ab-initio calculations of quantum particles coupled to photons. The electron-photon generalization of TDDFT in describing non-relativistic many-electron systems coupled to photon modes of mesoscopic cavities was introduced in Ref. [40]. Here we provide a general framework describing fully coupled electron-photon systems in most possible regimes/systems ranging from effective model Hamiltonians to strongly relativistic cases, which has been introduced in Refs. [41,42]. For clarity we divide the following presentation in two parts: We first demonstrate the basic ideas in a simple model system and then show how these concepts can be used in the case of general coupled matter-photon problems. A summary of all findings of the present work for the time-dependent density-functional description of QED at different levels of approximations, namely the basic variables, initial conditions and fundamental Kohn-Sham multicomponent equations is given in appendix F.
We start considering a simple model system for charged matter coupled to photons: the two-site Hubbard model interacting with one photonic mode. By employing density-functional ideas we show how one can solve this quantum-mechanical problem without the need to explic-itly calculate the complex many-body wave function. Instead, we derive equations of motion for a pair of reduced quantities from which all physical observables can by determined. We demonstrate that these equations have unique solutions and can be used to calculate the basic reduced quantities (here the basic pair of reduced quantities is the charge density of the particle and the potential induced by the photons) of the coupled problem. Therefore we here reformulate the coupled matter-photon problem in terms of an effective theory, that we call in the following a model of quantum electrodynamical densityfunctional theory (QEDFT). Since an explicit calculation of the coupled wave function is not needed, this approach allows to determine properties of the matter-photon system in a numerically feasible way. We introduce an new Kohn-Sham scheme to approximate the unknown functionals in the basic equations of motion and present results for a simple approximation. We compare these results to the exact Kohn-Sham functionals and identify shortcomings and indicate improvements.
Based on the ideas developed in the first part of this work we repeat the steps illustrated in our example but now we construct a density-functional reformulation for the full theory of QED [41,42]. We show that a straightforward approach based on the current and the potential leads to problems and that a consistent densityfunctional reformulation of QED has to be based on the polarization and the potential which is generated by the photons. This approach to the fully coupled QED problem we denote as relativistic QEDFT, and we present the corresponding Kohn-Sham construction and give the simplest approximation to the unknown functionals. In the following we then demonstrate how relativistic QEDFT reduces in the non-relativistic limit to its non-relativistic version of the corresponding non-relativistic Hamiltonian. By employing further approximations on the matter system or on the photon field a family of different approximate QEDFTs is introduced, which are consistent with their respective approximate Hamiltonians. At this level we recover the theory of Ref. [40]. In lowest order we rederive the model QEDFT of the first part of this work. Therefore, we demonstrate how all different flavours of QEDFT are just approximations to relativistic QEDFT in the same manner as different physical Hamiltonians are merely approximations to the QED Hamiltonian. Furthermore, by ignoring all photonic degrees of freedom, we find the standard formulations of TDDFT which are extensively used in the electronicstructure community [9,10].
Outline In Sec. II we investigate the QEDFT reformulation of a simple model of one particle coupled to one mode in great detail. The developed ideas are then employed in Sec. III to derive a QEDFT reformulation of QED. In Sec. IV we show how all different QEDFT reformulations are approximations to relativistic QEDFT. We conclude and give an outlook in Sec. V.

II. MODEL OF QEDFT
In this section, we introduce the basic formulation and underlying ideas of QEDFT. By employing a model Hamiltonian, we can almost exclusively focus on the density-functional ideas that allow a reformulation of the wave-function problem in terms of simple effective quantities. We first identify the pair of external and internal variables and then show that both are connected via a bijective mapping. As a consequence, all expectation values become functionals of the initial state and the internal pair. This allows for a reformulation of the problem in terms of two coupled equations for the internal pair. Then we introduce the Kohn-Sham construction as a way to find approximations to the unknown functionals, and show first numerical results.
To describe the dynamics of particles coupled to photons we solve an evolution equation of the form for a given initial state | Ψ 0 . Here ∂ 0 = ∂/∂x 0 with x 0 = ct and the standard relativistic (covariant) notation x ≡ (ct, r) (see also appendix A for notational conventions). The corresponding hermitean Hamiltonian has the general form where the depedence of the total Hamiltonian on t indicates an explicit time-dependence. Here the (timeindependent) HamiltonianĤ M describes the kinetic energy of the particles, i.e. how they would evolve without any perturbation, andĤ EM is the energy of the photonfield. The third term describes the coupling between the (charged) particles and the photons by the charge currentĴ µ and the Maxwell-field operatorsÂ µ (where the Einstein sum convention with the Minkowski metric g µν ≡ (1, −1, −1, −1) is implied and greek letters refer to four vectors, e.g. µ ∈ {0, 1, 2, 3}, while roman letters are restricted to spatial vectors only, e.g. k ∈ {1, 2, 3}). This term is frequently called the minimal-coupling term and arises due to the requirement of a gauge-invariant coupling between the particles and the photon field. The specific form of the operatorsĴ µ andÂ µ depends on the details of the physical situation. Finally, the last term describes how the particles interact with a (in general timedependent) classical external vector potential a µ ext and how the photons couple to a (in general time-dependent) classical external current j µ ext . While we usually have no control over how the particles and photons evolve freely or interact, i.e. the first three terms of the Hamiltonian (2), we have control over the preparation of the initial state | Ψ 0 and the external fields (a µ ext , j µ ext ). Therefore, all physical wave functions, i.e. found by solving Eq. (1), can be labeled by their initial state and external pair (a µ ext , j µ ext ), However, for any but the simplest systems the (numerically exact) solution of Eq. (1) is not feasible. Even if we decouple the matter part from the photons by employing the Coulomb-approximation (i.e. describing the exchange of photons by the respective lowest-order propagator) the resulting problem is far from trivial.
A. Two-level system coupled to one mode In this subsection, we introduce a simple model of charged particles coupled to photons. We discuss the basic concepts of a density-functional-type reformulation, identify the pair of conjugate variables and then deduce the fundamental equations of motion on which we base our QEDFT reformulation.
In order to demonstrate the basic ideas of a QEDFT we employ the simplest yet non-trivial realization of one charged particle coupled to photons: a two-site Hubbard model coupled to one photonic mode. The resulting Hamiltonian (see appendix E for a detailed derivation) reads aŝ where the kinetic energy of the charged particle is given byĤ and the energy of the photon mode readŝ Here t kin is the hopping parameter between the two sites, ω is the frequency of the photonic mode and (σ x ,σ y ,σ z ) are the Pauli matrices that obey the usual fermionic anticommutation relations. The photon creation and annihilation operators (â † andâ, respectively) obey the usual bosonic commutation relations. The current operator 1 is defined byĴ where l is a characteristic length-scale of the matter-part and λ is a dimensionless coupling constant.
1 To be precise,Ĵ is proportional to the dipole-moment operator, i.e. it is connected to the zero componentĴ 0 of the general fourcurrent operatorĴµ. To highlight the analogy in structure to the general case discussed in the later sections, we give it the units of a current and denote it byĴ The operator for the conjugate potential 2 is given bŷ where L is the length of the cubic cavity. Further, the current operator couples to the external potential a ext (t) and the potential operator to the external current j ext (t). These are the two (classical) external fields that we can use to control the dynamics. If we then fix an initial state | Ψ 0 and choose an external pair (a ext , j ext ), we usually want to solve Eq. (1) with the Hamiltonian given by Eq. (3). The resulting wave function, given in a site basis | x for the charged particle and a Fock number-state basis | n for the photons depends on the initial state and the external pair (a ext , j ext ). Thus, by varying over all possible combinations of pairs (a ext , j ext ), we scan through all physically allowed wave functions starting from a given initial state. Hence, we parametrize the relevant, i.e. physical, timedependent wave functions by | Ψ 0 and (a ext , j ext ). Since the wave functions have these dependencies, also all derived expressions, e.g. the expectation values for general operatorsÔ are determined by the initial state and the external pair (a ext , j ext ).
The idea of an exact effective theory like QEDFT is now, that we identify a different set of fundamental variables, which also allow us to label the physical wave functions (and their respective observables), and that we have a closed set of equations for these new (functional) variables, which do not involve the full wave functions explicitly. Such a functional-variable change is similar to a coordinate transformation, say from Cartesian coordinates to spherical coordinates. This can only be done if every point in one coordinate system is mapped uniquely to a point in the other coordinate system. For a functionalvariable change we thus need to have a one-to-one correspondence, i.e. bijective mapping, between the set of (allowed) pairs (a ext , j ext ) and some other set of functions (while we keep the initial state fixed). To identify the simplest new functional variables one usually employs arguments based on the Legendre transformation [43]. That is why these new functional variables are often called conjugate variables. We will consider this approach in the next sections where we investigate general QEDFT, and also show how one can determine the conjugate variables of this model system from more general formulations of QEDFT. For this simple model we simply state that a possible pair of conjugate variables is (J, A). In the next subsection we show that this functional variable-transformation is indeed allowed, i.e.
The main consequence of this result is, that from only knowing these three basic quantities we can (in principle) uniquely determine the full wave function. Accordingly, every expectation value becomes a unique functional of | Ψ 0 and (J, A). Thus, instead of trying to calculate the (numerically expensive) wave function, it is enough to determine the internal pair (J, A) for a given initial state. An obvious route to then also find a closed set of equations for these new variables is via their respective equations of motion. These equations will at the same time be used to prove the existence of the above change of variables, i.e. that the wave function is a unique functional of the initial state and the internal pair (J, A).
To find appropriate equations we first apply the Heisenberg equation of motion once and find Yet these two equations are not sufficient for our purposes: we need equations that explicitly connect (a ext , j ext ) and (J, A). Therefore we have to go to the second order in time, k = ω c and = 1 µ0c 2 . Here, Eq. (4) is the discretized version of ∂ 2 t n of standard TDDFT [44,45], and Eq. (5) is the inhomogeneous Maxwell equation for one photon mode [40].

B. Foundations of the model QEDFT
In the previous subsection we have stated that (J, A) and (a ext , j ext ) are the possible conjugate pair of the model Hamiltonian (3). In this subsection we want to demonstrate that indeed this holds true and that we can perform a variable transformation from the external pair (a ext , j ext ) to the internal pair (J, A). What we need to show is, that for a fixed initial state | Ψ 0 the mapping (a ext , j ext ) is bijective, i.e. if (a ext , j ext ) = (ã ext ,j ext ) then necessarily for the according expectation values (J, A) = (J,Ã).
To do so we first note, that in the above equations of motion every expectation value is by construction a functional of (a ext , j ext ) for a fixed initial state, i.e. they are generated by a time-propagation of | Ψ 0 with a given external pair (a ext , j ext ). Suppose now, that we fix the expectation values of the internal variables (J, A), i.e. we do not regard them as functionals but rather as functional variables. Then the above Eqs. (8) and (9) become equations for the pair (a ext , j ext ) that produce the given internal pair (J, A) via propagation of the initial state | Ψ 0 , i.e.
Obviously these equations can only have a solution, if the given internal variables are consistent with the initial state, i.e.
Here we have used the definition and every internal pair (J, A) that we consider is subject to these boundary conditions. Thus, the mapping (7) is bijective, if the corresponding Eqs. (10) and (11), which connect the internal pair (J, A) with the external pair (a ext , j ext ), allow for one and only one solution pair. Let us first note that for a given pair (J, A), Eq. (11) uniquely determines the external current j by Thus, the original problem reduces to the question whether Eq. (10) determines a ext (t) uniquely. The most general approach to answer this question is via a fixedpoint procedure similar to Ref. [46]. In the case of a discretized Schrödinger equation like Eq. (3) it should also be possible to apply a rigorous approach based on the well established theory of nonlinear ordinary differential equations [45]. However, for simplicity we follow Ref. [40] and employ the standard strategy of [47] which restricts the allowed external potentials a ext to being Taylor-expandable in time, i.e.
From Eq. (4) we can find the Taylor-coefficients of J (if they exist) by where the terms nÂ (α) and n (α) are given by their respective Heisenberg equations at t = 0 and only contain Taylor coefficients of a (β) ext for β < α. Now, assume that we have two different external potentials a ext (t) =ã ext (t). This implies, since we assumed Taylor-expandability of a ext andã ext , that there is a lowest order α for which For all orders β < α (even though the individual J (β) andJ (β) might not exist) it necessarily holds that But for α we accordingly find that provided we choose the initial state such that (n (0) = 0). Consequently, J(t) =J(t) infinitesimally later for two different external potentials a ext (t) =ã ext (t). Therefore, Eq. (10) allows only one solution and the mapping (a ext , j ext ) → (A, J) is bijective. As a consequence, since every expectation value of the quantum system becomes a functional of the internal pair (J, A), in the above Eqs. (10) and (11) we can perform a change of variables and find These coupled evolution equations have unique solutions (J, A) for the above initial conditions (12) and (13). Therefore we can, instead of solving for the manybody wave function, solve these non-linear coupled evolution equations for a given initial state and external pair (a ext , j ext ), and determine the current and the potential of the combined matter-photon system from which all observables could be computed. This is an exact reformulation of the model in terms of the current and the potential of the combined system only.
C. Kohn-Sham approach to the model QEDFT In the previous subsection we have derived a QEDFT reformulation in terms of the current and the potential. While the equation that determines the potential A is merely the classical Maxwell equation, and every term is known explicitly, the equation for the current contains implict terms. Therefore, to solve these coupled equations in practice, we need to give appropriate explicit approximations for the implict terms. Approximations based on (J, A) directly would correspond to a Thomas-Fermi-type approach to the model. As known from standard density-functional theory, such approximations are in general very crude and hard to improve upon. A more practical scheme is based on the Kohn-Sham construction, where an auxiliary quantum system is used to prescribe explicit approximations. However, the numerical costs of a Kohn-Sham approach compared to a Thomas-Fermi-type approach are increased.
The details of the Kohn-Sham construction depend on the actual auxiliary quantum system one wants to employ. The only restriction of the auxiliary system is that one can control the current and the potential by some external variables. Thus, one could even add further (unphysical) external fields to make approximations of the coupled quantum system easier. However, here we only present the simplest and most natural Kohn-Sham scheme, which is to describe the coupled quantum system by an uncoupled quantum system. To this end, we assume that we can find a factorized initial state that obeys the same initial conditions as the coupled problem (12) and (13). Especially, if the initial state of the coupled system is the same as in the uncoupled problem, then this condition is trivially fulfilled. In a next step we note, that for the uncoupled system subject to the external pair (a eff , j eff ) the equations of motion become (since λ = 0) + n([a eff , j eff ]; t)a eff (t), Now, obviously if one would choose (a eff , j eff ) = (a ext , j ext ), i.e. the external pair of the coupled problem, the uncoupled quantum system will in general lead to a different internal pair. However, can we find an effective pair that reproduces the internal pair (J, A) of the coupled system? The existence of such an effective pair can be based on equations for the uncoupled system similar to Eqs. (10) and (11). Note that before we were considering the question of uniqueness, i.e. can one have two external pairs leading to the same (J, A). Thus any internal pair (J, A) was apriori associated with an external pair (a ext , j ext ). If on the other hand we are given some internal pair (J, A), say from a different (coupled) quantum system, we do not apriori know that this internal pair can be represented by propagation of an initial state with some (a eff , j eff ). Thus this problem is equivalent to the existence of a solution to for a given pair (J, A) and | Φ 0 . As before, j eff is uniquely determined by simply rearranging Eq. (26) as while the existence of an a eff that reproduces (J, A) is less clear. Again, the most general approach to answer this question can rely on a fixed-point scheme similar to [46], or on mapping the problem to a special nonlinear Schrödinger equation [45,48,49]. Importantly, in the discretized case certain subtleties arise [49][50][51][52] that have to be treated with care [45,49]. Disregarding these more subtle points, we follow a simpler approach based on the assumption of Taylor-expandability in time of J. Then one can successively construct the Taylor coefficients of the effective potential from assuming that for the initial state | Φ 0 the expectation value n(0) = 0. Further assuming that this Taylor-series converges [44,53], we have constructed a pair Now, in order to actually predict the physical pair (J, A) via the Kohn-Sham system (and thus solve Eqs. (21) and (22)) we have to connect the coupled and the auxiliary system. To do so we make the composite mapping i.e. we employ the fact that (J, A) are unique functionals of the initial state | Ψ 0 and (a ext , j ext ). The definition of the resulting Kohn-Sham potentials and currents are then found by equalizing the functional Eqs. (21) and (22) with the according equations of the uncoupled auxiliary system. This leads to (now also indicating the appropriate dependence on the initial states) [40,42] Therefore, they are functionals of the two initial states, (J, A) and (a ext , j ext ), i.e.
With these definitions the coupled problem, starting from | Ψ 0 and subject to the external pair (a ext , j ext ), can be formally solved by the solution of an uncoupled, yet non-linear problem with initial state | Φ 0 and the Kohn-Sham pair (a KS , j KS ). The resulting equations are The self-consistent solutions of the Kohn-Sham Eqs. (29) and (30) by construction obey Eqs. (27) and (28), as well as equations of motion similar to Eqs. (25) and (26). By combining these equations we see that the solutions to the Kohn-Sham equations generate the solutions to the coupled Eqs. (21) and (22). We point out, that in the equation for the photonic mode we do not need any approximate functional. We merely need to solve a classical Maxwell equation. However, in practice it might be useful, especially when calculating non-trivial photonic expectation values, that one solves an actual (uncoupled) photon problem to have a first approximation to the photonic wave function.

D. Numerical example for the model QEDFT
In this section, we show numerical examples for our model system. We use the density-functional framework introduced in the previous sections and we explicitly construct the corresponding exact Kohn-Sham potentials. To illustrate our QEDFT approach, we focus mainly on two different examples: The first example treats a setup in resonance, where regular Rabi oscillations occur. We show results in a weak-coupling limit and in a strongcoupling limit. The second example includes the photon field initially in a coherent state. For this case, we study collapses and revivals of the Rabi oscillations.
The Hamiltonian in Eq. (3) is directly connected to the famous Jaynes-Cummings-Hamiltonian and the Rabi Hamiltonian [54][55][56][57], which is heavily investigated in where we transformed to the dimensionless external po- 2 a ext → a ext and the dimensionless external current 1 n 1 eωl j ext → j ext . Further, we also transform to a dimensionless time variable I t → t. To actually perform numerical calculations, we have to choose values for the free parameters. Here, we choose typically used values from the literature: t kin /I = 0.5, ω/I = 1, λ = (0.01, 0.1) and external fields which are set to zero j ext (t) = a ext (t) = 0. This set of parameters allow for a resonance situation, with no detuning between the transition energy of the atomic levels and the frequency of the field mode.
As discussed above the basic variables (densities) are the current operatorĴ and the operator for the field po-tentialÂ. In this two-level exampleĴ reduces toσ z and A reduces to â +â † .  If the rotating-wave approximation is applied to the Rabi Hamiltonian in Eq. (31), one recovers the Jaynes-Cummings Hamiltonian. This Hamiltonian is then analytically solvable. The rotating-wave approximation is only valid in the weak-coupling limit (λ ≈ 0.01). In the strong-coupling limit (λ ≥ 0.1), however, the rotatingwave approximation breaks down. Only recently, analytic results without the rotating-wave approxmation have been published [57]. Here we emphasize that the QEDFT approach presented in this paper is exact and does not rely on the rotating-wave approximation and hence also allows to treat strong-coupling situations.
In our first example we choose as initial state for both, the coupled many-body system and the uncoupled Kohn-Sham problem meaning the electron initially populates site one and the field is in the vacuum state. Therefore, no photon is present in the field initially. In Fig. 1, we show the inversion σ x (t), the density σ z (t) and the corresponding exact Kohn-Sham potential a KS (t) for the weak-coupling case. The atomic inversion σ x (t) shows regular Rabi oscillations. Rabi oscillations are also visible in σ z (t), where we observe the typical neck-like features [58] at t ≈ 150 and at later points in time.
To determine the exact Kohn-Sham potential for this case, we follow a fixed-point construction similar to [59]. As input for the fixed-point construction, we use the exact many-body densities. In addition, we also compare to an analytic formula for the Kohn-Sham potential for a one-electron two-site Hubbard model given in [45,50]. This expression gives an explicit formula for the dependence of the Kohn-Sham potential on the density. Such an explicit formula is only known in a few cases, while the fixed-point construction is generally valid. However, both methods yield in the present case the same results. A detailed discussion of the fixed-point construction for multicomponent systems of electrons and photons will be presented in a forthcoming work [60].
We emphasize that a propagation of the uncoupled Kohn-Sham system with the exact Kohn-Sham potential a KS (t) obtained in Fig. 1 reproduces by construction the exact many-body density (σ z (t) in the present case). However, as illustrated in Sec. IIC, if a Kohn-Sham propagation is used, the numerical expenses can be drastically reduced, since the Kohn-Sham construction effectively decouples the quantum system.
In practical calculations the exact Kohn-Sham potentials are normally not available and one has to rely on approximations. In the present case, the simplest approx- and nÂ ≈ n Â = nA. Then, from Eq. (27) we find the mean-field approximation to the Kohn-Sham potential The mean-field approximation is actually identical to the Maxwell-Schödinger approach, i.e. we treat the electromagnetic field as being essentially classical. Further, for λ → 0 and for λ → ∞ the mean-field approximation becomes asymptotically exact. In Fig. 2 and Fig. 3, we compare exact densities and exact Kohn-Sham potentials to densities and potentials, which were obtained by a self-consistend mean-field propagation. Already in the weak-coupling limit, Fig. 2, quite sizable differences between exact results and mean-field results become visible: Already at t = 0 the exact Kohn-Sham potential deviates from the mean-field potential. In the case of the densities, this leads to a frequency shift, where the meanfield density oscillates slower than the exact density. In the strong-coupling limit shown in Fig. 3, effects beyond the rotating-wave approximation are visible. In the exact Kohn-Sham potential, we see a non-regular feature at t = 30, which is also not coverd by the mean-field approximation. However, the mean-field approximation already covers at least some dynamical features of the propagation.
For the second example in this section, we start with the field initially in a coherent state. For a single field mode, coherent states [61,62] can be written as follows: In this example, we use as initial state for the many-body propagation and the Kohn-Sham propagation Here, the atomic state | g is the ground state of the electronic Hamiltonian (| g = 1 . For the field state we choose |α| 2 = â †â = 4. This example is in the spirit of the calculation in panel 3 in Ref. [55]. Hence, as shown in Fig. 4, we obtain a similar time-evolution of the inversion σ x (t) as in Ref. [55]. We see the Cummings collapse of Rabi oscillations at t = 250 followed by a quiescence up to t = 500 occuring. After t = 500, we see a revival of the Rabi oscillations. We also observe, as shown in [63], that the atomic dipole operator (here the density σ z (t)) continues to change during the interval of quiescence after the inversion collapse. As before, we show in the lowest panel the corresponding exact Kohn-Sham potential obtained via fixed-point iterations.
In Fig. 5, we show a comparison of the exact Kohn-Sham potentials and densities to the mean-field propagation. Here, we see that the mean-field approximation performs rather poorly. For this case the simple ansatz in Eq. (32) is not sufficient and more sophisticated approximations to the exact Kohn-Sham potential are necessary to reach a better agreement [64,65].
In summary, we have shown in this section the exact Kohn-Sham potentials which reproduce the dynamics of the exact many-body densities. In particular the coherent state example shows that there is a clear need for better approximations to the exact Kohn-Sham potential [40] that go beyond the mean-field level and that include correlation contributions. One possibility along these lines is provided by an approach based on the optimized effective potential (OEP) method [9,10,66]. We have already implemented such an OEP approach for the present model system and the corresponding results improve quite considerably over the mean-field approximation. The details of this general OEP approach to QEDFT are beyond the scope of the present paper and will be presented in a separate publication [67].

III. RELATIVISTIC QEDFT
After having presented the basic concepts of a QEDFT reformulation of a coupled matter-photon problem in a model system, we apply the very same ideas to the full theory of QED. While no new density-functional-type ideas have to be introduced, the intricacies of QED make the actual details more involved. A first subtlety is the gauge freedom of the photon field. In this work we choose Coulomb gauge to fix the superfluous degrees of freedom. This gauge has two distinct advantages over the other gauges: it reduces the independent components of the photon field to the two transversal (physical) polarizations, and it singles out the classical Coulomb interaction between the charged particles. Since we want to connect QEDFT to derived theories like cavity QED, where usually Coulomb-gauged photons are employed, and condensed-matter theory, where Coulomb interactions play a dominant role, the Coulomb gauge is for the present purpose the natural gauge to work in. However, we emphasize that also other gauges can be used as well [8,41,46].
We first present the standard approach to identify possible conjugate variables and introduce the basic equations of motions. While in the usual non-relativistic setting this route works just fine, in the fully relativistic situation the internal structure of the "Dirac particles", i.e. the electronic and positronic degrees of freedom, give rise to certain subtleties when performing a densityfunctionalization. Therefore, instead of using the current, we employ the polarization as a basic fundamental variable in relativistic QEDFT.

A. Equations of Quantum Electrodynamics
In the following we define the basic quantities of QED in Coulomb gauge and derive the equations of motion for the fundamental (functional) variables of the theory. We employ SI units throughout, since in the next section we perform the non-relativistic limit which is most easily done if we keep the physical constants explicit. A detailed discussion of quantizing QED in Coulomb gauge is given in appendix B.
The full QED Hamiltonian in Coulomb gauge (indicating explicit time-dependence of the Hamiltonians by t) is given bŷ HereĤ is the normal ordered (::) free Dirac Hamiltonian in the Schrödinger picture, whereψ andψ denote the Diracfield operators and γ k the Dirac matrices (see appendix B for definitions). The energy of the free photon field is given byĤ whereˆ E andˆ B are the (vector-valued) electric and magnetic field operators defined as in appendix B in terms of the Maxwell-field operatorsÂ k . We note that due to the Coulomb-gauge condition ∇ · A = 0 only the spatial components of the Maxwell field are quantized. The time component A 0 is given by the classical Coulomb field of the total charge density, which is the sum of the charge density of the Dirac field and the classical external current, and gives rise to the Coulomb term HereĴ 0 is the zero component of the Dirac current and j 0 ext is the zero component of a given external current j µ ext . In the Coulomb term the energy due to the Coulomb interaction of the external current with itself is elided. Since this term is purely multiplicative, i.e. it is equivalent to the identity operator times some real number, it does not influence the dynamics of the system and can be discarded. The rest of the coupling to the external fields is given bŷ Finally, the coupling between the quantized fields in Coulomb gauge reads aŝ Comparing to the Lorentz-gauge QED Hamiltonian [42] the main difference lies in the Coulomb term, that treats the zero component of the photon field explicitly. Without further refinements the above QED Hamiltonain is not well-defined, since it gives rise to infinities [1][2][3]. These infinities can be attributed, with the help of perturbation theory, to three divergent types of Feynman diagrams: the self-energy of the Fermions, the selfenergy of the photons (also called vacuum polarization) and the vertex corrections. These divergences vanish if we regularize the theory, e.g., by introducing frequency cut-offs in the plane-wave expansions of the fermionic as well as the bosonic field operators or by dimensional regularization [1]. Such procedures make the above Hamiltonian self-adjoint [68], but we have introduced a dependence on parameters that change the theory at smallest and largest length scales. Perturbatively one can remove these dependencies by renormalizing the theory, i.e. we first identify and then subtract the part of each of these three terms that diverges due to these parameters. The resulting three divergent counter-terms 3 can be recast as a renormalization of the mass and the field-operators of the Fermions (due to the self-energy), as a renormalization of the photonic field-operators (due to vacuum polarization) and a renormalization of the charge (due to the vertex corrections). We can do this to any order in perturbation theory due to the Ward-Takahashi identities [1]. Thus, the above QED Hamiltonian is at least perturbatively renormalizable. For simplicity, we assume in the following that one can fully renormalize the QED Hamiltonian (as has been shown for certain limits [69]) and interpret it as a bare Hamiltonian, i.e. we use the renormalized quantities 4 . That a full renormalization is possible has been shown, e.g., for the Nelson model of QED [70,71], where the divergent self-energy term shifts the spectrum of the Hamiltonian to infinity. Thus, subtracting this infinite shift, i.e. introducing a counterterm, makes the Hamiltonian well-defined (when removing the cut-offs), provided the energy of the system is below the pair-creation limit. The same condition, i.e. a stable vacuum, we need to impose also on our QED considerations as discussed in [8,41,42].
In a next step we identify the possible conjugate (functional) variables of the above QED Hamiltonian.
Here the physical, time-dependent wave function | Ψ(t) depends on the initial state and the external pair (a ext µ , j ext µ ), which is indicated by Thus, with ≡ T 0 dt d 3 r, the (negative) QED action [41,42] becomes a functional of these variables (T corresponds to an arbitrary time). Here we employed the definition of the QED Lagrangian of Eq. (B1) and defined the internal QED action with help of Eq. (B6) by Eq. (40) looks like a Legendre transformation between J µ ↔ a µ ext and A µ ↔ j µ ext . Since a Legendre transformation amounts to a change of variables, this indicates (for a fixed initial state) the possibility of transforming from (a µ ext , j µ ext ) to the conjugate variables (J µ , A µ ) 5 . If these variables would indeed be connected via a standard Legendre transformation the functional derivative with respect to a µ ext and j µ ext should give the respective conjugate variables. However, following derivations similar to [53] we find the appearence of extra terms, i.e. δÃ δa µ ext (x) These non-trivial boundary terms are due to the fact, that variations of the external fields give rise to nonzero variations of the wave function at the (arbitrary) upper boundary T (in contrast to direct variations of the wave function that are supposed to obey | δΨ(T ) = 0) [43]. These boundary terms are necessary to guarantee the causality of J µ and A µ [53]. Thus, Eqs. (41) and (42) show that a straightforward approach to demonstrate a one-to-one correspondence between (a µ ext , j µ ext ) and (J µ , A µ ) based on a Legendre transformation becomes difficult [42]. Nevertheless, usually this Legendretransformation arguments work well to identify the possible conjugate variables.
However, in the relativistic situation a further problem arises: the current has an internal structure due to the electronic and positronic degrees of freedom. The current J µ describes the net-charge flow of the negatively charged electrons and the positively charged positrons [2]. Therefore, the current expectation value can not differ between the situation of, e.g., the movement of two electrons and one positron or three electrons and two positrons. This fact, which is absent in the non-relativistic situation, will lead to problems when employing the ideas developed in subsection II B.
For the moment, however, we follow the above identification scheme and derive the basic equations of motion forĴ µ andÂ µ . Since d 3 r [Ĵ µ ( r),Ĵ 0 ( r )]f ( r ) = 0, where f ( r ) is any testfunction, the termĤ C commutes withĴ µ and the equation of motion for the four current is the same as in Lorentz gauge [42] i∂ 0Ĵ where the zero component is given by i∂ 0Ĵ 0 = −i ∇ ·ˆ J, i.e. the current obeys the conservation of charge. A different equation that determines the charge current J µ is found by the Gordon-decomposition [8], which is the evolution equation of the polarization where klj is the Levi-Cevita symbol and With the definition of bigger and smaller components of the Dirac-field operatorsψ † ( r) = φ † ( r),χ † ( r) we find that the current and the polarization are the real and imaginary part of the same operator The change of gauge only affects the equation for the photon-field operator which becomes and accordingly This is indeed the quantized Maxwell equation in Coulomb gauge.

B. Foundations of relativistic QEDFT
In this subsection we first reexamine the previous approach to relativistic QEDFT [41,42] and identify its shortcomings. We then show why physically the polarization is better suited as fundamental variable of the matter part and reformulate QED in terms of (P µ , A µ ). Already here we point out that both, a relativistic QEDFT based on the current or on the polarization, lead to the same density-functional-type theory in the non-relativistic limit.
A first restriction we impose is to fix a specific gauge for the external fields a µ ext . Since by construction external fields that only differ by a gauge transformation, i.e.ã µ ext = a µ ext + ∂ µ Λ, lead to the same current density (and polarization) 6 , the desired one-to-one corre-spondence can only hold modulo these transformations. Thus in principle we consider a bijective mapping between equivalence classes, and by fixing a gauge we take a unique respresentative of each class. For simplicity we impose a gauge condition similar to [72] a 0 ext (x) = 0.
In the following, any other gauge that keeps the initial state unchanged, i.e. the gauge function has to obey Λ(0, r) = 0, is also allowed [72]. This condition is necessary for our further investigations, since we will employ that the initial state is fixed (and thus the expectation values at t = 0), in accordance to the derivations of subsection II B. Further, we assume that the external current obeys the continuity equation ∂ 0 j 0 ext = − ∇ · j ext . This leaves the choice of the initial charge configuration j 0 ext (0, r). Since the photons only couple to moving charges this choice will not influence the dynamics of the system. Therefore, also in the case of the external currents we have an equivalence class (of possible zero components), and will therefore restrict to only prescribing j ext .
In [41,42] the one-to-one correspondence was based on the corresponding Ehrenfest equations wherê q k kin ( r) = − ecψ( r) γ k γ 0 γ · ∇ + γ · ← ∇ γ 0 γ k ψ ( r) and the D'Alembert operator reads as = ∂ 2 0 + ∂ k ∂ k . We can then reexpresŝ and thereforê If we then want to show a possible one-to-one correspondence we can follow the reasoning of Sec. II B and con-sider the uniqueness of solutions of the functional equations for given J k and A k 7 . As before in Sec. II B we can construct the external current uniquely. By defining the vector field we find from the Helmholtz decomposition of ζ = − ∇ξ + ∇ × Ξ and j ext = − ∇υ + ∇ × Υ that (52) Thus, we need to show that for given (J k , A k ) there can only be a unique a k ext 8 . To show this we first define formally construct the respective Taylor coefficients and consider two external potentials a k ext =ã k ext that differ at lowest order α. Accordingly we find in this order where 7 Note that in correspondence to the freedom of the external variable a k ext the freedom of the internal variable J k is also restricted, since J 0 is fixed by the initial state and the continuity equation for all times. Similarly the freedom of the external current j k ext is in correspondence to the freedom of the internal field A k . 8 We note at this point, that one can also fix the zero components j 0 ext and A 0 , respectively. Since A 0 is given by Eq. (B3) its first order time-derivative is, due to the continuity equation, given in terms of jext and J. Thus, to determine A 0 we only have to choose an initial condition. However, due to Eq. (B3) this choice also fixes automatically the initial condition for j 0 ext for the continuity equation ∂ 0 j 0 ext = − ∇ · jext.
While before we could conclude that the difference between the currents is necessarily non-zero provided n (0) = 0, here we find that this is not sufficient. Actually, we need to restrict the allowed potentials a ext to those that are perpendicular to n (0) . If we do this, then Eq. (54) makes the currents necessarily different and we can conclude that we have a one-to-one correspondence. This aspect was not taken into account in previous work [41,42], which is restricting effectively the one-to-one mapping to a smaller set of potentials and currents in these proofs. Still, it seems possible to find a different way to show the bijectivity of the complete mapping (J k , A k ) ↔ (a k ext , j k ext ). However, the true drawback of a relativistic QEDFT based on the current is found if we try to reproduce a given pair (J k , A k ). If we choose a current that obeys then the resulting equation that defines the Taylorcoefficient of the external potential reads by employing Eq. (53) and following the same strategy as in subsection II C This equation does not have a solution and therefore any current that obeys the above form cannot be reproduced by the respective quantum system. This does also call into doubt the possibility of exactly predicting the current of a coupled system by an uncoupled one, i.e. the Kohn-Sham construction of [41,42]. Of course we can remedy this problem by adding terms to the QED Hamiltonian that break the minimal-coupling prescription of the Lagrangian. Such procedures could then be alternatively used to provide a Kohn-Sham scheme to describe the fully coupled QED problem. The advantage of such an approach is, that still the equation for the vector potential is known explicitly in terms of the internal pair (J k , A k ). This is not the case, when we use a different basic variable for the matter part of the QED system, as we will do in the following.
To avoid the problems with the relativistic current, we will in the following base our considerations on the polarization P k . While the current describes the flow of the charge of the system (which is conserved), the polarization depends on the actual number of particles and anti-particles (which is not conserved). Therefore, the polarization can actually differ between a local current produced by, e.g., two electrons and one positron or three electrons and two positrons, in contrast to the current. To show now that for a fixed initial state | Ψ 0 we actually have (a k ext , j k ext ) we demonstrate that for a given internal pair (P k , A k ) the two coupled equations allow only for a unique solution (a k ext , j k ext ). Here we used the definitionŝ These coupled equations can only have a solution if the pair (P k , A k ) obeys the initial condition enforced by the fixed initial state | Ψ 0 , i.e.
Since the current J k is now a functional of (a k ext , j k ext ) the previous explicit construction of j k ext is no longer valid. However, if we assume (a k ext , j k ext ) both to be Taylorexpandable we find for the lowest order α on the one hand that provided P (0) 0 ( r) = Ψ 0 |P 0 ( r)|Ψ 0 = 0, which corresponds to the (local) total number of particles and antiparticles. On the other hand we have Otherwise there would exist a current j ext ( r) = 0 that fulfills However, due to the definition of the Coulomb potential as the Green's function of the Laplacian (see Eq. (B3) and (B4)) we find from the divergence of Eq. (62) that ∇ · j ext = 0 and thus the only possible current that fulfills Eq. (62) is j ext = 0. Thus, the mapping (55) is bijective (at least for Taylor-expandable external pairs (a k ext , j k ext )). Therefore we can, instead of solving the fully coupled QED problem for the (numerically infeasible) wave function | Ψ(t) , determine the exact internal pair (P k , A k ) from the coupled non-linear equations for the initial conditions (58) and (59). In order to solve these equations simultaneously we need to find approximations for the unknown functionals. The only drawback in this more general approach than the ones used in [41,42], is that now we also have an unknown functional in the classical Maxwell equation, i.e. J[P k , A k ].

C. Kohn-Sham approach to relativistic QEDFT
In this subsection we provide the adopted Kohn-Sham construction based on the internal pair (P k , A k ) and give the simplest approximation for the Kohn-Sham potential and current. As in Sec. II C, we choose our auxiliary Kohn-Sham system to be an uncoupled system. While different Kohn-Sham constructions are possible, this approach is the numerically least demanding.
In a first step, in accordance to Sec. II C, we first construct an uncoupled system that can reproduce a given internal pair (P k , A k ) of the fully coupled QED system. To do so, we first need an initial state | Φ 0 that fulfills the initial condition (58) and (59) of the full QED system. This allows that the coupled equations can only have a unique solution. Obviously, for the case of the uncoupled problem we can use a construction similar to Eqs. (52) to determine the unique j k eff . To show the existence of a solution to Eq. (65) we perform the standard Taylor-expansion construction and assume that the series converges [42,44,72]. A more general approach would be to follow a fixed-point procedure [46]. The respective Taylor-coefficients of the effective potential are given by This construction makes plausible that there exists an uncoupled system subject to the effective external fields (a k eff , j k eff ) that reproduces a given pair of a fully coupled QED problem. The above construction actually resembles the mapping for a given pair (P k , A k ). Now, to predict the internal pair (P k , A k ) of the full QED problem we again introduce a composite mapping The resulting Kohn-Sham potential and Kohn-Sham current are then given by the functional equations This allows to solve an uncoupled system instead of the fully coupled QED problem. However, as also pointed out in [42], we can only fully decouple the matter from the photon part if also the initial state is of product form, i.e. | Φ 0 = | M 0 ⊗ | EM 0 . And if we further assume that | M 0 is given in terms of a Slater-determinant we can actually map the whole problem to solving a Dirac equation with the above Kohn-Sham potential a k KS and simultaneously a classical Maxwell equation with j k KS . The mean-field approximation recovers the approximation introduced in [42] and reads as Since for simplicity we used a gauge where a 0 ext = 0 while for the photon field we employed Coulomb gauge, we have to perform an according gauge transformation to have the mean field a µ MF in either the one or the other gauge completely. This approximation is similar to the Maxwell-Schrödinger approach, that assumes the photon field to behave essentially classically.

IV. NON-RELATIVISTIC QEDFT
While for the sake of generality we have been considering the full QED problem in the previous section, we are actually mainly interested in the behaviour of condensedmatter systems or atoms and molecules that interact with photons. In such situations the external fields are usually small compared to the Schwinger-limit, i.e. we do not have pair-production is such situations. Further, we want to investigate systems, where the quantum nature of the photons becomes important. Most prominently this happens for the case of a cavity, where different boundary conditions for the Maxwell field have to be considered. These quantum-optical situations also naturally restrict the available photonic modes. Such physical situations are then well described by models of non-relativistic particles interacting with a quantized electromagnetic field, such as the Pauli-Fierz Hamiltonian (see e.g. [73,74]) or the Nelson model [70,71]. In the lowest order of approximations we find the situation of a two-level system interacting with one photonic mode, similar to the one presented in Sec. II. This simplest of models is the prime example of a quantum-optical problem.
We realize at this point, that all the conditions we had to impose in order to make our starting QED Hamiltonian well-defined, are naturally met in the situations we aim at investigating. Actually, we even do not need to adopt a field-theoretical treatment for the particles in the first place and usually only need to take into account a few photonic modes. Such an approach would avoid a lot of unpleasent problems in connection with renormalization and regularization of these theories. However, one would then need to introduce a new QEDFT approach for every new type of model Hamiltonian. Therefore, in this section we want to demonstrate how naturally all lower lying QEDFT reformulations are just approximations to the fully relativistic QEDFT that we presented in the previous sections. In lowest order we then recover the two-site Hubbard model coupled to one mode of Sec. II.

A. Equations of motion in the non-relativistic limit
In this subsection we derive the non-relativistic limit of the basic equations of motion, on which the QEDFT reformulations are based. We show how approximations in the Hamiltonian correspond to approximations in the basic equations of the corresponding QEDFT approaches.
Let us first start with the non-relativistic limit of the fully coupled QED Hamiltonian in Coulomb gauge. From the Heisenberg equation of motion, defininĝ and α k = γ 0 γ k , we find the quantized Dirac equation (in the Heisenberg picture) and accordingly forψ † . We see that the electronic componentsφ of the four-spinor are mixed with the positronic componentsχ. Of course, for small energies only the electronic component of the four-spinor is important, and therefore we would like to find an equation based solely onφ. So naturally we would like to decouple the upper componentφ from the lower componentχ. A possible way would be to find a unitary transformation of the Dirac Hamiltonian that does this, at least perturbatively. A possible expansion parameter for auch a perturbative transformation would be (mc 2 ) −1 , since we know that the energies involved in non-relativistic processes are small compared to the rest-mass energy. This energy also represents the spectral gap between the electronic and positronic degrees of freedom, which effectively decouples the dynamics of the particles and anti-particles for small enough energies. The resulting unitary transformations are known as the Foldy-Wouthuysen transformations [75] and are routinely used to generate the non-relativistic limits of the Dirac equation to any order desired. Here, we employ an equivalent but different procedure to decouple the electronic from the positronic degrees of freedom. To do so, we first rewrite Eq. (71) componentwise where we defined And thus we (formally) find that If we assume non-relativistic energies, the main contri-bution to the energy of the system stems from mc 2 , i.e. i c∂ 0 ≈ mc 2 . Accordingly, from the Neumann series of the resulting operator we find the inverse operator to lowest order as D (x) +mc 2 −1 ≈ 1/2mc 2 , and consequentlŷ At this level of approximation to the full QED problem we find the Pauli-Fierz Hamiltonian (already transformed back to the Schrödinger picture), where the non-relativistic kinetic energy reads aŝ the energy of the electromagnetic field is given as before, the Coulomb energy is given bŷ and the non-relativistic current is defined bŷ Here we used the definition of the paramagnetic current the magnetization densitŷ and the zero component of the current By construction the current obeys the continuity equation ∂ 0Ĵ0 (x) = − ∇ ·ˆ J(x). We note here, that due to the non-relativistic limit the physical current defined in Eq. (75) becomes explicitly time-dependent [5]. Further, we point out that the result of the above (formal) derivations is the same as the result obtained by first performing the non-relativistic limit of the classical Hamiltonian H QED (t) (constructed from the classical Lagranian density of Eq. (B1)) and then canonically quantizing the Schrödinger field, as shown in Fig. 6 (a). In the non-relativistic limit the resulting Hamiltonian commutes with the particle-number operatorN = d 3 rφ † ( r)φ( r), as can be seen directly from the continuity equation. Accordingly we do not need to employ a field-theoretical description for the electrons and all matter-operators can be expressed in first-quantized notation (while still being a many-particle problem). Nevertheless, we can still encounter infinities due to the interaction between the non-relativistic particles and the quantized Maxwell field [73,74]. While we do no longer have vacuum polarization (no electron-positron pairs are possible) and vertex corrections, we still have an infinite self-energy [73]. To first order in the coupling the ground-state energy (for a ext = j ext = 0) diverges as where Λ is the ultra-violet cut-off for the photon modes. By subtracting the infinite self-energy of the groundstate, which is equivalent to introducing a renormalized mass, we can renormalize the Hamiltonian perturbatively. In the following we assume, that the Pauli-Fierz Hamiltonian can be fully renormalized. For instance, in the limit of only scalar photons (the Nelson model) we know that we can perform a full renormalization of the Hamiltonian by subtracting the self-energy (provided that the kinetic energy of the problem is smaller than mc 2 ) [70,71]. Therefore, we interpret the electron mass in the Hamiltonian as a bare mass, i.e. we subtract the infinite self-energy. Now, the equation of motion forĴ k can be either found by the Heisenberg equation with the Pauli-Fierz Hamiltonian or by the non-relativistic limit of Eq. (43) (see appendix (C)). We have explicitly checked both ways of performing the non-relativistic limit as schematically indicated in Fig. 6 (b). After some calculation we find (omitting the spatial and temporal dependences) wherê is the usual momentum-stress tensor and is the interaction-stress force (the divergence of the interaction-stress tensor) [5,9,76]. If we would have started with an uncoupled problem, we would find a similar equation with the replacementˆ A tot → a ext and W k → 0. Further, the equation for the electromagnetic field does not change, except that we now have to employ the non-relativistic current (see appendix C).
In a next step we perform the non-relativistic limit for the equation of motion of the polarization, i.e. Eq. (44). We find to order 1/mc 2 Thus, at this level of approximation the polarization does not change in time.

B. QEDFT for the Pauli-Fierz Hamiltonian
In this subsection we derive the basic formulation of non-relativistic QEDFT for the full Pauli-Fierz Hamiltonian. We show how the Gordon-decomposition, i.e. the equation of motion for the polarization P k , makes the current J k a unique functional of (a k ext , j k ext ) and thus becomes the basic variable for the matter part in this limit. Further we demonstrate how the non-relativistic limit of the above Kohn-Sham construction produces the Kohn-Sham construction for the Pauli-Fierz Hamiltonian. A comparison of this level of approximation with relativistic QEDFT and with other approximations is presented schematically in appendix F.
We start by performing the non-relativistic limit of Eq. (60). Irrespective of the difference between a ext and ã ext (note, that we again employ the a 0 ext = 0 gauge for the external potentials as explained in Sec. III B) the equation in this limit is always zero. However, by employing Eq. (77) we can rearrange the non-relativistic limit to which is non-zero provided the density obeys J Hamiltonian in Eq. (74) accordingly. Thus, e.g., by assuming a negligible magnetic density M l (x) ≈ 0, i.e.
the according Hamiltonian as well as the defining Eqs. (76) and (77) change. Actually, all termsM l and q M l vanish in these equations for this approximation. We again find due to the corresponding Eq. (78) that we have (a k ext , j k ext ) and we can consider the corresponding coupled Eqs. (80) and (81). The Kohn-Sham current becomes accordingly j k KS = j k ext +J k and the Hxc potential in this limit reduces to On the other hand we can also restrict the allowed photonic modes. For instance we can assume a perfect cubic cavity (zero-boundary conditions) of length L 9 . Then, with the allowed wave vectors k n = n(π/L) and the corresponding dimensionless creation and annihilation operatorsâ † n,λ andâ n,λ (see appendix (E) for more details) we find where the mode function S is given in Eq. (D1). If we further restrict the modes by introducing a squaresummable regularization function f EM ( n) 10 , e.g. f EM = 1 for | n| < mcL/(2π ) (energy smaller than rest-mass energy) and 0 otherwise, the resulting regularized field makes the coupled Pauli-Fierz Hamiltonian self-adjoint without any further renormalization procedure [74]. Such a restriction is assumed in the following. These approximations are then directly reflected in the Hamiltonian and the derived equations of motion. While the basic Eq. (77) does not change, and thus J k is the basic mattervariable, the basic equation of motion for the potential A k has to reflect the restriction to specific modes. By multiplying Eq. (81) from the left by k ( n, λ)S( n · r) and integrating we find the mode expansions whereq n,λ =â n,λ +â † n,λ and we use the definition j ext n,λ (t) = d 3 r ( n, λ) · j ext (x) S( n · r).
The Coulomb part vanishes since we employ a partial integration and the fact that ( n, λ) · n = 0. Of course, one finds the same equations by a straightforward calculation of the Heisenberg equation of motion for the Maxwell-field (85) with the according Pauli-Fierz Hamiltonian (74). From the restriction to specific modes the field A k is restricted in its spatial form and therefore the photonic variable changes from A k to the set of mode expectation values This change in basic variable is also reflected in the conjugate external variable which is given via Eq. (86) as Thus we accordingly find j k ext (x) → j ext n,λ (t) , and the conjugate pairs become (a k ext , j ext n,λ ) Thus we have to solve the mode Eqs. (86) together with the according equation of motion for the current. Correspondingly also the Kohn-Sham scheme and the mean-field approximation for a Hxc change to its modeequivalents.
If we then also employ the dipole-approximation, i.e. we assume that the extension of our matter system is small compared to the wavelengths of the allowed photonic modes, we find This only changes the definition of effective currents that couple to the modes, i.e. j ext n,λ (t) = d 3 r L 3/2 ( n, λ) · j ext (x), but leaves the structure of the QEDFT reformulation otherwise unchanged. If we assume the magnetization density M l to be negligible, we have from first principles rederived the QEDFT formulation presented in [40]. In this work the situation of only scalar external potentials, i.e. a ext = 0 and a 0 ext = 0, has been considered as a second case. In this situation, the gauge freedom is only up to a spatial constant, which is usually fixed by choosing a 0 ext → 0 for | r| → ∞. Since a 0 ext couples to the zero component of the current, i.e. the densityĴ 0 , the conjugate pair becomes (a 0 ext , j ext n,λ ) To demonstrate this mapping, the first time derivative ofĴ 0 is obviously not enough, since this amounts to the continuity equation and no direct connection between the two conjugate variables of the matter part of the quantum system is found. Therefore, one has to go to the second time derivative ofĴ 0 [40]. If we then further simplify this physical situation (see appendix (E) for a detailed derivation) we find the model Hamiltonian of Sec. II. In a similar manner, by imposing the restrictions on the corresponding equations of motion, we can rederive the model QEDFT of Sec. II B 11 . Finally, for a simple overview, we have collected the different QEDFTs that we have explicitly considered in this work in appendix F.

V. CONCLUSION AND OUTLOOK
In this work we have shown how one can extend the ideas of TDDFT to quantized coupled matter-photon systems. We have first explained the basic ideas of QEDFT for a model system of a two-site Hubbard model coupled to a single photonic mode. By rewriting the problem in terms of an effective theory for a pair of internal functional variables and proving the uniqueness of solutions for the resulting non-linear coupled equations, we have demonstrated how an explicit solution for the coupled photon-matter wave function can be avoided. Further we have discussed how an auxiliary quantum system, the so-called Kohn-Sham system, can be used to construct approximations for the implicit functionals appearing in the effective equations. The Kohn-Sham construction gives rise to effective fields and effective currents, which are termed Kohn-Sham potential and Kohn-Sham current, respectively. By numerically constructing the exact Kohn-Sham potential and Kohn-Sham currents, we have 11 We note, that one could have also derived the model QEDFT by employing the gauge of Eq. (47) for the external potentials. By applying the dipole approximation also to the external vector potential the conjugate variable becomes the density (dipole moment). For clarity of presentation, though, we have chosen to start from the scalar-potential case.
illustrated the capability of this new approach to exactly describe the dynamics of coupled matter-photon systems and contrasted these exact fields with the mean-field approximation.
In the following, instead of reformulating every possible approximate treatment of coupled matter-photon systems seperately, we have shown how these QEDFTs for approximate Hamiltonians are merely approximations to relativistic QEDFT, which itself is based on QED. To avoid problems with the Kohn-Sham construction, we have based relativistic QEDFT on the expectation value of the polarization and the vector potential of the quantum system. By then performing the nonrelativistic limit of QEDFT we have demonstrated that the resulting theory is the QEDFT reformulation of the Pauli-Fierz Hamiltonian. The non-relativistic limit automatically makes the (non-relativistic) current the basic variable for the matter system. Accordingly, the nonrelativistic limit of the Kohn-Sham potentials and currents leads to the corresponding Kohn-Sham fields for the Pauli-Fierz Hamiltonian. By performing further approximations for non-relativistic QEDFT, e.g. assuming the magnetic density negligible, we have shown how other QEDFTs ( that reformulate the corresponding approximate Hamiltonians) can be derived. Depending on the level of approximation, the basic internal functional variables change, e.g. if we confine the electromagntic field with a cavity, the (allowed) mode expectation values become the new internal variable of the photons. In a final step we restricted to a two-site model coupled to only one mode, recovering the model QEDFT of the beginning.
We point out, that at every level of QEDFT we recover the corresponding (standard) time-dependent densityfunctional reformulations [47,72] if we assume the quantized nature of the photons negligible, i.e. the charged particles interact via the classical Coulomb interaction only. This will be the case in most standard situation of condensed matter theory, e.g. when investigating dynamics of atoms or molecules in free space. However, we expect that interesting effects happen when the boundary conditions for the Maxwell field are changed, e.g. for atoms in a cavity. Thus we have a potential tool that can treat complex electronic systems in the setting of quantum optics. Also, we can investigate the explicit interplay of photons with molecules or nanostructures, e.g. in nano-plasmonics. However, for this theory to be practical we are in need of reliable approximations to the basic functionals. In [67] functionals based on an optimized effective potential approach [9,10] are constructed, which provide good results even in the situation of strongly-coupled systems. Although the currently available approximations have only been tested for simple model systems, the hierarchy of QEDFT approximations allows to simply scale up these functionals to more complex situations. Thus we can develop approximations for simple systems, e.g. only one mode couples to the matter system, and then extend these approximations to more involved problems, e.g. considering more modes.
In this way we can easily control the validity of our approximations. In this respect, we are also working on a fixed-point approach in the spirit of [46,59], which allows us to construct the exact Kohn-Sham potentials and compare the approximate potentials to the (numerically) exact expressions. Details of this approach will be part of a forthcoming publication [60]. On the other hand, the fixed-point approach is also a way to extend the validity of QEDFT beyond Taylor-expandable fields. A different way, especially for discretized matter systems, is the non-linear Schrödinger equation approach introduced in [45,49]. Certain theoretical and mathematical details of the model QEDFT of Sec. II, that are beyond the scope of the current manuscript, will be discussed in [77].
Finally, since we are aiming at investigating quantum optical settings, we also need to discuss the cavity and the problem of open quantum systems. In this work we focused on closed systems and on a perfect (cubic) cavity. It is straightforward (but tedious) to extend the current work to an arbitrary shape of the perfect cavity. We have to use an expansion of the photon field in the according eigenfunctions of the cavity, such that these modes obey the Coulomb-gauge condition. However, in actual quantum optical experiments, the cavities are not perfect but rather an open quantum system, which allows for an exchange with the environment. To take care of this channel of decoherence there are several possible ways. One can employ the current formulation of QEDFT and derive a master equation, as has also been done for standard TDDFT [78,79]. Also extensions to stochastic equations [80][81][82] are possible. On the other hand, one can couple further bosonic degrees of freedom to the system and prescribe a bath spectral density, making these degrees of freedom a bath for the system [40]. Since the present framework allows for a consistent treatment of interacting fermionic and bosonic particles, the inclusion of a bath and coupling to other fields, e.g. phonons, will be the subject of future work.

Appendix A: Conventions
In this work we employ the standard covariant notation x µ = (ct, r) with greek letters indicating four vectors, e.g. µ ∈ {0, 1, 2, 3}, and roman letters indicating spatial vectors, i.e. k ∈ {1, 2, 3}. To lower (or raise) the indices, i.e. going from contravariant vectors to covariant vectors (or vice versa), we adopt the convention for the Minkwoski metric. We denote spatial (contravariant) vectors with the vector-symbol, i.e. A k ≡ A, and the derivatives with respect to the space-time vectors x µ by ∂ µ = ∂/∂x µ . With these definitions the divergence can be written as ∂ k A k = ∇ · A, where we also adopt the Einstein summation convention. Further we note that With the help of the Levi-Civita symbol ijk we can write the curl as ijk ∂ j A k ≡ − ∇ × A and the multiplication of Pauli matrices becomes σ k σ l = Further, for notational simplicity we only point out in the text (when necessary), whether we are in the Schrödinger or Heisenberg picture, and do not explicitly indicate the picture used in the operators. In the Schrödinger picture, operators which are not explicitly time-dependent only carry a purely spatial dependence, e.g.Â k ( r). We indicate explicit time-dependence in the Schrödinger picture by either carring the full space-time dependence, e.g.Ĵ k (x) (for the Pauli-Fierz current density) of Eq. (75), or by a dependence on t, e.g.Ĥ(t) in Eq. (74). In the Heisenberg picture, every operator also depends on time, e.g.ψ(x) in Eq. (71).

Appendix B: Quantum Electrodynamics in Coulomb gauge
In this appendix we give a detailed derivation of QED in Coulomb gauge. We start from the (classical) coupled QED Lagrangian with external fields a ext µ (x) and j ext µ (x) given by [2] Here we use the standard definitions for the (classical) Dirac fields, i.e.
is a Dirac four-spinor with the two-component (spin) functions φ(x) and χ(x), the Gamma matrices are given by with σ i the usual Pauli matrices,ψ =ψγ 0 and is the conserved (Noether) current. Further we use the Minkowski metric g µν = (+, −, −, −) to raise and lower the indices. For the (classical) Maxwell field we have where is the electric field tensor and A µ (x) is the vector potential. Now we employ the Coulomb gauge condition for the Maxwell field, i.e. ∇ · A(x) = 0. Then it holds that where ∆ is the Laplacian. If we impose squareintegrability on all of R 3 12 the Green's function of the Laplacian becomes ∆ −1 = 1/(4π| r − x |) and therefore Since the zero component of the four potential A µ (x) is given in terms of the full current, it is not subject to quantization. The conjugate momenta of the photon field (that need to be quantized) are the same as in the current-free theory and thus the usual canonical quantization-procedure applies [2], i.e.
whereÊ k is the electric field operator, δ ⊥ kl ( r − r ) = (δ kl − ∂ k ∆ −1 ∂ l )δ 3 ( r − r ) is the transverse delta-function and k, l are spatial coordinates only. Equivalently we can define these operators by their respective plane-wave expansionŝ ( k, λ) â k,λ e i k· r −â † k,λ e −i k· r , 12 If we consider the situation of a finite volume, e.g. due to a perfect cavity, the boundary conditions change. These different boundary conditions, in principle, change the Green's function of the Laplacian and thus the instantaneous interaction. We ignore these deviations from the Coulomb interaction in this work for simplicity.
If we further define the magnetic field operator by cˆ B = ∇ ×ˆ A, the Hamiltonian corresponding to L E is given in Eq. (35). We used normal ordering to get rid of the infinite zero-point energy in this expression. Also for the Dirac field, the coupling does not change the conjugate momenta. Therefore we can perform the usual canonical quantization procedure for Fermions which leads to the (equal-time) anti-commutation relations [2] {ψ α ( r),ψ β ( r )} = γ 0 αβ δ 3 ( r − r ). The Hamiltonian corresponding to L M thus becomes the one of Eq. (34), where we used r · y = −x k y k .
It is straightforward to give the missing terms of the QED Hamiltonian due to the coupling to the external fields as well as due to the coupling between the quantized fields. If we apply the definition of the quantized current J µ of Eq. (37) to Eq. (B4) we find (using normal ordering, i.e. rearranging the annihilation parts of the operators to the right, to discard the corresponding zero-point energy) : J 0 ( r)+j ext 0 (x) A 0 (x) := 1 2c +2Ĵ 0 ( r)j 0 ext (x )+ :Ĵ 0 ( r)Ĵ 0 ( r ) : .
Adding and subtracting on the right hand side the term eφ † σ k i c∂ 0 −D χ and employing Eq. (72) we find With the help of the definition [...] = [D + mc 2 ] −1 this can be rewritten as Now, if we employ the approximation [...] ≈ 1/2mc 2 (also in the Coulomb terms) we end up with which is just the equation of motion for the nonrelativistic current (75) with the Pauli-Fierz Hamiltonian. For the Maxwell field the non-relativistic limit of Eq. (49) is with the help of Eq. (75) straightforward. It is only important to see, that this does agree with the equation of motion forÂ k due to the Pauli-Fierz Hamiltonian (74). The main difference to the fully relativistic derivation is, that now we have a term of the form e 2mc 2 d 3 rĴ 0 (x) Â k (x) + a k ext (x) Â k (x) + a ext k (x) .
This term does not change anything in the first order equation, i.e. ∂ 0Âk = −Ê k . In the second order we find due to Eq. (B5) that Now, with the above definition for ∆ −1 used in Eq. (B4) we find that these commutators lead to the terms of the equation of motion for the Maxwell field in the nonrelativistic limit. The rest of the derivation is similar to the relativistic situation.

Appendix D: Mode expansion
If we restrict the allowed space for the photonic modes we also need to impose according boundary conditions. Let us first start with a cubic cavity of length L with periodic boundary condition. We then find with the allowed wave vectors k n = n(2π/L) and the corresponding dimensionless creation and annihilation operatorsâ † n,λ and a n,λ , which are connected to their continuous counterparts by lim L→0 L 3/2â n,λ →â k,λ , that A k ( r) = c 2 0 L 3 n,λ k ( n, λ) √ 2ω n â n,λ e i kn· r +â † n,λ e −i kn· r .
where t kin is the kinetic (hopping) energy, l is the vector connecting two sites, and a 0 ext (t) corresponds to the potential difference between the sites. To highlight the general structure of the photon-matter Hamiltonian (and bring it to the form used in Sec. II) we also redefine the external current, the external potential, and the photon field as follows ∂ t j ext (t) → ωj ext (t), After implementing the above redefinitions in Eq. (E4) and neglecting irrelevant constant terms we arrive at the following Hamiltonian We note, that the same model Hamiltonian could have been derived by assuming an external vector potential in a gauge such that a 0 ext = 0 and a ext = 0. Then by the dipole approximation the corresponding Hamiltonian to Eq. (E2) we would have terms of the form a ext · ∇, a 2 ext and mixed terms of internal and external vector potential. By going into length gauge also for the external potential and performing the same steps as above, one ends up with the same two-site one-mode Hamiltonian.

Appendix F: Overview of QEDFTs
Here we give an overview of the different QEDFTs that we have discussed explicitly. We employ for the Kohn-Sham scheme an uncoupled auxiliary quantum system with an initial state | Φ 0 = | M 0 ⊗ | EM 0 . For the different levels of approximation the prerequisites for this initial state change, i.e. we might have different initial conditions that have to be fulfilled. Further we use the notational convention that the super-index s refers to the (uncoupled) Kohn-Sham quantity, e.g. P 0 [Φ 0 , P k , A k ] = P s 0 . are coupled to the charged quantum particles, which effectively also leads to a coupling between the photons. This can be most easily seen in the non-relativistic limit, where the inhomogeneity contains terms like Ĵ 0Âk . Since the current of the auxiliary Kohn-Sham system is by construction equal to the exact current (at least for the non-relativistic limit), this coupling between the photons is also present in the Kohn-Sham Maxwell equation. The term J s 0 a Hxc of the Kohn-Sham current contains these non-trivial couplings as functionals of the initial states and internal pair. When restricting the photons to a cavity, the Kohn-Sham current is then responsible to couple the different photon modes. The coupling terms in the Kohn-Sham current are specifically relevant in the context of, e.g. nano-plasmonics, where the electromagnetic fields are enhanced due to the presence of the plasmons, or in the optical control of currents in solids [84].