Existence, Uniqueness, and Construction of the Density-Potential Mapping in Time-Dependent Density-Functional Theory

In this work we review the mapping from densities to potentials in quantum mechanics, which is the basic building block of time-dependent density-functional theory and the Kohn-Sham construction. We first present detailed conditions such that a mapping from potentials to densities is defined by solving the time-dependent Schr\"odinger equation. We specifically discuss intricacies connected with the unboundedness of the Hamiltonian and derive the local-force equation. This equation is then used to set up an iterative sequence that determines a potential that generates a specified density via time propagation of an initial state. This fixed-point procedure needs the invertibility of a certain Sturm-Liouville problem, which we discuss for different situations. Based on these considerations we then present a discussion of the famous Runge-Gross theorem which provides a density-potential mapping for time-analytic potentials. Further we give conditions such that the general fixed-point approach is well-defined and converges under certain assumptions. Then the application of such a fixed-point procedure to lattice Hamiltonians is discussed and the numerical realization of the density-potential mapping is shown. We conclude by presenting an extension of the density-potential mapping to include vector-potentials and photons.


General overview
The Schrödinger equation [1] (together with its variants, e.g., [2,3]) is ubiquitous in physics, chemistry, material science, and biology as it describes in detail the interactions between electrons and atomic nuclei which are the building blocks of atoms, molecules, and solids. These interactions completely determine the physical and chemical properties of atomic, molecular, and condensed matter systems and a deep understanding of these properties therefore requires the solution of the Schrödinger equation. This is, however, a very difficult problem for realistic systems due to the Coulomb interaction between the electrons which prohibits the decoupling of the many-electron Schrödinger equation [4,5] into single-particle problems, which can be solved efficiently on modern computers (see e.g. [6,7]). Therefore, in principle, we have to treat the huge number of degrees of freedom of an interacting manybody system explicitly, which can only be done for simple, i.e., small systems. This exponential increase of complexity with the number of interacting particles is known as the exponential wall [8]. Several approaches (see e.g. [4,5,9]) have been developed that try to avoid this exponential scaling by considering reduced quantities instead of the full many-body wave function. For the description of electronic systems in their ground state one of the most successful of these approaches [10] is density-functional theory (DFT) [11,12], which allows to determine the exact ground-state observables by only knowing the oneparticle density. The foundation of ground-state DFT is the (time-independent) density-potential mapping that was first established in the seminal paper of Hohenberg and Kohn [13]. By applying the Rayleigh-Ritz minimal principle of quantum mechanics, they could show that there exists a one-to-one correspondence, i.e., a bijective mapping, between the set of ground-state densities and their respective external scalar potentials. Densities which are connected via the solution of a Schrödinger equation to an external potential are termed v-representable. The second cornerstone of DFT is the Kohn-Sham construction that allows to determine the density of an interacting quantum system by considering an auxiliary non-interacting system [14] (and thus decoupling the problem into single-particle problems). The Kohn-Sham construction actually employs a composition of the Hohenberg-Kohn mapping of an interacting and a non-interacting system. To be able to do so, one has to assume that the set of interacting and non-interacting ground-state densities is the same, i.e. that every interacting density is also non-interacting v-representable [11,12]. That this is true has been shown under certain restrictions [15,16,17,18], but in the most general situation [19] this is still an open issue b. However, what has been proven so far [15,16,17,18] already provides a sound theoretical basis for numerical implementations of the Kohn-Sham method which is gratifying since it is currently the most widely used electronic structure method in solid state physics and quantum chemistry [10].
The discussion so far was on ground-state properties. However, for the description of dynamical properties of electronic systems an extension of the ground-state formalism of DFT is required. The first such extension of DFT to dynamical systems was discussed by Peuckert [20]. This work assumed the existence of a time-dependent density-potential mapping, but did not present a formal proof of a bijective mapping between a set of time-dependent external scalar potentials and their respective time-dependent one-particle densities. The formal justification for time-dependent density-functional theory (TDDFT) [21,22] was presented in [23], where Runge and Gross showed a one-to-one correspondence between time-analytic potentials and their respective time-dependent densities based on the (divergence of the) local-force equation of quantum mechanics [24]. Similar to the ground-state case the construction of a corresponding Kohn-Sham scheme requires that the set of interacting v-representable densities is the same as the set of non-interacting vrepresentable densities. A first proof of the existence of a time-dependent Kohn-Sham scheme [25] has also been based on the local-force equation (and the assumption of time-analytic potentials as well as densities). Again, the mapping from a timedependent interacting to an auxiliary non-interacting system is the composition of two density-potential mappings and demands that the set of densities does not depend on the specific interaction potential. We further like to mention that an alternative proof of the Runge-Gross theorem exists in the linear response regime. One can prove the invertibility of the density response function for perturbations from a non-degenerate ground state [26,27]. This relaxes the constraint of having Taylor expandable external potentials and only requires that the Laplace transform in time of the potentials exist. However, since in this review we want to deal with the density-potential mapping in its most general context we will not consider this more specific case.
Since TDDFT is a younger field of research its mathematical foundations are not as established as that of ground-state DFT. Also the required mathematical proofs are of a different nature. Ground-state DFT singles out a specific state, namely the ground state, which can be obtained from a minimum principle. Many of the proofs are therefore based on certain properties (such as convexity) of energy functionals which are peculiar to the ground-state problem. The basis of TDDFT on the other hand is formed by the time-dependent Schrödinger equation which describes an initial-value problem, i.e. for a given initial state its time-evolution is to be sought. Mathematical proofs must therefore be based on evolution equations. The ground state does not play any special role except that it can be chosen as the initial state of a time-evolution.
In the following we will give an extensive review of the density-potential mapping in TDDFT, which, as explained in the preceding short historical introduction, forms the basis of TDDFT and the time-dependent Kohn-Sham approach. We will not discuss any approximations that are needed to perform TDDFT in practice, but will rather concentrate on what is known about the exact properties of the one-to-one correspondence between densities and potentials. In the following we will first set the stage for the closer investigation of the density-potential mapping. We highlight how the density appears as a fundamental and useful quantity in time-dependent quantum mechanics in Sec. 1.2. Assuming the existence of a density-potential mapping, in Sec. 1.3 we demonstrate how Kohn-Sham TDDFT can help to avoid the exponential wall and approximately determines properties of large quantum systems. Before we consider the direct mapping of potentials to densities we first illustrate in several examples the intricacies of the time-dependent Schrödinger equation in Sec. 2.1. This section will illustrate that we carefully need to take into account the unboundedness of the Schrödinger Hamiltonian if we want to discuss the density-potential mapping. In Sec. 2.2 we give a short overview of the history of the mathematical treatment of explicitly time-dependent initial-value problems, before we define exactly what we mean by a solution of the time-dependent Schrödinger equation in Sec. 2.3. Next we give precise conditions for the existence of such solutions in Sec. 2.4. We will then properly define the potential-density mapping in Sec. 2.5. After the definition of v → n we will give an analytical example of the inverse mapping n → v in Sec. 3.1 and discuss certain implications. In Sec. 3.2 we present an iterative procedure to construct the density-potential mapping in the most general case. For this iteration to be possible we need that a certain operator is invertible, which will be discussed in Sec. 3.3. Based on these considerations we will present the famous Runge-Gross theorem that establishes the existence of a density-potential mapping for analytic potentials in Sec. 3.4. Then in Sec. 3.5 we extend these result to more general potentials and give conditions such that the iterative construction of the density-potential mapping converges.
How one can apply the iterative procedure in the case of lattice systems (and hence circumvent some of the mathematical problems of the continuum case) is shown in Sec. 4. Further, in Sec. 5 we present a numerical scheme that is able to construct the potential for a given density and initial state for interacting and non-interacting systems. Next we show how the ideas developed for the Schrödinger equation can be extended to systems with vector potentials and photons in Sec. 6. Finally we conclude with a summary and outlook in Sec. 7.
1.2. The time-dependent many-particle problem and the density-potential mapping Before we go in much more detail into the properties of the density-potential mapping we want to describe here in general terms the main ideas and motivations for the consideration of such a mapping. We start with a discussion of the time-dependent Schrödinger equation (TDSE). The TDSE of an N -electron system is a partial differential equation of the form where x is the collection of spatial and spin variables of the system, i.e. x = (x 1 , . . . , x N ) where x j = r j σ j is a space-spin variable. The wave function Ψ 0 (x) is the prescribed state at the initial time t 0 , which is often taken to be an eigenstate of a time-independent Hamiltonian. The Hamilton operatorĤ(t) contains all the information on the system and has the structurê whereT represents the kinetic energy operator,V (t) the time-dependent external potential andŴ the two-body interactions. Their explicit form in first quantization is (in atomic units) The two-body interaction w is in realistic applications almost always the Coulomb potential w(r) = 1/|r| (for the mathematical considerations in this review it should at least have the property that w(r) = w(−r)). The external potential v(r, t) is typically the sum of a static and a time-dependent part and depends on the physical system of interest. For instance, for a general molecule with M fixed atomic nuclei it is given by where Z α and R α are the charge and position of atomic nucleus α and v ext (r, t) is an externally applied field, for example a laser pulse (which can be approximately described by a scalar potential in some suitable approximation). The many-body problem is now completely defined. The task is to solve the TDSE and once the wave function is obtained we can calculate all physical observables of interest. The central problem which remains in practice, is to do this in an preferably efficient manner for realistic systems of physical interest. This is also often referred to as the quantum many-body problem.
An important observation is now that when, for instance, studying different molecules with the same number of electrons we only change the potential v but always keep the functional form of the two-body interaction and the kinetic energy operator the same. Any system of interest in this setting is fully specified by the external potential and the initial state, i.e., the wave function that we calculate from the TDSE is a functional of v for a given initial state Ψ 0 . We can therefore talk about a mapping v → Ψ[v] from potentials to wave functions for a given initial state. This idea immediately raises the question for which class of potentials the initial-value problem of the TDSE has an acceptable solution. One can, for instance, imagine that for singular enough potentials (for instance making the Hamilton operator unbounded from below) a solution may not exist or develop undesirable properties such as infinite expectation values for certain physical observables such as the kinetic energy. This is a question about the proper domain of the mapping v → Ψ[v] (see Sec. 2.4). For the moment it will suffice that we consider a domain of physically reasonable potentials and refer for a more detailed discussion to later sections of this review. If for a given set of potentials we can construct the many-body state Ψ[v] then we can also construct any observable of interest as a functional of v. In particular ifÂ is any operator representing a physical quantity then its expectation value is a functional of the external potential v. It is clear that A([v], t) is a functional of the potential v but the converse is, in general, not true as A ([v], t) could be an object of lower dimensionality than the potential v(r, t). For instance, there could be many different potentials v(r, t) that all generate the same time-dependent dipole moment d(t) of some molecule. For a given initial state Ψ 0 the knowledge of the function A([v], t) is, in general, not sufficient to determine the potential v. There is, however, a natural variable for which such an inversion is possible, namely the time-dependent density. The reason is that this observable is related in a special way to the potential. Let us explain this in more detail. We can rewrite (4) aŝ is called the particle density or simply density. This is a physical quantity for which n(r, t)dr gives the probability to find a particle in a small volume dr around point r.
More explicitly this can be written as where the volume element dx N −1 refers to a integration of the spatial coordinates and summation over the spin coordinates of N − 1 particles. Which of the coordinates these are is not relevant if we take |Ψ(x, t)| 2 properly symmetrized under interchange of particles. The coordinates that we integrate over we can label by x 2 , . . . , x N and the space-spin-coordinate not included in the volume element dx N −1 we will call x = (r, σ). Due to special form of (5) an equation of motion, known as the local-force equation, can be derived that directly relates the density n(r, t) and the potential v(r, t). Given this equation it can be proven that under certain conditions there exists a one-to-one mapping between densities and potentials. This equation will be discussed in detail in Sec. 2.5 and therefore we restrict ourselves at this point to an intuitive argument. We can imagine, at least for slow enough temporal variations of the potential, that by making the potential more attractive in some region of space we will increase the probability of finding particles in this region and therefore its particle density. From this intuitive picture we may conjecture that to a given density profile n(r, t) there corresponds a unique potential v(r, t) that produces it for a given initial state. It is one of the main topics of this review to specify in which sense this statement is true. As a first step we note that if we change the potential by a purely time-dependent function v(r, t) → v(r, t) + C(t) then the wave function solving the TDSE will only change by a time-dependent phase factor, i.e. Ψ → e iα(t) Ψ where α(t) = − t t0 dt C(t ), and therefore will not change any of the physical observables. Such a change of the potential is just a gauge and we will regard potentials that differ only in a gauge as physically equivalent. Therefore the more correct conjecture is that to a given density profile n(r, t) there corresponds a unique equivalence class of potentials that produces it for a given initial state. In order to specify a unique potential within this class we can always choose a particular gauge. For instance, if the class of potentials that we consider remains spatially constant at |r| → ∞ we can choose a gauge such that v(r, t) → 0 (|r| → ∞). Assume that we have made a particular gauge choice, then our conjecture implies that there exists a mapping n → v which, for a given initial state Ψ 0 , maps the density to the potential that generates it. This means that, rather than parametrizing the wave function by the potential, we could parametrize it by the density and write Ψ[n]. This implies that, in particular, the expectation value of an observableÂ will be a functional of the density and we can write These considerations form the basis of TDDFT as we will summarize in the next section and in which it will be shown that the existence of a density to potential mapping is used to derive the Kohn-Sham equations. From a more mathematical point of view the possible existence of a density to potential mapping immediately raises the new question for which prescribed density profiles such a potential can be found. In other words, what is the domain of the density-potential mapping? It would, for instance, be very useful to know if for every density profile n with certain continuity or differentiability properties there exists a potential v with certain other specified properties that generates it. Densities which are produced by a potential v by solving the TDSE are called v-representable and the problem of the characterization of the set of v-representable densities is often referred to as the v-representability problem. The elucidation of this problem is one of the main topics of this review.

Summary of time-dependent Kohn-Sham density-functional theory
This section gives a summary of the Kohn-Sham (KS) approach to TDDFT with a focus on conceptual issues. The KS equations are derived under the assumption of the existence of the density-potential mapping which will be discussed in detail later in this review. The KS formalism is heavily used in applications and the practitioner of TDDFT is in this way already introduced to a familiar framework before going deeper in the more fundamental issues of existence and uniqueness of the density functionals.
The central idea of KS DFT is to associate with an interacting system an effective noninteracting system, known as the KS system, which has the same time-dependent density as the true system. Since the KS system is a noninteracting system it is much easier to treat in numerical applications than the fully interacting system that we started out with. The price we pay for this simplification is that the one-body potential of the KS system, known as the KS potential, is an in general unknown functional of the density. Despite this difficulty it has turned out that practically useful approximations for this potential can be devised. Let us describe the KS approach in more detail. The existence of a density-potential mapping n → v for a fixed initial state Ψ 0 , which we also denote by v[Ψ 0 , n], is assumed not to depend on the chosen two-body interaction for a physically reasonable class of interactions (we will be more specific later). Specifically this means that we have a density-potential mapping for interacting as well as noninteracting systems. For the case of a noninteracting system this mapping is called v s [Φ 0 , n] (the subscript s is usually not explained in the density-functional literature but we can assume that it refers to "single particle" as the potential often appears in effective single-particle equations as explained below). Since in this case we have no two-body interactions the Hamiltonian is then simply given bŷ Let us assume that this is a closed-shell system of 2N electrons and that the initial state Φ 0 of the noninteracting system is chosen to be a single Slater determinant of orbitals ϕ j,0 (r). This allows us to reduce the TDSE for the noninteracting system to single-orbital equations of the form where the initial conditions are given by ϕ j (r, t 0 ) = ϕ j,0 (r). For a given density profile n(r, t) we have to construct a proper initial state Φ 0 producing the initial density n(r, t 0 ) (this can be done in practice, for example, using the so-called Harriman construction [28]) and the initial current [29,30]. Once this is done our assumption guarantees the existence of the potential v s [Φ 0 , n]. If the density that we prescribe happens to be the density of an interacting system then we have succeeded in reproducing this density within a noninteracting framework. However, this is not relevant in practice, since the (8) and (9) do not allow us to predict the density of the interacting system as they obviously do not have any information on which interacting system we would like to solve. To make a predictive scheme we have to connect the true and the KS system. To do this we introduce the KS potential where v ext it the external potential of the interacting system of interest. Note that here we make a careful distinction between v s and v KS as they have different functional dependencies. This is usually not done in the density-functional literature where both quantities are often denoted by the same symbol v s which can lead to confusions. If we assume the full knowledge of the functionals v[Ψ 0 , n] and v s [Φ 0 , n] then the set of equations does have a unique solution for a self-consistent density n sc . By definition of v s [Φ 0 , n] this self-consistent density is exactly obtained whenever which according to (10) is precisely satisfied when In turn, this is exactly true when n sc is equal to the density produced by the potential v ext in the interacting system with initial state Ψ 0 , which is precisely the density that we are interested in. The procedure that we outlined here therefore comprises a predictive computational scheme to calculate the density of an interacting system of interest with the single-orbital equations (11) and (12). To make the scheme useful in practice we need an approximation to the functional v s [Φ 0 , n] − v[Ψ 0 , n] appearing in (10). Usually this functional is split into two pieces as follows where v H is the Hartree potential defined as v H ([n]; r, t) = dr w(r − r ) n(r , t), which describes a mean-field classical electrostatic potential between the electrons, and v xc is the so-called exchange-correlation (xc) potential. By this redefinition we have shifted all difficult functional dependencies to the xc potential. Putting everything together we recover the standard KS equations as they appear in textbooks although most textbook discussions are less precise and do not indicate the initial state dependencies of the xc potential explicitly. The main obstacle for applying these equation to the calculation of electronic properties is obtaining a good approximation for the xc potential. Several useful approximations for this quantity have been developed. The discussion of such approximations is, however, not the main aim of this review. For a comprehensive overview of approximate xc potentials and their applications we refer to a recent textbook on TDDFT [22].

Overview of the time-dependent Schrödinger equation
The very foundation of the density-potential mapping is that for a given initial state the solution of the TDSE for two different external potentials (differing more than a gauge) leads to two different densities, i.e. the mapping from potentials to densities is injective. However, to investigate the basic properties of this mapping, we first need to properly define it. To do so, we need to specify the domain of the mapping, i.e., the set of potentials, the codomain or range of the mapping, i.e., the set of densities, and the rule how the potentials are mapped to their respective densities. It is obvious that we want to have a domain such that for every potential in this set we can actually solve the TDSE uniquely. Therefore we want to investigate under which conditions the TDSE has a unique solution. Although in general this is tacitly assumed, for a proper investigation of the density-potential mapping we need to know specifics. Thus, in Sec. 2.3 we discuss what we mean by a solution to the TDSE and set the stage for a precise presentation of the density-potential mappings in TDDFT. We introduce the notion of classical and non-classical solutions to the initial-value problem. These non-classical solutions arise since the Hamiltonians in quantum mechanics are usually unbounded operators and thus can lead to infinities.
To make sense of these generalizations we also need to present some (in physics often ignored) details about self-adjoint operators. We briefly discuss the idea of a selfadjoint domain and give conditions on the two-body interactions and the potentials such that the resulting Hamiltonians are self-adjoint on a common domain. We then present conditions (on the external potentials) for the existence of unique solutions to the TDSE. Regularity properties, e.g., under which conditions we have classical solutions, are discussed as well. In a next step we then investigate physical quantities derived from the wave functions. We give exact conditions such that the density obeys the continuity equation, and discuss which restrictions we need to impose on the potential in order to obey the fundamental equation of TDDFT. Since the rigorous discussion of the TDSE involves some abstract concepts we will first illustrate these concepts by discussing a simple but very common physical situation, namely the free propagation of a wave packet.

Free propagation of a wave packet
Rather than formally discussing the TDSE at this point let us first point out some issues one stumbles upon when trying to solve the equation in practice. These problems then will force us later to adopt a more careful and rigorous approach. We will start to consider the simple case of one particle in one dimension enclosed in a box of length L. The corresponding Hilbert space H is given by the square integrable functions in the box, or more formally H = L 2 ([0, L]) with standard inner product between functions Ψ and Φ. Using the inner product we can assign to any function Ψ in the Hilbert space the norm Ψ = Ψ|Ψ . Let us now start by considering a simple free evolution of a single particle wave packet. This means that we want to calculate a state Ψ(t) ∈ H at time t from a given initial state Ψ(0) ∈ H. The TDSE (in atomic units) for this problem is given by and we specify the initial state Ψ(0) at time t = 0. We have not specified yet any properties of the initial state Ψ(0) but it is reasonable to assume that it is twice differentiable such that the action of the HamiltonianĤ on it is well-defined (further conditions will follow soon). A formal solution of (20) is given by where the exponent of an operator is formally defined by its Taylor series. This, of course, assumes that the infinite series converges in some norm sense, which we did not check at this point (and in fact turns out to be false in general). One sees immediately that problems arise with (21) whenever the function Ψ(x, 0) is only a finite times differentiable so let us assume that Ψ(x, 0) is infinitely differentiable on the real line (more technically Ψ(x, 0) ∈ C ∞ ([0, L])). It is, of course, already suspicious that we have to demand this infinite smoothness condition for the initial state when the TDSE only contains second spatial derivatives, but let us ignore this issue for the moment and simply continue. To be specific we take the initial state to be [31] Ψ(x, 0) = exp which describes a wave packet localized in the interval [a − b, a + b] where we take a and b to be positive real numbers such that the wave-packet is properly localized within [0, L] (for instance we can take a = L/2 and b = L/100). One can check that this function is infinitely differentiable and that all derivatives are zero for |x − a| ≥ b.
Let us now see what we get if we insert this initial state into the formula of (21). For |x − a| ≥ b the formula then tells us that at any time Ψ(x, t) = 0 and therefore that the wave packet does not spread and never leaves the interval [a − b, a + b]. This is very much in disagreement with our intuition that free wave packets do spread. What has gone wrong? To understand this it is useful to talk about the exponent of a linear operatorÂ in a more abstract sense. Let us try to derive some conditions under which the definition makes sense. First of all when we act with with the exponential operator on a state Ψ we see thatÂ n Ψ must be well-defined for all n. This is certainly the case if the domain of the operatorÂ is the whole Hilbert space, i.e. D(Â) = H, since then any square integrable function is mapped to another square integrable function and we can then apply the operator repeatedly. Let us therefore assume thatÂ is defined on all of H. Then if we act with eÂ on a state Ψ then every term in the sum (23) is well-defined. However, this does not mean that the infinite sum converges. To guarantee this we must have that eÂΨ ∈ H or equivalently eÂΨ < ∞. We therefore want to make sense of the sum as a sufficient condition. For the right hand side to give a finite sum the terms Â n Ψ should not grow too fast with n. An estimate can be made for so-called bounded operators for which there exists a positive number M such that for all states Ψ in the Hilbert space. In particular, repeated use of this inequality implies that Â n Ψ ≤ M n Ψ . If we use this in (24) we have This means that (23) is well-defined for bounded operators. The problem with the HamiltonianĤ in our example is that it is not a bounded operator. Indeed a more careful analysis of the norms Ĥ n Ψ for our wave packet [31] shows that these grow so fast with n that the rightmost infinite sum in (24) diverges if we takeÂ = −itĤ and Ψ to be our initial wave packet (moreover there is also pointwise divergence for any |x − a| < b in (21) [31]). As a consequence (21) does not present the solution to our initial-value problem. Suppose now, however, that we discretize the TDSE of (20) on a finite spatial grid, i.e. we replace the differential operator by a finite matrix. In that case our Hilbert space is finite dimensional andĤ becomes a bounded operator and therefore the exponential exp(−itĤ) is well-defined by the series expansion. This immediately raises the question what happens when we make our grid spacing finer and finer and take a continuum limit. It will be instructive to do this calculation. The grid points x j = j∆x are labelled by an integer j. The second spatial derivative of Ψ(x, t) in grid point x j can be approximated by the finite difference formula where j = 1, . . . , N . We see that the determination of the second derivative of Ψ in N grid points requires the knowledge of the Ψ in N + 2 grid points. We use this feature to include the hard-wall boundary conditions. The first grid point is taken to be x 0 = 0 and the last will be x N +1 = (N + 1)∆x = L where we demand Ψ(x 0 , t) = Ψ(x N +1 , t) = 0 for all times, leaving N remaining points in between in which we have to determine Ψ(x j , t). Our Hilbert space will then be N -dimensional. By this discretization the Hamiltonian becomes an N × N matrix acting on the timedependent vector (Ψ(x 1 , t), . . . , Ψ(x N , t)) with the explicit form while the TDSE becomes an ordinary differential equation of matrix form In this case, since the Hamiltonian is now a bounded operator, we can take the exponential of a matrix and we find that It remains to calculate the action of H m on the initial state. To do this we expand the initial state in the eigenvectors of H. Since H is a symmetric matrix it has N real eigenvalues (l) and the eigenvectors ϕ (l) (x j ) are orthogonal. These are determined from the equation The eigenvectors turn out to be real as well for the case of symmetric matrices. If we therefore normalize the eigenvectors such that then we find that where l = 1, . . . , N . Let us now expand the initial state Ψ(x j , 0) in terms of the eigenstates of H. We have where due to the orthonormality condition (32) we easily find that the coefficients c l are given by With these preliminaries we can continue the evaluation of the infinite sum in (30). Inserting (35) into this equation we have If we insert the explicit form of the coefficient c l of (36) then we can write this as Or in other words Now that we have obtained an exact result for our discretized problem we can take the continuum limit. We let N → ∞ for a fixed value of L = (N + 1)∆x for the length of the box. This means that ∆x → 0 and the sum over k in (37) becomes a Riemann integral. We then get where We recover the well-known eigenfunctions and energies for the particle in a box. We can rewrite (39) as where formally the propagation kernel Strictly speaking the sum is not defined until after integration over x as in (39) (as at t = 0 it becomes the delta distribution δ(x − x )) and moreover it can depend in very complicated manner on the space and time arguments even for simple systems such as the particle in a box [32]. We can now define by (42) an evolution operator U (t) with the propertŷ We write this evolution operator by definition asÛ (t) = exp(−itĤ). So rather, than defining the evolution operator by a Taylor series we define it by a spectral representation involving the eigenfunctions and eigenvectors as in (42) [33,34]. It also follows from (39) that This means thatÛ (t) is a bounded operator which preserves the norm (in other words: it is unitary). This is an obvious requirement from a physical point of view as the total probability of finding a particle should be conserved in time. A more puzzling property is that thatÛ (t) is defined on any square integrable function, including non-differentiable functions on which the action of the HamiltonianĤ is not defined. Before we go into these issues let us now go back to our example of the localized wave-packet. If we now take Ψ(x, 0) as in (22) and apply (39) we will find that the wave-packet correctly spreads in the box (see Fig. 1). One can in fact prove that a wave packet enclosed in a bounded region of space at time t = 0 will have tails reaching all over space for almost all times t = 0 [35]. This is not difficult to understand from a physical point of view since the Fourier components of the initially localized wave function have momenta of arbitrarily high value allowing the particle to move arbitrarily fast. This also implies that the localized wave packet in an enclosed big box with hard walls will feel the presence of the boundary immediately. If we would have used other boundary conditions, such as the periodic one Ψ(0, t) = Ψ(L, t), dΨ(0, t)/dx = dΨ(L, t)/dx then immediately after t = 0 the time-evolution will be different. Clearly the formal series in (21) which is just specified by the differentiation rule has no information on such boundary conditions as they were neither encoded in the initial state nor in the exponential form of the time evolution operator. It is therefore no surprise that (21) can not be used to predict the time-evolution correctly. However, not all hope is lost in applying (21). Clearly we can apply (21) to a finite linear combination of eigenfunctions of the form since for this initial state we obtain from (21) that which is easily checked to be a valid solution of the TDSE of (20) with the right initial conditions. So why did the formal approach work in this case? We first note that in this case the function Ψ(x, t) is a real-analytic function of x and t, i.e., it has a Taylor series with non-zero convergence radius around any points x and t in its domain. Thus around any (x 0 , t 0 ) within its radius of convergence. In fact the function is so nice that it has infinite convergence radius and hence can be extended to the whole complex plane in x and t. Could it be that the formal expression of (21) would work for all real-analytic functions? Since the example of our localized wave packet of (22) is non-analytic at points x = a ± b (it has a Taylor series with convergence radius zero at these points) this would then explain in another way the failure of (21). One can prove that a function that is exactly zero on an interval of the real line can not be real analytic unless it is the zero function and therefore any initially localized wave packet fails to be real analytic. Let us give an example which shows that the requirement of analyticity is not sufficient to make (21) work. Rather than take the interval [0, L] we take the free propagation of a wave packet on the real line, i.e., our Hilbert space will be L 2 (R). The Schrödinger equation will again be given by (20) and as initial state we take a Lorentzian function This is a real analytic function on the whole of the real axis with a convergence radius of at least 1 for any Taylor expansion of Ψ(x, 0) in powers of x − x 0 around x 0 . If we insert this initial state into (21) we obtain the series From simple convergence criteria we see that this is a divergent series for any value of x and t. Therefore real analyticity is not a sufficient criterion to be able to apply (21).
To get a sufficient condition we follow the classical derivation given by Kowalevskaya [36]. Define the two time-dependent functions for which we will assume that they are real analytic. Those are new initial values if we exchange the meaning of x and t. This is necessary because [36] is concerned with initial-value problems where the time derivative appears in highest order. Then the general solution of (20) can be written as a formal power series as can be checked by insertion of this expression into (20). Let now the Taylor expansions of Ψ (0) (t) and Ψ (1) (t) be given by Then in terms of the coefficients c ν and c ν the expansion (51) at time t 0 attains the form where we defined Since we assumed Ψ (0) (t) and Ψ (1) (t) to be real analytic functions they have finite convergence radii R (0) and R (1) . Let R be a positive radius smaller than R (0) and R (1) . Then since both series (52) and (53) converge there exists a positive number g such that for all ν These conditions imply that for the coefficients b ν and b ν in (54) that One sees from standard convergence criteria that this implies that Ψ(x, t 0 ) as a power series in (x − x 0 ) has an infinite convergence radius. Now we understand what went wrong when we chose as initial state the Lorentzian function of (48). The function is real analytic but the radius of convergence is not infinite. The solution of the TDSE will then not be time-analytic, i.e., series of the form as in (52) and (53) do not exist for this initial state. Note that the requirement of infinite convergence radius by itself is not enough, one really needs to satisfy the constraints (57). For example the function has an infinite convergence radius but does not satisfy the constraints (57). If we now go back to the initial state of (46) we can check that it not only has an infinite convergence radius but that it also satisfies the constraints and therefore the direct exponentiation worked. Another physical example in which direct exponentiation is allowed is that of an initial Gaussian wave packet [37] since also for this case the conditions (57) are satisfied. The issue of time non-analyticity has been raised in some papers in connection with initial states that are not differentiable at cusps [38]. However, the analysis already carried out by Kowalewskaya shows that the situation is more severe. A Taylor expansion in time for Ψ(x, t) does not even exist for a large class of real-analytic initial states without cusps. We see that by putting rather stringent conditions on the initial state we can give the expansion in (21) a meaning. However, by means of the evolution operator defined in (42) we can give meaning to time-evolution for a much larger set of initial wave functions. Let us therefore forget again about analyticity and return to this more general definition of the evolution operator. In fact the expression is defined on any square integrable function, even on functions that are not in the domain of the Hamiltonian operator. Let us illustrate this with an example for the particle in a box again. We apply the evolution operator in (42) to the normalized initial state given by This initial state can be expanded in the eigenstates ϕ (l) with expansion coefficients which only gives a non-zero value when l is odd. The time-evolution is then given according to (39) as It turns out that this function for a given value of x is continuous but nowhere differentiable with respect to time. Moreover at almost all times it is continuous but nowhere differentiable as a function of x [39,40]. A snapshot of |Ψ(x, t)| 2 is displayed in Fig. 2. . This result seems puzzling at first sight, how can we get a nowhere differentiable function as a solution of a dynamics governed by a partial differential equation? One would expect that any solution would at least be once differentiable in time and twice with respect to the spatial coordinates. The problem is that the initial state is not in the domain D(Ĥ) of the Hamiltonianc. What one can show is that if the initial state is in the domain of the Hamiltonian then a well-behaved time-evolution, i.e., remaining in the domain, is guaranteed. However, since this domain is dense in L 2 (meaning that any element in L 2 can be approached to arbitrary accuracy with an element of D(Ĥ) as measured by the L 2 -norm) and the evolution operator is bounded, the domain of the evolution operator can be extended by continuity to all of L 2 , such that a time-evolution is well-defined for any square integrable initial state. Thus Schrödinger dynamics still exists in the so-called mild sense (see also Appendix B) but to define this properly the Schrödinger equation needs to be transformed to integral form A solution of this equation is called a mild solution (but later for time-dependent Hamiltonians we need to generalize the definition of mild solutions again so it will only appear here in this form). It turns out that when we integrate Ψ(x, s) of (61) between zero and t the resulting solution becomes twice differentiable with respect to x (see Fig. 2) and is in the domain of the Hamiltonian such that the last term in (62) is well-defined (note that integration and differentiation cannot be interchanged in these cases). This is, in fact, a general feature that can be proven using semi-group theory [41]. So our non-differentiable solution (61)  for all Ψ and Φ in the domain ofÂ. Let us give an example for the momentum operator for our particle in a box. We define the momentum operator bŷ and let its domain be the one times continuous differentiable functions Ψ on the interval [0, L] with boundary conditions Ψ(0) = Ψ(L) = 0 . We have by partial integration.
So clearlyp is a symmetric operator. Note, however, that we can extend this definition to larger domains. This is easily seen from our example. Since Ψ(0) = Ψ(L) = 0 we do not need to put any conditions on the functions Φ at the boundary to make the boundary term vanish. It will for example suffice that they are simple square integrable and differentiable to make (65) valid. The operatorp acting on Φ can therefore be regarded as an adjoint operator acting on a larger domain. Let us make this statement more precise with a definition. For a given operatorÂ with domain D(Â) we take D(Â † ) to be the set of all Ψ ∈ H for which there exists a χ ∈ H such that for all Φ ∈ D(Â). This defines the adjoint operatorÂ † Ψ = χ on D(Â † ). Clearly it follows from this definition that for symmetric operators D(Â) ⊆ D(Â † ), i.e., the domain ofÂ is equal to or a subset of the domain ofÂ † . In case D(Â) = D(Â † ) the operator is called self-adjoint. Such operators can be diagonalized (or more precisely have a spectral representation [33,34]) and their eigenvalues are real which are key features used in quantum mechanics. WhenÂ is a self-adjoint operator a unitary time-evolutionÛ (t) = exp(−iÂt) (defined by its spectral representation) can be defined on all of the Hilbert space (L 2 ([0, L]) for our example) and if Ψ(0) ∈ D(Â) the function Ψ(t) =Û (t)Ψ(0) ∈ D(Â) satisfies the evolution equation This important statement is also called Stone's theorem [33]. The issue of selfadjointness is thus very important for establishing the solvability of evolution equations such as the TDSE. Let us go back to our example of the momentum operator for the particle in the box. We already noted that the action ofp on Φ in (65) could be defined on a bigger domain. This means thatp is not self-adjoint. The adjoint operatorp † has the same appearance as a differentiation rule but the domain ofp † is larger, D(p) ⊂ D(p † ) and thereforep =p † . A quick calculation shows that the operator p that we defined does not have eigenfunctions satisfying zero boundary conditions, as one would expect ifp were self-adjoint. Furthermore one can see that the timeevolution equation withÂ =p has no solution either for any initial state (since the formal solution is Ψ(x, t) = Ψ 0 (x − t) and will be ill-defined as soon as it hits the wall of the box). So this is a simple example how self-adjointness is important for the solvability of evolution equations. It is clear from our example that if we put less restrictions on the functions in the domain of the momentum operatorp then we reduce the domain ofp † as we need more restrictions on Φ to make the boundary term in (65) vanish. We might therefore expect that for a sufficiently large extension of the domain ofp we have D(p) = D(p † ) and the momentum operator becomes self-adjoint d. The idea is to join the end points of the box and form a ring. The most general condition under which the boundary terms in (65) vanishes is [33] Ψ(0, t) = e iγ Ψ(L, t), where γ is a real number (note that this is indeed an extension of our earlier domain).
To make the boundary terms disappear these conditions must apply both to Ψ and Φ in (65), unless Ψ(0, t) = 0. If we the take as domain the space of functions such that both Ψ and dΨ/dx are square integrable (more technically called the Sobolev space H 1 ([0, L])) then with these boundary conditions the operatorp becomes selfadjointp =p † . For each γ we therefore find a different self-adjoint domain. The case γ = 0 corresponds to the standard periodic boundary conditions (i.e. the momentum operator of a quantum particle on a ring of circumference L). If we useÂ =p in the evolution equation (67) with this periodic domain (more precisely requiring Ψ and its first derivative to be square integrable and satisfying periodic boundary conditions) then Stone's theorem guarantees a solution for an initial state which is in the domain ofp (note that the self-adjoint "Hamiltonian"p is unbounded from below). Let us now go back to our original TDSE problem of (20) and discuss the self-adjoint domain ofĤ. A straightforward calculation shows that d We note that on a Hilbert space of finite dimension any symmetric operator is self-adjoint. If we discretize the momentum operator on a finite grid with zero boundary conditions then one finds that eigenfunctions exist but that they become ill-defined in the continuum limit while its eigenvalues diverge in the same limit [42].
It is clear that the Hamiltonian is symmetric on the domain C 2 ([0, L]) of twice continuously differentiable functions that satisfy Ψ(0, t) = Ψ(L, t) = 0. Since the derivatives dΨ/dx do not need to vanish at the boundary we see that also for Φ we need to require the boundary conditions Φ(0, t) = Φ(L, t) = 0 in order to make the boundary term in (69) vanish. It is still not clear, however, that the specified domain would makeĤ self-adjoint. Indeed the domain can be extended without violation of the boundary conditions. It is easily imagined that we could have a series of smooth functions ϕ n satisfying the boundary conditions which converges in the L 2 -norm to a function ϕ such that the seriesĤϕ n would converge in L 2 -norm to another function ξ. The function ξ need not be continuous but merely square integrable (in which case ϕ is not in C 2 ([0, L])). If this is the case then we can extend the domain of our original operator by definingĤϕ = ξ. Clearly for functions obtained by this limit procedure the inner products in (69) are finite since due to the Cauchy-Schwarz inequality we have | Ĥ Ψ|Φ | ≤ Ĥ Ψ Φ and bothĤΨ and Φ are square-integrable. Moreover, both inner products on both sides of the equality sign are identical since the series ϕ n preserves the boundary conditions. The domain of the Hamiltonian is by this procedure the space of functions Ψ for which the norm is finite and where Ψ satisfies the zero boundary conditions. This space of functions is equal to the Sobolev space where the superindex 2 refers to the second derivative and the subindex 0 to the zero boundary conditions [41,43]. One can now show thatĤ is self-adjoint on H 2 0 ([0, L]) [34]. Based on all these preliminaries we can finally conclude with the important statement that a solution of the TDSE exists with the property but still in L 2 then a trajectory Ψ(t) in Hilbert space exists but Ψ(t) will be outside this domain for all t. The solution then solves the mild version of the TDSE. Our example of a constant function as initial state can now also be understood better. This function can can be approached as close as one wants with C ∞ -functions ϕ n that satisfy zero boundary conditions (since they are dense in L 2 ). But to do so the first and second derivatives near the endpoints at 0 and L become very large in such a way that Ĥ ϕ n diverges. Our initial state therefore had infinite expectation value of energy and was outside H 2 0 and this lead to our strange non-differentiable solution. All these things discussed here can also be extended to the case of many-particles systems with self-adjoint and time-independent Hamiltonians. The main problems then is to know whether, for instance, the molecular Hamiltonian with Coulombic potentials actually is self-adjoint. This is answered by a famous theorem of Kato. For the case of time-dependent Hamiltonians the situation is much more complicated. First of all, due to a lack of time-translational invariance of the Hamiltonian the unitary evolution operatorsÛ depend not only on the length of the time-interval of propagation, but on the initial time t i and final time t f of the time-propagation, i.e. one writesÛ (t f , t i ). These issues are discussed in more detail in the following sections. The main aim of this section was to serve as an introduction to these section as they illustrate by a simple example the subtleties of the TDSE.

History of explicitly time-dependent initial-value problems
Before we start with the precise mathematical definition of the TDSE, let us first give a short historical overview of such explicitly time-dependent Cauchy problems. According to Kato [44] the first investigations of such initial-value problems with timedependent operators date back to Phillips [45] in 1953. In this work the perturbative expansion of the full evolution, i.e., the sum of all possible combinations of free propagation and the interaction with a time-dependent potential, is shown to converge under certain conditions. In contrast to this, one can also consecutively propagate along a direct path for short time intervals, assuming an evolution operator with a time-constant potential, then making those time-intervals smaller to get the desired evolution operator in the limit. This was achieved by Kato [46] and others, where the Hamiltonian was supposed to be maximally dissipative at all times, a property that self-adjoint operators automatically exhibit. However, conditions on the relations of the Hamiltonian at different times can exclude typical cases of external potentials. A more up-to-date summary of these techniques can be found in Pazy [47, ch. 5] where the notion of stable families of infinitesimal generators is used. The method of Kato [46] is later used in the books of Reed and Simon [48, Th. X.70], where the authors develop a theory specifically for the Schrödinger case for one particle, allowing also singular Coulombic potentials. This approach, which we follow in this review by and large, can be extended to the N -particle case with methods given in the same book.e Other approaches [49,50] use extended Hilbert spaces that involve time. But all these results apply to abstract operators, not beneficially taking into account any special structure of a Schrödinger Hamiltonian. On the other hand, Wüller [51] treats the special case of the Schrödinger Hamiltonian for a single quantum particle with certain moving Coulombic potentials. This specific approach unitarily transforms the equation to a static singular potential and then uses results from Tanabe [52] which in turn rely on Kato [46]. The study in Yajima [53] resumes the perturbative approach of Phillips [45] and combines it with specific properties of the Schrödinger Hamiltonian, e.g., Strichartz-type estimates to relate the freely evolved wave function with the initial state. The results apply to arbitrary spatial dimensions and thus allow to consider multiple particles in three-dimensional space. While the conditions on the temporal behaviour of the applied potentials are very mild (less restrictive than those employed in Reed and Simon [48]) the conditions on the spatial behaviour of the potentials become restrictive for more than one particle in three dimensions (Coulombic potentials are excluded for instance). Since we will use certain results of this approach, a brief introduction to [53] can be found in Appendix B. A slightly more general approach than in [53] was presented by D'Ancona et al. [54] just relying on a fixed-point procedure to show existence and uniqueness of Schrödinger dynamics rather than explicitly constructing the evolution operator as a Neumann series.

Classical and non-classical solutions to the initial-value problem
In this section we will give a general discussion of the initial-value problem, also known as the Cauchy problem, of the TDSE and introduce the notion of classical and non-classical solutions. In (1) we defined the TDSE for a many-electron system. The many-body wave function Ψ(x, t) is a function of both the space and the spin variables and therefore the inner product with another wave function Φ(x, t) is given by an summation over spin variables as well as an integration over spatial variables, i.e., where we used the notation that r = (r 1 , . . . , r N ) is the collection of spatial variables and that σ = (σ 1 , . . . , σ N ) is the collection of spin variables. For a discussion of the analytical properties of the wave function it is convenient to deal separately with the space and the spin variables. The function Ψ can be expanded as a finite expansion in terms of a product of space and spin functions. More precisely The requirement that the wave function is anti-symmetric under simultaneous interchange of space and spin variables implies that the functions Ψ j and θ j transform under an M -dimensional representation of the symmetric group S N of permutations of N elements upon permutation of the particle labels [55,56]. The spin functions θ j are typically chosen in such a way that the wave function is an eigenfunction of the total spin operators S z and S 2 and the choices of the eigenvalues determines the dimension of the representation M . Explicit techniques for construction of such spin functions are described in [55,56]. Since the spin functions θ j are linearly independent the initial value problem of (1) is equivalent to M separate initial-value problems Mathematically all these initial-value problems are equivalent, such that we can forget about the sub-index j again. The corresponding Hilbert space is simply the space of square integrable functions on R 3N , or in case we wish to discuss an N -particle system enclosed in a volume Ω ⊂ R 3 , the square-integrable functions on the configuration space Ω N . Mathematically this Hilbert space is denoted as H = L 2 (Ω N ) with inner product where the integration is over Ω N or R 3N depending on the situation of interest. This is the setting in which we will discuss the solvability and properties of the TDSE. When we discuss the properties of observables such as densities and currents we naturally have to remember to properly take into account the spin-structure of the many-particle wave functions.
Since the time-variable is a parameter in quantum mechanics and not a coordinate (i.e. there is no generic time operator) we will often suppress the spatial coordinates when writing the TDSE. The initial value problem (73) is then written as The HamiltonianĤ(t) is a self-adjoint operator on H parametrically dependent on times in the time interval [0, T ]. The time T is positive and can be arbitrarily large but we take it to be finite since it appears in estimates in Sec. 3.5. The solution will be a trajectory Ψ(t) in Hilbert space, or more mathematically a mapping from the interval [0, T ] to H. One of the most basic requirements that we can put on such a mapping is that it is continuous. Since the space H is a Hilbert space continuity is defined with respect to the Hilbert space norm. This means that Ψ(t) is defined to be a continuous mapping when the L 2 -distance of two "snapshots" of wave functions at time ∆t apart goes to zero as the time-difference goes to zero, or more precisely that The space of such continuous functions will be denoted as C 0 ([0, T ], H). In physical applications it is reasonable to further demand that the energy expectation value is finite. To guarantee this property we need more than continuity of the mapping. As we will discuss in more detail later this is guaranteed if also Hilbert space trajectories with this property are called continuously differentiable mappings and to indicate that Ψ(t) has this property we write Ψ ∈ C 1 ([0, T ], H) where the super-index 1 refers to the first-order derivative with respect to time. A solution of the TDSE which is such a C 1 -mapping is called a classical solution. We are, however, not always able to find a classical solution to the TDSE. This happens, for instance, when the initial state Ψ 0 is normalizable but whenĤΨ 0 is not, in which case the initial state is not in the domain of the Hamiltonian and the expectation value of the energy might be infinite. This was, for example, the case for the constant initial state of (59) for the example of the particle in a box discussed in Sec. 2.1. However, if we generalize the notion of a solution, even such problems can be solved uniquely [48,53]. Such generalizations of solutions to the TDSE can be found in different ways (see Appendix B for time-dependent Hamiltonians). The simplest way to do so is by first considering time-independent HamiltoniansĤ. In this case the Hamiltonian gives rise to a unitary evolution operatorÛ (t) = exp(−iĤt) defined by its spectral representation as discussed in Sec. 2.1. Although the Hamiltonian is not defined on all of H, its evolution operator is a bounded operator and can therefore be uniquely extended to all square-integrable wave functions. Consequently, we have a unique generalized solution to the TDSE Ψ(t) = exp(−iĤt)Ψ 0 even if the initial state Ψ 0 is not in the domain of the Hamiltonian. The mapping Ψ(t) regarded as a trajectory in Hilbert space is in that case not differentiable but still continuous, i.e. Ψ ∈ C 0 ([0, T ], H). Such a trajectory will be called a non-classical solution. It does not solve (75) but it might solve an equivalent but less stringent version of the TDSE, such as (62) which we discussed in Sec. 2.1. Since unitary evolution conserves the norm it follows in the case of time-independent Hamiltonians that ifĤΨ 0 is not normalizable then alsoĤΨ(t) is not normalizable and we therefore can differ between classical and non-classical solutions by their initial state. However, in the case of explicitly timedependent HamiltoniansĤ(t) we cannot straightforwardly use the same construction, since its spectral representation changes with time.
The discussion in this section has been very general as we did not discuss the actual properties of the Hamiltonian. This will be the topic of the next section in which we will address this issue in more detail and discuss for which class of external potentials and interactions the TDSE of the many-electron system is guaranteed to have a solution.

Existence and uniqueness of solutions to the Schrödinger equation
In this section we investigate under which conditions (and in which sense) we can define an evolution operator for a time-dependent Hamiltonian of a many-electron system. Before we discuss the existence of solutions of the TDSE we want to guarantee that the HamiltonianĤ(t) is a self-adjoint operator for all times t. This is not only a basic requirement for any observable as dictated by the mathematical structure of quantum mechanics, but we have also seen in Sec. 2.1 that self-adjointness is an important property which is closely connected to the solvability of evolution equations such as the TDSE. We first give conditions such that the time-independent HamiltonianĤ 0 is self-adjoint, then include also the time-dependent partV (t) and finally give conditions for the existence of an evolution operator.
The time-independent partĤ 0 of the Hamiltonian usually consists of two parts, the kinetic energy of the particlesT and the interaction between the particlesŴ . The kinetic-energy operator is given by the mapping Here the spatial derivative is meant in the weak sensef which is defined by integration against smooth test functions. The defining equation simply uses the basic equation of partial integration for the product of functions. For example, we say that a function Ψ defined on an open domain Ω ⊆ R n has the weak derivative Φ = ∇Ψ (n-dimensional gradient) when the following equation is valid for any function ϕ(r) which is infinitely differentiable and which is only nonzero on a bounded region (has compact support in mathematics language). The advantage of talking about weak derivatives is that one can talk about the derivatives of functions which do not have derivatives in the classical sense. For example, the weak derivative of the function Ψ(x) = |x| in one dimension is the equivalence class of functions Φ(x) which are equal to 1 for x > 0, equal to −1 for x < 0 and take an arbitrary value in x = 0. The concept of weak derivative considerably simplifies the mathematical treatment of partial differential equations. For an extensive discussion of these issues we refer to [43,41].
Let us now go back to (77). Obviously, not every square integrable wave function is again mapped to another square integrable function. Therefore, the kinetic energy operator is only defined for a restricted set of functions in the Hilbert spaceg, which is called its domain and which we will denote by D(T ). To be a proper domain this set of functions has to be dense in the whole Hilbert space, which means that we can approximate every normalizable wave function arbitrarily close (in the L 2 norm) with functions of the domain. Without this condition a unique adjoint operator can, for instance, not be defined. Now, there are two necessary conditions to make the kinetic-energy operator self-adjoint: the domain has to be equal to its adjoint domain andT has to be symmetric (see (66) and (63) for definitions). However, for particles restricted to a general volume Ω in three-dimensional space there are many domains that make the mapping of (77) a self-adjoint operator. We can construct those different self-adjoint domains by choosing different boundary conditions, such as periodic-or zero-boundary conditions. Depending on these conditions the properties of the associated operators can change dramatically. This is clear physically as, for instance, the Hamiltonian for a free particle in a box or for the free particle on a ring have different energy eigenvalues and eigenstates. Therefore, it is usually not enough to just prescribe the rule of an operator, such as the differentiation rule of (77), but one also needs to fix the domain and with it the boundary conditions. Only then we have the unique definition of a self-adjoint operator. The full three-dimensional space Ω = R 3 is an exception. In this case there is only one self-adjoint domain for the kinetic energy operator of (77) which is the Sobolev space D(T ) = H 2 (R 3N ) (see Appendix A for further details on Sobolev spaces). By defining the self-adjoint kinetic-energy operatorT , we have also chosen a specific set of functions that are guaranteed to have finite kinetic energyh. However, functions that have finite kinetic energy in one self-adjoint realization ofT might not have finite kinetic energy in another realization. For example the expectation value of the kinetic energy for the initial state (59) is infinite for hard wall boundary conditions but finite for periodic ones. After having discussed the kinetic energy operator we turn our attention to the twobody interactions and consider the static partĤ 0 =T +Ŵ of the Hamiltonian. The corresponding new rule for mapping wave functions is given by where two-body interaction w(r) is a real scalar function defined on R 3 which is typically taken to be Coulombic, i.e. w(r) = 1/|r|. In the following we will take it always to be a function of the inter-particle distance |r|. We want to ensure that the operatorĤ 0 is self-adjoint on the same set of (physical) wave functions as the kineticenergy operator, i.e., D(Ĥ 0 ) = D(T ). Using the theory of Kato perturbations [57] we can find rather simple conditions for this to hold i. For the case that we discuss particles in the whole three-dimensional space R 3 the operatorĤ 0 defines a selfadjoint operator with the same domain as the kinetic energy operator when w can be written as the sum w = w 1 + w 2 , one square integrable and the other bounded. This class of potentials is also known as the class of Kato perturbations. In a more mathematical notation we can write w 1 ∈ L 2 (R 3 ) and w 2 ∈ L ∞ (R 3 ). The space of functions L ∞ (R 3 ) is the set of functions w(r) for which there is a positive number M such that |w(r)| < M for all r (technically speaking this has to hold almost everywhere, meaning up to a set of measure zero). The space L ∞ has a norm but no inner product and is therefore not a Hilbert space, instead it is called a Banach space (see Appendix A for a further discussion). The class of Kato potentials on R 3 is written as and is again a Banach spacej. In the case h The most general set of finite kinetic-energy states is usually bigger, since it only needs to ensure that the expectation value is finite.
i In this context a Kato perturbation of a self-adjoint operatorÂ is a symmetric operatorB with From this it is obvious that every bounded operator, e.g., a multiplication with a bounded interaction potential w, is automatically a Kato perturbation. j This Banach space has norm w K = inf{ w 1 2 + w 2 ∞|w = w 1 + w 2 , w 1 ∈ L 2 (R 3 ), w 2 ∈ L ∞ (R 3 )} in which · 2 and · ∞ are the norms on L 2 and L ∞ respectively.
that we discuss particles restricted to a finite volume Ω we just have K(Ω) = L 2 (Ω). An important potential which is included in the class of Kato potentials K(R 3 ) is the Coulomb potential since it can be written as where θ(x) = 1 for x > 0 and is zero otherwise. The first term after the equal sign is square integrable and the second term is bounded.
In the final step we now add an explicitly time-dependent external potential to the time-independent HamiltonianĤ 0 to build the full HamiltonianĤ(t) = H 0 +V (t). Again applying the theory of Kato perturbations we find that the HamiltonianĤ(t) is a self-adjoint operator on the domain D(T ) for all times whenever the potential v(r, t) belongs to the class of Kato potentials. We note that important physical models such as the harmonic oscillator or the dipole fields are not included in the Kato class if we consider the full three-dimensional space R 3 . SinceT and H(t) have the same domain it also follows thatĤ(t)Ψ is normalizable wheneverT Ψ is normalizable. For a Hamiltonian in which the external potentials and two-body interactions are in the Kato class the total energy expectation value is therefore finite whenever the kinetic energy expectation value is finite. Now that we have identified the class of external potentials and two-body interactions for which the Hamiltonian is self-adjoint on the domain of the kineticenergy operator we can start to discuss the solvability of the initial-value problem for the TDSE.
There are now several ways of investigating the existence and uniqueness of solutions to the TDSE (see Sec. 2.2). For the time being we restrict ourselves to an approach similar to the one presented in [48]. Certain details and a comparison to a different approach based on purpose-build Banach spaces of potentials and wave functions [53] are discussed in Appendix B. From the previous considerations we have seen that if we take w(r) and v(r, t) to be Kato perturbations, the resulting HamiltonianĤ(t) has the same domain D(T ) at every time. The task is now to prove that a well-defined time-evolution exists. The idea of the proof is to divide the time- . In each such time interval the HamiltonianĤ(t i ) defines a self-adjoint operator and we know (by Stone's theorem already employed in Sec. 2.1) that a well-defined evolution operator exists for s, t ∈ [t i , t i+1 ]. For arbitrary times t and s in [0, T ] we can define the evolution operatorÛ k (t, s) by glueing together the evolution operators in different where in the product the operator with the latest time is always ordered to the left. It can then be shown [48] that .k Let us elaborate on this condition. If we view the potential as a trajectory v(t) in the space of Kato perturbations then (83) is valid when v(t) is a continuously-differentiable mapping with respect to the norm of K. The evolution operatorÛ (t, 0) is then unitary and therefore defines a unique continuous trajectory in Hilbert space for any normalizable initial state Ψ 0 , or more precisely Ψ ∈ C 0 ([0, T ], H). Therefore we have existence and uniqueness of a (generalized) solution for an important class of time-dependent potentials. Such potentials include for example molecular potentials of the form of (4) provided v ext (t) is a continuous differentiable mapping to the Kato-class. We note that the differentiability condition with respect to time on the external potential excludes a sudden switch-on, but by using the technique presented in [53] we can show existence and uniqueness of a generalized solution also for such situations (see Appendix B for more details). Further, if the initial state Ψ 0 is in the domain of the kinetic-energy operator D(T ) then also Ψ(t) ∈ D(T ) for every time t ∈ [0, T ] and therefore the expectation value of the kinetic and total energy are finite. In that case one can has Ψ ∈ C 1 ([0, T ], H) which is therefore a classical solution to the TDSE. Let us summarize the most important results of this section. We can establish a well-defined time-evolution if the potential is continuously differentiable in time with respect to the norm of K. We define this set of allowed potentials as Now we have properly defined (including the domains) the mapping from potentials to wave functions that we discussed in Sec. 1.2. For a given potential v ∈ V we can solve the TDSE for a normalizable initial state. There are now two cases to consider. Either the initial state is in the domain of the Hamiltonian or it is not. In the latter case the time-evolution of the initial state defines a continuous trajectory Ψ(t) and the trajectory regarded as functional of v is given by a map In the case the initial state is in the domain of the Hamiltonian the time-evolution of the initial state defines a continuous differentiable trajectory Ψ(t) and the trajectory regarded as functional of v is given by a map We have therefore established well-defined potential to wave function mappings. Clearly the mapping of (86) is the most relevant for physical applications as it guarantees the finiteness of the expectation values of the kinetic energy and the oneand two-body interactions.
k Considering the proof of [48, Th. X.70] it seems possible that Lipschitz-continuity in time with respect to the norm of K is enough.

From the wave function to observable quantities
After having established a well-defined mapping from potentials to wave functions we can continue with the discussion of the mapping from wave functions to physical quantities such as densities and currents. To calculate these quantities properly we have to take into account the correct spin structure of the wave function of (72) and use the inner product of (71). The solvability properties of the TDSE discussed in the previous section for each of the spatial parts of the wave function immediately imply the same solvability properties of the TDSE for full anti-symmetric space-spin function. Let us now consider the calculation of an arbitrary physical observable. If such an observable is described by a (time-independent) self-adjoint operatorÔ then we want to evaluate This expectation value is well-defined if the domain ofÔ contains the domain of the Hamiltonian. We already established that for the potentials V given by (84) the kinetic energy as well as the two-body interaction energy are finite. We therefore have well-defined functionals defined on the set of potentials in V. However, not all physical quantities are defined by self-adjoint operators on a Hilbert space. The most important ones for us are the density and the current density. The density was already defined in (6). If the initial state is further in the domain of the kinetic-energy operator like in the case of classical solutions (86) then the time-derivative of the wave function is well-defined and the density obeys the continuity equation where Here we used the same notational convention as in (6) and ∇ is the gradient corresponding to the (non-integrated) coordinate r. Depending on whether we wish to consider particles in the whole space R 3 or in a finite volume Ω the spatial integrations in this expression are restricted to R 3(N − 1) or Ω N −1 . In the following we will just use Ω for the volume with the understanding that possibly Ω = R 3 . We shall mention it explicitly whenever the distinction is relevant. The density n(r, t) has the obvious property that it is positive and that its integral is given by the number N of electrons. It therefore belongs to the space of L 1 (Ω) functions (see Appendix A for a definition). The density n(t) can thus be regarded as a continuous differentiable trajectory in this space. The potential to density mapping is therefore given by In the following we then want to investigate under which conditions and restrictions such a mapping between v and n is invertible. Obviously L 1 (Ω) is too general for the space of densities since it contains also negative functions and moreover the finiteness of the kinetic energy implies that ∇ n(r, t) exists and is square integrable [19]. The fact that n and v have the same degrees of freedom does at least give some hope that an inverse map may exist. To investigate this in more detail we start by considering an equation that connects n and v more directly. Formally, such an equation can be derived by combining the above continuity equation (90) with the so-called local-force equation of quantum mechanics, i.e., the time-derivative of (91) where the components of the vector Q are given by and ∂ k ≡ ∂/∂r k as well as summation over multiple indices is implied. The momentum-stress tensor is defined by (suppressing the dependence of the wave function on the different variables and again using the notational convention of (6)) The interaction-force density is defined by The combination of these two equations formally leads to the fundamental equation of TDDFT where It is obvious that this equation cannot hold for every possible Ψ[v] ∈ C 0 ([0, T ], H), since already the continuity equation might not be well-defined for generalized solutions like in (85). We therefore need extra conditions on the potentials and initial states that guarantee that the trajectory n(t) in L 1 (Ω) is twice differentiable with respect to time, i.e. more precisely n ∈ C 2 ([0, T ], L 1 (Ω)). To find those we analyse each constituent of (97) in detail. First, for a classical solution of the TDSE we know that T kl is at least in L 1 (Ω) since Ψ(t) is in the domain of the kinetic energy operator and thus twice differentiable. To guarantee that then the kinetic part of the operator q[v] is integrable, we restrict ourselves to initial states that obeyT 2 Ψ 0 ∈ H too and to potentials that stabilize this condition, i.e.,T 2 Ψ(t) ∈ H for all times t ∈ [0, T ]. This holds, for instance, if we impose periodic boundary conditions on the kinetic-energy operator and restrict ourselves to infinitely-often differentiable (in space and time) interactions and external potentials with the same boundary conditions [58]. Since in the given reference the question of stabilising an arbitrary number of derivatives is considered, we expect that for our case weaker conditions are sufficient. Based on the conditions for the stability of the domain D(T ) under time-evolution in the proofs of [48,53] we conjecture that it is enough that v, ∇v and ∇ 2 v are in C 1 ([0, T ], K(Ω)).
We point out, that physically such conditions are quite reasonable. If we, for instance, assume that the external potential is due to a L 2 (Ω) charge distribution and hence determined from the Poisson equation (see for instance (217) where Coulomb gauge on R 3 is employed), then v(t), ∇v(t) and ∇ 2 v(t) are in L 2 (Ω) ⊂ K(Ω). For the case Ω = R 3 such a potential is of the form v(r, t) = − dr ρ(r , t) |r − r | where ρ(r, t) is a square-integrable time-differentiable charge distribution (which excludes the case of external point charges which are described by delta distributions).
For example, such potentials arise when in molecules when the atomic nuclei are represented by finite charge distributions (for a more extensive discussion see [59]). For the interaction w it suffices that it is in the Kato class of potentials, since the derivatives with respect to r in the definition of W k [v] can be expressed as derivatives with respect to r 2 and thus by partial integration be absorbed by the wave function.
For the interactions the important case of pure Coulomb potentials is therefore allowed.
These conditions make q[v] integrable but note that since we rely on Hilbert space techniques in the later discussion on the inversion problem in Sec. 3.3 a L 2 (Ω) condition arises there. Finally we consider under which conditions the external-force term ∇ · (n[v]∇v) is at least integrable as well. By the product rule we can express the external-force expression in two terms ∇n[v] · ∇v and n[v]∇ 2 v. Under the above assumptions of ∇v(t) and ∇ 2 v(t) being in K(Ω) it holds that they are individually integrable l.
We will in the following denote such a set of potentials v for which an initial state with the propertyT 2 Ψ 0 ∈ H impliesT 2 Ψ([v], t) ∈ H by V. Then we have a mapping where N is the set of densities n[v] generated by all possible potentials in V. These conditions guarantee that the internal and external forces (as well as their divergences) are finite. On the basis of these domains we can discuss the bijectivity of the mapping v → n.
3. The density-potential mapping

Exemplification
So far we have discussed the potential-density mapping v → n, where v ∈ V is the set of potentials for which we have a unique (possibly generalized) solution of the TDSE. However, in order to invert this mapping we need to also show injectivity, i.e., that every density n has at most one potential associated. Then and only then the potential-density mapping is bijective which allows to define its inverse, the densitypotential mapping n → v.
As pointed out in the previous section, we will employ the fundamental equation of TDDFT (97) to establish injectivity and thus bijectivity. However, in order to do so we need to restrict the set of allowed potentials to the smaller set V ⊂ V that guarantees that the individual terms of the fundamental equation are all well-defined.
l For instance, for Kato perturbations ∇ 2 v ofT it holds that for Ψ in the self-adjoint domain But this does not imply that we could not have a bijective mapping (and can thus invert the mapping) for a more general set of potentials. For instance, if we restrict to all those potentials that guarantee a classical solution of the TDSE we can for a one-dimensional noninteracting spin-singlet problem with periodic boundary conditions construct the density-potential mapping explicitly [60].
The example we want to invert is where the θ is the anti-symmetric singlet spin-function with the explicit form So for this particular two-electron case the expansion in spin functions of (72) has only one term. For the spatial part of the wave function we impose periodic boundary conditions on the interval Ω = [0, L]. This system is routinely investigated in TDDFT due to its simplicity. The TDSE for two (noninteracting) particles can then be rewritten into the single-orbital TDSE with an initial state in the domain of the self-adjoint kinetic-energy operator with periodic boundary conditions that thus obeys ϕ 0 (0) = ϕ 0 (L) and ∂ x ϕ 0 (x)| x=0 = ∂ x ϕ 0 (x)| x=L . We can rewrite the initial state in its unique polar representation provided n 0 (x) > 0 everywhere for a unique phase S 0 . Therefore the density n 0 (x) obeys the above periodicity conditions and where m ∈ Z. Now, every (classical) solution ϕ([v s ], x, t) gives rise to a temporally continuously-differentiable density n([v s ], x, t) = 2|ϕ([v s ], x, t)| 2 that obeys the same boundary conditions as the initial density. Further, the resulting phases S([v s ], x, t) obey the same boundary conditions as the initial phase in (106). A sudden jump from m to m is not allowed, since it needs an external potential that is proportional to a delta-distribution in time, which is not part of V. Under these conditions we can explicitly invert the potential-density mapping v s → n (at least for some finite time-interval [0, T ]). The construction of the density-potential mapping in this case is now based on the continuity equation, which in the above polar representation reads as We can interpret the continuity equation as a Sturm-Liouville equation for the phase S for a given time-dependent density n supplemented with the boundary conditions of (106) for a fixed value of m. The unique solution for a periodic density n(x, t) > 0 therefore becomes [60,61] S([m, n], x, t) = L 0 dy K t (x, y)∂ t n(y, t) , with θ the Heaviside function and η(xt) = 1 2 x 0 dy n(y, t) + x L dy n(y, t) .
Consequently, we have infinitely many orbitals ϕ[m, n] = n/2 exp(−iS[m, n]) which are labelled by m and correspond to realizations of the same time-dependent density n starting from different initial states. If we interpret the periodic system as a quantum ring of length L, the different realizations correspond to different rotations of the ring [60]. Finally we can invert (104) for v and employ the continuity equation, which allows us to express the potential as a functional of the initial state and the density v s ([m, n], x, t) = 1 2 Thus, we have inverted the ususal map from potentials to densities and explicitly constructed a density-potential map n → v s .
Let us consider a few consequences of the existence of a density-potential mapping with the help of this explicit example. First of all, such a mapping implies that the wave functions and thus also all observables are functionals of the initial state and the density by the composite mapping n → v → Ψ. This is one of the main implications of the famous Runge-Gross result [23] and the very foundation of TDDFT. It allows us to determine any physical quantity by only knowing the density and the initial state (at least in principle). In our example we can give an explicit realization of the Runge-Gross result since by the above construction we can express the wave function Φ[m, n] as a functional of the initial state (labelled by m) and the time-dependent density n. Any observableÔ inherits the functional dependence by For instance the kinetic energy as a functional of the initial state and the density reads as where the first term is the Weizsäcker energy functional and the second term is an initial-state dependent velocity contribution.
A further detail of the density-potential mapping is the initial-state dependence. For every different possible initial state we have a different density-potential mapping and consequently different wave functions that generate the same density in time. This makes the construction of universal time-dependent density-functionals a lot more complicated, since in principle we would need to incorporate the initial-state dependence as well. Therefore, one usually employs ground-state DFT to get rid of the initial-state dependence [22,21,62,63]. By the Hohenberg-Kohn theorem, for any ground-state density there is (usually) a unique ground-state wave function that minimizes the Hohenberg-Kohn functional, i.e., the kinetic and interaction energy. Thus by restricting to ground-states as the only allowed initial states one can ignore the initial state dependence and have a "pure" density-functional. This restriction excludes initial densities that have nodes or densities with ∂ t n(x, t)| t=0 = 0. As we can then also see in our explicit example, the Hohenberg-Kohn functional corresponding to (112), i.e., has a unique minimum for m = 0 which singles out the ground state corresponding to the chosen density n 0 (x). While from a purely formal point of view it seems desirable that we "only" need to approximate the dependence of functionals on the density, the initial-state dependence can also be an advantage in practice, e.g., in the case of charge-transfer problems, which are briefly discussed later in this section.
The main approach to perform practical TDDFT calculations is the timedependent KS scheme discussed in Sec. 1.3. By employing two different densitypotential mappings we can determine the time-dependent density of a quantum system by solving an auxiliary non-linear problem. While usually this is done to determine the density of an interacting many-particle problem by a non-interacting auxiliary system, the KS construction allows to connect any two different systems. We can, for instance, connect two different non-interacting systems with two different initial states. In that case we will still have an xc potential but this then, in the absence of two-particle interactions, purely generated by initial state dependence. In our example at hand we can determine the density found by solving (101) for a fixed v s (x, t) starting from an initial state characterized by ϕ 0 = n 0 /2 exp(−iS 0 [m, n]), by solving a non-linear auxiliary problem of the form with a different initial state ϕ 0 = n 0 /2 exp(−iS 0 [m , n]). The xc potential in this casem is determined by m We note that since we look at two non-interacting problems the Hartree term is zero by construction.
Note, that with help of the current and the continuity equation we can express the terms ∂ t n time-locally by the orbital ϕ. Thus we do not need any further information to uniquely solve the KS equation than the chosen potential of the original problem v s (x, t) and the initial states. If the exchange-correlation potential would depend on higher-order derivatives, we might need further information to determine the unique solution.
Even though the above density-potential mapping v s [m, n] and the example of an exchange-correlation potential v xc [m, m , n] do only depend on the instantaneous density (and its time-derivatives), this time-local behaviour is not a generic feature. Actually, the density potential-mappings usually depend not only on the initial state but also on the density at previous times. This property is termed memory [62,63,22,21] and for fixed initial state Ψ 0 it is formally expressed by where the inverse linear-response kernel χ −1 [n] is assumed to obey Here  ([n], r, r , t − t ) [22,21]. This form also shows most clearly why memory is a necessity of most density-potential mappings, especially in the context of the KS construction. If we Laplace-transform (related to the Fourier-transform with a step-function) the (inverse) linear-response kernel from t − t to the frequency ω we find that the linear-response kernel has poles at the eigenfrequencies of the (timeindependent) system that has the initial density n 0 as its ground-state density. The inverse linear-response function has zeros at these frequencies. If we then want to simulate the linear-response of an interacting system (starting from its ground state) by the linear-response of a KS system [22,21], we see that the (linear-response of the) exchange-correlation potential needs to cancel the poles of the auxiliary system and generate the poles of the original interacting system (see [60] for an analytic example of these properties and [64] for a numerical reconstruction). Thus, for any KS construction where we want to simulate an interacting by a non-interacting system, the exchange-correlation potential necessarily will have memory (since the spectrum of their respective Hamiltonians are very different). This argument also illustrates, why we get away with such a simple (time-local) exchange-correlation potential in our example above: we simulate one system by another system with a similar Hamiltonian (and thus with a similar spectrum). This result together with the knowledge that memory and initial-state dependence are closely related [62], i.e., it is possible to replace v([Ψ 0 , n], x, t) by v([Ψ(t ), n], x, t) for all 0 < t < t using the functional relation Ψ(t ) = Ψ([Ψ 0 , v], t ), can be used to simplify the KS scheme and make it more reliable. Firstly, by choosing a KS system that already has some of the properties of the interacting many-body system one wants to simulate, the complexity of the exchange-correlation potential can obviously be reduced. The more similar the original system and the KS system are, the less the density-potential mappings will differ and the exchange-correlation potential will become less signficant. A further simplification is possible, if one employs the initialstate dependence and chooses an initial KS state that incorporates some of the physical properties of the interacting system [60]. For instance, one of the major challenges the current (usually time-local) approximations to the exchange-correlation potential face, is the proper description of charge-transfer reactions [65,66,67]. If the simulations are started from the ground state it is well-known [68,69,70,71,72] that one needs time-nonlocal (frequency-dependent) exchange-correlation functionals to get the transfer process right. However, recent results [73,74] show that if one starts from an excited state, the usual simple (time-local) approximations reproduce a charge-transfer reaction reasonably well.
Finally, let us make use of our analytical expression v s [m, n] from (110) and give an example of the exchange-correlation potential for two different values of m. Since we have been briefly discussing the challenge of charge-transfer problems within TDDFT, we will consider a very simple toy model of such a process. The interacting reference system is a two-electron system on a ring Ω = [0, 12] which has the periodic interaction w(x) = cos(2πx/12)/2 and starts in the ground state of the time-independent potential which is displayed in red in Fig. 3. The first two terms in this expression describe wells around x = 4 and x = 8 while the last term reduces the depth of the well around x = 8. In the same figure we also have given the ground-state density n 0 which is localized in the left well. Now, we prescribe a time-dependent density profile evolving from the ground state density n 0 in which we split the two-particle density into two parts and move half of the density to the right well in T = 20 as is displayed in Fig. 4.   To determine the exchange-correlation potential for these two different initial states we finally also need to know the external potential v[Ψ 0 , n] of the interacting reference system with interaction w that generates the same density via time propagation (see (15)) starting from the interacting ground state Ψ 0 . While we do not have an analytical formula in this more complex situation (except for the initial time when it is equal to v 0 ) we can determine this potential from the numerical procedure that will be discussed in Sec. 5. The results are then displayed in Fig. 6.  Here we see the external time-dependent potential v ext (where we have subtracted the time-independent part, i.e., v ext = v[Ψ 0 , n] − v 0 ) that forces the interacting twoparticle system to obey the above prescribed rigid charge transfer. The Hartree potential v H is the same for both initial states, since it only depends on the instantaneous density. All the memory and initial-state dependence is found in the exchange-correlation potentials given by Obviously simple time-local approximations to the exchange-correlation potential cannot capture the rich structure of v xc in this case.

The local-force equation approach to the density-potential mapping
Let us now consider the general case of a density-potential mapping. In Sec. 2.5 we have seen that the local-force equation makes a connection between the density of a many-electron system and the potential v that generates the density n by timepropagation of the TDSE. Here we will outline how this equation can be used to show invertibility of the mapping v → n for a given initial state, i.e., the existence of the density-potential mapping, and how it can be used in an iterative way to calculate the potential v that generates a given density n.
We first rewrite (97) as Here both n([v], r, t) and q([v], r, t) are functionals of the potentials. Suppose now, however, that we fix the density n(r, t). Then we have a non-linear equation for the potential v, i.e., If we have prescribed a density that we have generated by time propagation of the initial state Ψ 0 with an external potential v, we can use this equation to ask whether a density-potential mapping exists. Namely, if we can show that the only potential that solves this equation is v, we can show that the mapping v → n is injective and hence invertible. To answer whether v is the only solution, we linearize the above nonlinear equation by an iterative procedure On the other hand, if we prescribe a density n for which we do not know a priori that it is generated by solving the TDSE we can try to employ (123) to find an appropriate v. So we start by making a guess for an initial potential v 0 (r, t) on a time interval [0, T ]. Then we propagate the TDSE with this potential and a given initial state (compatible with the initial density) to calculate q[v 0 ]. Then we can calculate a new potential v 1 (r, t) by solving (123). With v 1 we can repeat the procedure to find a new potential v 2 , etc. In this way we have constructed the mapping which defines a series {v k } of potentials. The goal is now to show that under certain assumptions v k → v (in some norm sense) for k → ∞, i.e., that v is a fixed point of (123). However, if we can invert −∇ · (n∇) then at t = 0 we are already converged after the first iteration, i.e., v 1 (r, 0) = v(r, 0) since q([v 0 ], r, 0) is given in terms of the initial state only and therefore independent of v 0 . Under which conditions we can invert this operator will be discussed in detail in Sec. 3.3. Assuming this for the moment we can conclude that for small enough times already v 1 will be close to the exact v. Therefore it seems likely that the next step in the iteration will be even closer to v unless the q[v 1 ] is very different from the exact q[v]. This, however, can only happen if the internal-force densities generated by two potentials that are arbitrarily close differ strongly. Such a situation would be quite unphysical, since arbitrarily small changes in a potential would lead to totally different dynamics within very short times (making any prediction for real system impossible). Although one could imagine such situations (for instance in the case of non-classical solutions to the TDSE discussed in Sec. 2.3), we exclude them from our considerations and will discuss this in more detail in Sec. 3.5. Summarising we therefore conclude that for short enough time intervals the iteration scheme based on (123) is expected to converge fast, which is an important feature that is also used in our numerical implementation of Sec. 5. However, even if we assume that v k → v, does this guarantee that v really generates the prescribed density by propagation of the initial state Ψ 0 ? The potential that we found is by construction a fixed point of (123) but does also obey its own local-force equation (97). If we subtract both equations and denote ρ = n[v] − n the difference between the density generated via propagation of v and the prescribed density n, we find [75,61] Now, if we assume that the prescribed density obeys the minimal restrictions (here the current operator is defined byĵ ) the above equation is a linear evolution equation with initial conditions ρ(r, 0) = ∂ t ρ(r, t)| t=0 = 0. Thus the question whether v generates the prescribed density via propagation reduces to the question whether (125) has only ρ = 0 as unique solution for the above initial conditions. While for analytic potentials v one can rigorously show by the classical considerations of Kowalevskaya [36,43] that this is true, we are not aware of a general proof for this statement. However, due to the fact that we can impose further conditions on ρ [75,61], e.g., dr ρ(r, t) = 0 for all times, we presume in the following that it holds true as it indeed seems highly probable. Under this assumption our iteration scheme, if it converges, does indeed reproduce the prescribed n.
Finally we note that the above iteration procedure in terms of v k also gives rise to iterations in Ψ k = Ψ[v k ]. Thus one could make the TDSE part in the above iteration explicit by introducing for the initial state Ψ 0 , where the individual potentials inV ([Ψ k ], t) are determined by (123). If we converge then Ψ k → Ψ and the the resulting wave function solves the non-linear TDSE where accordingly the v[Ψ] are determined by (122). We therefore see that the fixed-point procedure is equivalent to the nonlinear-TDSE approach to TDDFT [76,77,63,78].

The Sturm-Liouville operator and its invertibility
We see from (122) that, when we denote the right-hand side by ζ(r, t), that we need to solve an equation for v of the form − ∇ · [n(r, t)∇v(r, t)] = ζ(r, t) for a given inhomogeneity ζ with the property Ω dr ζ(r, t) = 0 The fact that the inhomogeneity integrates to zero is due to the fact that q is a divergence and that the total number of particles is conserved. From the previous considerations we have seen that the invertibility of this Sturm-Liouville equation (130) (usually this terminology is only used in the one-dimensional case but we will employ it for also for higher dimensionality) is fundamental to the construction of a densitypotential mapping. This equation appears in the iterative sequence of potentials for a given density n of (123) as well as in the definition of the (equivalent) non-linear TDSE approach (which itself gives rise to a iterative sequence of wave functions). We therefore provide in this subsection a detailed discussion about the properties of the Sturm-Liouville operator −∇ · [n(r, t)∇] and conditions for the invertibility of the Sturm-Liouville (130). The discussion will be general in the sense that we will not use the explicit form of the inhomogeneity ζ in terms of the density and the divergence of the local forces. Possibly stronger results may be obtained by taking into account this specific structure but we will leave this for future work. For simplicity we start with the one-dimensional version of (130) To discuss this equation in a general setting we will in the following regard the Sturm-Liouville operator as a linear operator in the Hilbert space of square integrable functions on either some finite interval [a, b] or one whole real line. In this setting it is then important that both v and ζ are in the Hilbert space and therefore squareintegrable. The issue of a unique inversion up to a pure gauge is then equivalent to the question of the (possible) self-adjoint domains which have the purely time-dependent function as the unique square integrable eigenfunction with zero eigenvalue. In the following we therefore want to say something about the eigenvalues and eigenfunctions of the Sturm-Liouville operator (and whether they exist) for certain types of boundary conditions. Let us start by a few simple manipulations to discover some general features. By integration of (132) we find Clearly we can find an equation for ∂ x v(x, t) provided we are allowed to divide by the density. This is allowed everywhere except at the boundaries of our interval where the density may go to zero. Let us, however, assume that we are allowed to divide by the density. Then we can do another integration to write v(x, t) = v(a, t) + n(a, t)∂ x v(a, t) In case that ζ = 0 we see that the general form of the solution is v(x, t) = c(t) + d(t) x a dy 1 n(y, t) .
This corresponds to the eigenfunction φ 0 (x, t) of the Sturm-Liouville operator with zero eigenvalue. We see that the zero eigenfunction is more general than just a pure gauge v(x, t) = c(t). We can always add φ 0 (x, t) to a particular solution of (132) and it will be another solution. To make the inversion unique up a gauge we therefore have to impose boundary conditions such that the second term in (135) vanishes. It turns out that which boundary conditions we can choose and which eigenspectrum we can obtain depends very much on the behavior of the function x a dy/n(y, t) where a is a boundary point. The simplest case is when n(x, t) ≥ ε > 0 for some positive number ε. In this case division by the density is no problem and the mathematics is the simplest. The most relevant physical case in which this happens is the case of a system with periodic boundary conditions, such as is the case for particles on a ring. In this case it is then natural to also impose periodic boundary conditions on the solutions of the Sturm-Liouville equation, and a quick calculation then shows that the only possibility for the zero eigenfunction in (135) to be periodic is to demand that d(t) = 0. In this case we can choose periodic boundary conditions also on the potentials, making the Sturm-Liouville operator self-adjoint with a purely discrete spectrum [79,80,61], i.e., Since ζ by assumption (131) is perpendicular to φ 0 (x, t) = c(t) the (pseudo-) inverse operator is well-defined and bounded, i.e., where D = λ −2 1 . Alternatively, one can show boundedness also from (109) for m = 0 (see also [61]), since is the spectral form of the Green's function of (109)n. This shows that for every ζ which is perpendicular to φ 0 = d we have a well-defined v, which is obtained by the action of the inverse operator in (139) on ζ. This brings us to the more difficult case in which the density can become zero at one of the boundary points. To be more definite we take zero boundary conditions on the interval [a, b] where n(a, t) = n(b, t) = 0. We consider a class of densities such that n(x, t) ∼ (x − a) 2 close to the boundary (and similarly in point b) which is the most generic case for particles in a box. For this case the integral x a dy/n(y, t) diverges and we have to use so-called singular Sturm-Liouville theory [80]. Within this theory the boundary points of our problem are so-called limit-point endpoints of the Sturm-Liouville equation [79,80,61]. In this case we only have one possible self-adjoint domain (for details we refer to [80]). Any twice-differentiable potential v with a behaviour at the boundaries which is less singular than v(x) ∼ (x − a) −1/2 will be in this domaino. The self-adjoint Sturm-Liouville operator then (possibly) has also n To be precise, since for Ω bounded the spectral form of the inverse is defined on L 2 (Ω) ⊂ L 1 (Ω) we can define it as the restriction of the Green's function of (109) onto L 2 (Ω). [−∂x(n∂x)] −1 is a bounded operator on L 2 (Ω) as well as on L 1 (Ω). We further note, that the specific number of the bound D might depend on whether we consider the operator on L 2 (Ω) or on L 1 (Ω). o We point out, that although from the condition on the differentiability of n∂xv potentials as singular as v(x) ∼ (x − a) −1 are possible, they are no longer in L 2 (Ω) and thus outside of the self-adjoint domain.
a continuum in the spectrum and its spectral representation is In our case of limit-point endpoints the unique zero eigenfunction is given by φ 0 (x, t) = c(t) and since n(x) ∼ (x − a) 2 at the boundaries, the continuum is gapped away from zero by some δ > 0 [61] and we can define a (pseudo-) inverse which determines v in terms of the inhomogeneity ζ up to the physical gauge freedom. Again, the inverse is a bounded operator on the space perpendicular to the constant function with bound D = δ −2 < ∞. This is no longer the case, however, if n(x) ∼ (x − a) p where p > 2 at the boundaries [61]. Nevertheless, since we only consider potentials that are less singular than v(x) ∼ (x−a) −1/2 , we expect that the generated density goes at most as n(x) ∼ (x−a) 2 and thus such densities allow for a unique inversion (up to a gauge). In the case that we consider the Sturm-Liouville problem on the whole real axis we are again in the case of limit-point endpoints, and there is again only one self-adjoint realization of the Sturm-Liouville operator. Unfortunately it is not known under which conditions on the density the self-adjoint operator has a spectral gap around zero and thus allows for a (pseudo-) inverse of −∂ x (n∂ x ). However, from considerations similar to [23] one can presume that for most non-vanishing densities a unique inversion should be possible.
So far we have seen that the one-dimensional case is already quite involved. Turning back to (130) in more dimensions we pose again the question of invertibility to solve for v in a certain class of potentials, respecting the appropriate boundary conditions. By turning to a weak formulation of the problem, these immediately arise. To this end we adjoin a scalar field u, thought of being from the same class as the potential v, by means of the standard L 2 (Ω) inner product.
− u|∇ · (n∇v) = u|ζ (141) The possible time dependence of all quantities is now suppressed, the equation is to hold at every instant. Now if the class of potentials is assumed to have zero or periodic boundary conditions partial integration defines a symmetric bilinear form Q by We employed the specific boundary conditions to have a vanishing boundary term after partial integration. Note that in the case of periodic potentials also the density n and ∇v have to obey this periodicity. Another option would have been to demand n = 0 at the border like in the original Runge-Gross proof [23].
The theorem of Lax-Milgram [81] now gives a direct and positive answer to the question of existence and uniqueness of a solution v of (142). Moreover the solution depends continuously on the given data ζ thus the inverse operator is bounded. For this theorem to hold, the bilinear form has to fulfil for all u ∈ H Q and thus Q is automatically continuous on the Sobolev space H 1 (Ω) with the additional boundary conditions and norm u 1,2 = u 2 + ∇u 2 . To show coercivity we need to rely on Poincaré's inequality [82, 6.30] u 2 ≤ c(Ω) ∇u 2 that is true on bounded domains Ω (or domains that can be fully enclosed between two parallel hyperplanes) and again zero or periodic boundary conditions. Thus a unique solution of the Sturm-Liouville problem can be guaranteed to lie in the given Sobolev space.
The strategy for more general densities is similar but more involved, a detailed account is given in [83]. We construct a density-adapted weighted Sobolev space with norm u 1,2,n = u 2 + √ n∇u 2 and one shows that for this space the bilinear form Q is naturally continuous and coercive if the density fulfils n −s ∈ L 1 (Ω) with s > d 2 , d being the dimensionality of Ω.p Note that the restrictions to a bounded domain Ω (real boundedness because one does not only rely on Poincaré's inequality) and zero or periodic boundary conditions are still active. A further practical consequence is that the weighted Sobolev space is compactly embedded in L 2 (Ω). Nevertheless the condition n −s ∈ L 1 (Ω) seems harsh, especially if the problem is not periodic, and will not be fulfilled by natural densities with zero boundary like a particle in a box.
While the above considerations give us conditions on n and ζ such that we can uniquely solve the Sturm-Liouville equation, we did not discuss whether the resulting potential is regular enough to generate a well-defined q [v]. However, in order to set up an iterative scheme as proposed in Sec. 3.2, we need to guarantee that q [v] exists. Details about conditions on the initial state, interaction and potentials will be given in Sec. 3.5. For the moment assume that for an appropriate v 0 we can determine via propagation a well-defined q[v 0 ]. If we are given a density n for which the Sturm-Liouville operator allows for an inversion we can construct v 1 uniquely by In a next step we then can construct v 2 accordingly and we finding from the boundedness of [−∇ · (n(t)∇)] −1 by a time-dependent constant D t that We point out that this inequality does not necessarily need to use the same norms for v and q [v]. Thus if we denote the norm on v by · and the norm on q[v] by · , we have alternatively where the constant D t depends on the norms used. This inequality will play an important role later in the fixed-point procedure.
p Note that the restriction s > d 2 was missing in the main theorem of [83]. Still everything is correct in the typical case d = 3 where the smallest applicable integer value is s = 2 as stated there.

Proof based on Taylor expansion
In the previous section we have discussed the Sturm-Liouville problem. With this we can approach the issue of proving the existence of a density-potential mapping, which is basically equivalent to showing that the mapping v → n is injective for an initial state Ψ 0 and some well-defined set V. While to show the uniqueness of a fixed point of the approach presented in Sec. 3.2 is quite involved, there is a very clever trick to establish the injectivity of v → n without too much ado. The trick was first presented in the seminal work [23] and forms the basis of the Runge-Gross theorem and its many variants for different physical situations [84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104].
We first assume an initial state that is spatially twice differentiable, i.e., Ψ 0 ∈ D(T ), and has an initial density that is zero at most at the boundaries. This is fulfilled, for instance, if we start with the ground-state of a quantum system. Then, for any two different potentials v, v ∈ V (the set defined for the mapping (100)) the fundamental equation (97) gives rise to a well-defined Sturm-Liouville equation. Since at t = 0 both systems have the same density, one can subtract both equations and finds Now, if both potentials differ by more than a gauge constant, the left hand side of the equation is non-zero. If Ψ 0 is at least four-times differentiable, i.e., Ψ 0 ∈ D(T 2 ), and the interaction square-integrableq such that Ψ 0 |q(r)Ψ 0 is well-defined and due to q( , r, 0) = 0. (149) As a consequence, the densities will be different infinitesimally later in time. Therefore, any two potentials that differ by more than a gauge at t = 0 will generate different densities and thus injectivity of the mapping v → n is established. Note that if one of the two potentials would not obey the conditions imposed by the invertibility of −∇ · (n∇), then we could not make this conclusion. In this case the application of the Sturm-Liouville operator to the difference v −v would not be defined and the equation would not exist r. Such a problem appears, for instance, in the example of [84]. There a difference function of the form v(r, 0) − v (r, 0) ∼ exp(|r|) / ∈ L 2 (R 3 ) is used, for which the Sturm-Liouville operator is clearly not defined. The initial resolution of this alleged counter example [105], i.e., to only allow for potentials that are generated by a finite charge distribution (and are thus square-integrable), is in clear accordance with the necessary conditions to ensure invertibility of the Sturm-Liouville equation.
The trick of Runge and Gross is now based on the consecutive application of the above result to the Taylor expansion in time of the different potentials. Obviously this is a first restriction, since not all possible potentials for the solution of the TDSE need to be infinitely-often differentiable with respect to time (with appropriate boundary conditions). If we assume that the wave functions are infinitely-often differentiable in time as well, then we formally find [25,26,22,21] ∇ · n (0) (r)∇v (k) (r) = n (k+2) (r) − q (k) (r) q We point out that by partial integration all the derivatives of the interaction can be shifted to derivatives on the wave function in the definition of q[v]. r Remember that if the potential is outside of the domain of the operator, then by construction where we write the k-th derivative in time at t = 0 of the different functions by as Since the different q (k) and n (k) can be calculated from only knowing v (l) with l < k (due to the Heisenberg equations for these operators), we find that if two potentials differ in order k (while they are the same for l < k) we have Thus the two potentials will necessarily lead to different densities. We point out that in order for this conclusion to be made we need to ensure that all the functions in (150) exist. If this would not be the case, we could not subtract (150) for two different potentials and rearrange them in the form of (151). Sufficient conditions are that for the initial state we have Ψ 0 ∈ D(T n ) for all n ∈ N, such that the kinetic-energy operator can be repeatedly applied from which infinite differentiability follows, and that the interaction w is in the Kato class of potentials, which includes the Coulomb interaction. This constraint on the initial state excludes initial states with cusps in the density. Such initial states occur, for instance, if we solve for the ground state of the static Schrödinger equation with external Coulomb potentials generated by point charge nuclei. However, if we soften the external potential by using finite nuclei [59] in an infinitely differentiable way the cusps in the ground state density will vanish and by using the corresponding ground state as an initial state the Runge-Gross proof is valid without changing essential physics (in fact finite nuclei are more realistic than point nuclei).
A final loophole we have to close is that there are still infinitely-often differentiable potentials which are different but all their derivatives are the same at one point, e.g., v(r, t) − v (r, t) = f (r) exp(−t −2 ) at t = 0. So we cannot conclude for these type of potentials that they will necessarily lead to different densities. To overcome this problem we restrict to only those potentials that have a converging Taylor expansion for some finite time t > 0, i.e., v(r, t) = k v (k) (r)t k /k!. Therefore, if we assume an appropriate initial state, the mapping from the set of Taylor-expandable potentials (with the appropriate boundary conditions) to densities is invertibles. This is the statement of the famous Runge-Gross theorem and provides the foundation of TDDFT.
The Runge-Gross result enables us to perform a density-functionalization of timedependent quantum mechanics. Instead of solving the full TDSE for a given initial state Ψ 0 and an external potential v we can self-consistently solve the equivalent non-linear evolution equation with the initial conditions n(r, 0) = Ψ 0 |n(r)Ψ 0 and ∂ t n(r, t)| t=0 = − Ψ 0 |∇ ·ĵ(r)Ψ 0 . However, there are two related problems we encounter at this point. Firstly, we do not know the set of v-representable densities we are allowed to vary over in search for the s All these conditions are fulfilled, for instance, in the case of an infinitely-often differentiable initial state (with n 0 ≥ > 0) on a torus, i.e., a periodic system. In this case the mapping from all spatially smooth and temporally Taylor-expandable potentials to their respective densities is invertible.
(existing) self-consistent solution. We would need a precise specification of the set of densities associated with the given set of potentials. The second problem is that while this equation would be in principle enough to do (orbital-free) TDDFT, similar to the minimization of the energy functional in ground-state DFT, it is extremely challenging to find good approximations to the operator q[Ψ 0 , n]. Especially the kinetic part of the operator (see (98)) is notoriously hard to approximate in terms of the density and initial state only. Therefore one usually wants to use an approximation to the q operator based on an auxiliary quantum system. Both problems are related to the question of v-representability. If we know the set of v-representable densities, then we can solve (153) by varying over this set, and if we can show that two different Hamiltonians have the same set of v-representable densities, then we can connect both systems by a KS construction. This allows us to approximate the (divergence of the) internal forces of an interacting problem by the q s of a non-interacting problem. In this case, the so-called Hartree-exchange- , n] would be defined (assuming that both systems have initial states Ψ 0 and Φ 0 with the same initial density and first time-derivative of the density) by Now, can we learn something about the set of v-representable densities in the Runge-Gross approach? First of all, the densities generated under the above assumptions are infinitely often differentiable in time and space at t = 0. From (150) we even have the Taylor coefficients of the density in terms of v and Ψ 0 . However, it is not clear whether the series converges and thus that the density is analytic. It does definitely not hold in general, since one can find counter examples (see Sec. 2.1 and [106,107]). On the other hand, using (150) we can also construct the Taylor coefficients of the associated potential from a time-analytic density n and the initial state Ψ 0 . Again, we cannot guarantee that the resulting series converges, i.e., that the density is v-representable by a time-analytic potential. However, if we assume that (under certain conditions) a Taylor-expandable density gives rise to a Taylorexpandable potential (and vice versa) we have a specific characterization of a set of v-representable densities. This forms the basis of the extended Runge-Gross approach presented in [25], which also has been applied to different physical situations [91,92,93,94,95,96,97,98,99,100,101,102,103,104]. In this case we can vary over all time-analytic densities to solve (153) and we can explicitly construct the Hartreeexchange-correlation potential by where n (k) is determined by v (l) and v (l) Hxc for l < k [108]. In this way we can construct the Hartree-xc potential from its Taylor expansion in time This proof therefore gives a construction of the xc potential and forms the theoretical basis for the existence of a KS system corresponding to a given interacting system. Note that the construction can also be carried out for connecting systems with two different interactions [25]. Again it will be important that at least these interactions are in the Kato class of potentials.

Proof based on a fixed-point scheme
We have started our general discussion of the density-potential mapping in Sec. 3.2 with a quite intuitive approach. Under the assumption that the Sturm-Liouville operator −∇ · (n∇) is invertible and that if two external potentials are close then their respective internal forces are close, we employed the fundamental equation (97) for a fixed time-dependent density to generate an iterative sequence of external potentials. These potentials were supposed to reproduce the prescribed density better and better via time propagation as we proceed in the iterative procedure. In Sec. 3.3 we presented conditions such that the first part of the assumption, i.e., that the operator −∇·(n∇) is invertible, holds. These considerations were enough to show the existence of a densitypotential mapping for Taylor-expandable potentials, i.e., to prove the classical Runge-Gross theorem. While the original work [23] did not have this iterative approach in mind, the result is directly applicable to this sequence of potentials. It tells us that if the sequence converges to a potential, then the fixed-point is unique within the set of Taylor-expandable potentials. Of course there could still be a second fixed-point outside of this set. While these results provide the existence of a density-potential mapping, they do not give us precise conditions for a density to be v-representable. However, for a solution of (153) or the KS approach to TDDFT we need a characterization of these densities. Put differently, we would like to know under which conditions a fixed point of the iterative procedure introduced in Sec. 3.2 exists. In this section we will discuss the existence (as well as an extension of the uniqueness results) of a fixed point with the help of a combination of (147) and an inequality, which makes the statement that two potentials which are close generate similar internal forces, precise.
Let us start by first discussing under which conditions two external potentials generate similar internal forces. From our considerations in Sec. 2.5 we know that in order to have a well-defined divergence of the internal-force density q[v] the corresponding wave function Ψ[v] needs to be at least four-times differentiable. Consequently we impose this condition on the initial state Ψ 0 and only allow for external potentials v that stabilise this property when propagating with the associated evolution operatorÛ ([v], t, 0). We denoted this set of potentials by V. While from the discussion preceding (100) we presume that potentials which are twice differentiable in space and which have some regularity with respect to the time variable will fulfil this condition, we only know it explicitly for infinitely-often differentiable potentials in space and time [58]. For simplicity we impose the same conditions on the interaction, although since we can shift the derivatives with respect to w in the definition of q[v] to derivatives of the wave function, we expect that any square-integrable interaction (for instance the Coulomb interaction) should be possible as well. Under these conditions the different potentials v are bounded functionst on Ω and give rise to well-defined (divergence of) internal-force densities q[v]. Now, for bounded potentials the wave functions Ψ[v] are Fréchet differentiable with respect to the external potential v (see Appendix B or [109]). This implies that changing the potential slightly by ∆v will affect the evolution of the initial state even less. Further, in this case the fundamental theorem of calculus for Banach spaces holds. Provided the wave functions are regular enough we can apply the fundamental theorem also to q[v] by where ∆v 12 = v 2 −v 1 and v λ = v 1 +λ∆v 12 . The linear-response kernel can be expressed explicitly in the Heisenberg picture of quantum mechanics by a commutator of the form δq([v], r, t)/δv(r , t ) = −i Ψ 0 |[q(r, t),n(r , t )]Ψ 0 [75,61]. With an estimate for the linear-response kernel on the connecting line from v 1 to v 2 , we have with the constant C t [v 2 , v 1 ] < ∞ depending on v 2 , v 1 and time. This constant stands for the bound of the linear-response kernel along the straight line from v 1 to v 2 . Inequality (157) shows that if two potentials are close in · norm over time, then their respective internal forces are close in · norm as well.
In the following we combine (147) and (157) to provide (under certain conditions) uniqueness and existence of a fixed point of the iterative sequence introduced in Sec. 3.2. We first employ a well-known trick from the theory of differential and integral equations [110,111,112]. By using the so-called Bielecki norm for an 0 ≤ α < ∞ that changes the balance of the norm towards earlier times where we denote (147) can be rewritten in terms of the Bielecki norm as We note that all α norms are equivalent since and thus they all define the same Banach space of potentials.
Assume now that we would have two fixed points of our iterative procedure F[v] = v and F[u] = u which differ by more than just a gauge. Then by taking However, this can only be true if u − v α = 0 and therefore u = v. As a consequence the iterative procedure F[v k ] = v k+1 has at most one fixed-point in V, which is equivalent to state that the mapping v → n is invertible on this domain. Finally we consider the existence of a fixed point. While before we already knew that the density n[v] is associated with a potential, since we wanted to ensure the uniqueness of a v-representable density, now we do not a priori know that the respective density is generated from a TDSE. In the case of a KS scheme, for instance, we want to ensure that a potential that was generated by a different TDSE, i.e., with a different interaction and initial state, can also be reproduced by a non-interacting system. The most basic conditions we need to impose on the density is that in accordance to Sec. 3.2 the density obeys the initial conditions Further, the initial state and the density have to be chosen such that the Sturm-Liouville operator −∇ · (n∇) is self-adjoint and invertible (up to a gauge) for the whole time interval [0, T ] (see Sec. 3.3 for details). A further minimal condition is that n ∈ L 1 (Ω) is at least four-times differentiable in space (with the appropriate boundary conditions) and two-times differentiable in time, since it should correspond to a density that is generated from a Ψ[v] which obeys the fundamental (97). Now, to set up a well-defined iterative sequence, we need to ensure that every v k ∈ V gives rise to an appropriate q[v k ] from which we can construct F : v k → v k+1 as defined in (124) and which is again in the same set. For instance, if the initial state, v k and n are infinitely-often differentiable (in space and time) then q[v k ] is infinitely-often differentiable [58], and consequently v k+1 has the same properties. For a more general set V we do not know whether q[v k ] is regular enough to guarantee v k+1 ∈ V. However, from (110) we see that a v k twice differentiable in space and continuous in time generates a temporally continuous ∂ 2 t n[v k ] and q[v k ], which would in turn generate a v k+1 having the same properties. Also the proof strategy of [48, Th. X.70] makes it plausible that potentials which are twice differentiable in space and continuous in time generate a continuous q[v k ] which would be enough to guarantee that v k+1 obeys the same conditions.
Under the above assumptions we can apply (160) to the iterative sequence. The constants C[v k , v k−1 ] are all estimates on the lines connecting two successive iteration points v k−1 and v k . If we assume α big enough such that α ≥ 2DC[v j , v j−1 ] for all iterations j ≤ k we discover that its distance from the starting point v 0 is indeed limited by the α-norm of the very first step.
The next constant C[v k+1 , v k ] can be chosen as a maximal estimate of the linearresponse kernel over a set of potentials big enough such as to include the connecting line from v k to v k+1 . From (164) we see that both v k and v k+1 will surely be included in the convex set − v 0 α } and this is true for all further steps if we assume that there exists a constant and take α = 2DC. Note that this assumption is still needed here because such a supremum does not necessarily exist on a ball in a infinitely dimensional Banach space such as the space of potentials. Then v k is a Cauchy sequence in the Banach space equipped with the Bielecki norm and therefore the resulting fixed point v = lim k→∞ v k solves the local-force equation and by the argument of (125) generates the prescribed n by propagation of the initial state Ψ 0 .
In summary, for the possible v-representability of a density via the fixed-point procedure we need that the initial state Ψ 0 is at least four-times differentiable and that the prescribed density obeys the initial conditions of (162) and (163). Further, the associated Sturm-Liouville operator has to be invertible, i.e., n is strictly positive except possibly at the boundary, ∂ 2 t n is at least continuous in time and ∇ 4 n is integrable, since it should be associated with a Ψ[v] which is four-times differentiable as well.

The density-potential mapping in lattice systems
In the course of this review we have seen that a lot of subtleties arise due to the standard formulation of quantum mechanics in terms of unbounded operators on an infinite-dimensional Hilbert space. The same subtleties are of course also found in the density-potential mappings. In order to avoid the resulting problems of the standard formulation, we can consider an approximate treatment of quantum mechanics from the start. There are two different (but closely related) approaches to do so: We can either stay within an infinite-dimensional Hilbert space of quantum states but restrict ourselves to bounded operators only (which allows to consider the C * -algebra of bounded and self-adjoint operators [33,113]), or we can consider quantum mechanics on a finite dimensional Hilbert space (which makes all operators automatically bounded). Both approximation strategies are quite general, and their explicit realizations depend strongly on the actual physical situation we want to model [5]. Most clearly we see the relation between both approximation strategies if we consider discretized models, so-called lattice systems. For instance, we can divide R 3 into (infinitely many) small boxes and approximate the square-integrable functions by their mean-value in these boxes. By then approximating the kinetic-energy operator by a finite-difference expression in terms of these mean values we have a bounded operator on the appropriate (infinite-dimensional) sequence space of square-summable functions (see for instance [15]). Also the interaction energy can be expressed as a bounded operator in this sequence space and the external-energy operator is given in terms of a (bounded) multiplication operator (also called the on-site potential). By then restricting further to a finite part of the lattice we have a finite-dimensional approximation. The resulting TDSE (without loss of generality we again disregard the spin-degrees of freedom) for N particles on the lattice with M sites for every particle reads [114,115,116] i∂ t Ψ( z 1 , ..., z N , t) = − N n=1 yn T ( z n , y n )Ψ(..., y n , ..., t) + N n=1 v( z n , t)Ψ( z 1 , ..., z N , t) where the hopping rate obeys T ( z n , y n ) = T ( y n , z n ) and T ( z n , z n ) = 0 can be assumed (since it amounts to a constant shift in the on-site potential, i.e., a gauge transformation). The existence of a unique solution Ψ ∈ C 1 ([0, T ], H d ) on the discretized Hilbert space H d for every Ψ 0 ∈ H d can then be based on the Picard-Lindelöf theorem of ordinary-differential equations (even for the case of infinitely many sites). The Picard-Lindelöf theorem (which itself is an application of the Banach fixedpoint theorem) implies that an iterative mapping converges to a unique solution provided that f is continuous in its first argument and there is a constant L < ∞ such that for all s ∈ [0, T ] and Ψ 1 , Ψ 2 ∈ H d . Since in our case f (s, Ψ(s)) =Ĥ d (s)Ψ(s) andĤ d is the bounded discretized Hamiltonian of (167), the Lipschitz constant L is the highest eigenvalue of the Hamiltonian. Equivalently one could use similar ideas as presented in the Sec. (2.4) to show the existence of a unique evolution operator. Thus, for the set of continuous and bounded on-site potentials V d we have a mapping for every initial state Ψ 0 . Accordingly we can then define mappings for observable quantities like the (on-site) density the density matrix or the link current J( z, y, t) = 2 {T ( z, y)ρ( z, y, t)} provided we have anti-symmetrized the wave function appropriately. The question whether we can establish a density-potential mapping for this discretized problem will again be based on an equation that connects n and v explicitly. In analogy to the continuum we first determine the discretized version of the continuity equation [117,116,118] ∂ t n( z, t) = − y J( z, y, t).
By construction every Ψ[v] obeys this equation (which is not automatically true for the continuum case since there we are dealing with unbounded operators). By then calculating the equation of motion for the link-current J and then combining it with the continuity equation we find the fundamental equation of lattice TDDFT where and Here we have defined the two-particle density matrix by Equation (175) is the lattice equivalent to (97) of the continuum version. We point out, that a sufficient condition for this equation to hold is to to guarantee that Ψ ∈ C 2 ([0, T ], H d ). This can be shown to be true (at least) if v is one-times continuously differentiable in time, which restricts the set of allowed on-site potentials.
Can we now define (maybe without all the mathematical trouble we face in the continuum case) a lattice density-potential mapping based on this equation? Before we answer this question, we point out that it was realized by different authors [115,114] that there is an evident lack of v-representability for certain timedependent densities on a lattice. These issues can be demonstrated most easily for a simple two-site model. The Hamiltonian for this problem iŝ where {σ x ,σ y ,σ x } are the Pauli matrices, T kin > 0 a hopping parameter, and the Hilbert space is simply H d = C 2 . The function v(t) is actually the potential difference between site 1 and 2, thus we have effectively fixed a gauge. The conjugate observable to the potential difference is the density (difference) operatorσ z (t) and therefore we would like to investigate whether we can define a σ z → v mapping with σ z (t) = Ψ(t)|σ z Ψ(t) . We do so by considering the equivalent equations to (174) and (175).
Provided v ∈ C 1 ([0, T ], R) any solution of the above two-site TDSE obeys these two equations. Obviously, if we want to realize a density-potential mapping for a fixed initial state Ψ 0 , we need to obey the initial conditions that are implied by these two equations. However, these equations also immediately give other conditions on vrepresentable densities. Firstly, from the lattice continuity equation we see that only those densities are possible to achieve, which obey a maximality condition of the form |∂ t σ z (t)| ≤ 2T kin .
(182) The hopping parameter T kin restricts the maximal change of density at each site. Further, if we want to follow the Runge-Gross idea by restricting to the set of Taylorexpandable v, we need to ensure that σ x (0) = 0 for the initial state (which excludes for instance states of the form (1, 0) or (0, 1)). In the general case of (175) this conditions is equivalent to guarantee the invertibility of the M × M symmetric matrix K(0) with entries k( x, y, 0) (in a space perpendicular to the constant function, i.e., the gauge freedom of the on-site potential). To circumvent some of these issues, in [117,116] time-dependent link-current density-functional theory was developed. However, by restricting to initial states that obey the invertibility condition onK(0) a time-dependent density-functional theory on the lattice can be formulated [118].
We start by choosing an initial state Ψ 0 for whichK(0) is invertible (up to a constant function), e.g., the ground state of a connected lattice [118]. In the above two-site example this amounts to demand σ (185) In a next step we could now express the dependence on Ψ in terms of v and the initial state Ψ 0 (which is merely a functional variable change), which would lead to the lattice equivalent of (123). A self-consistent solution of the resulting functional equation merely in terms of v would amount to the lattice version of the fixed-point approach of Sec. 3.5. We will investigate this scheme at the end of this section.
Before we do so, we follow [118] and use (185) to set up a non-linear TDSE (at least until t * ) by expressing This non-linear TDSE is the lattice equivalent to the non-linear TDSE introduced in [76,77,63,78] and also discussed in Sec. 3.2. Now, remember the Picard-Lindelöf theorem that we used to guarantee existence and uniqueness of solutions to the original TDSE (and thus the potential-density mapping). We only needed to show that the right-hand side of the TDSE obeys (169). Provided that we are within B ⊂ H d , the inversion ofK(t) perpendicular to the constant on-site potential is possible and thus the right hand-side of (187) is bounded. This local boundedness is enough for a local version of the Picard-Lindelöf theorem, which guarantees that existence and uniqueness of a solution to (187). The only restriction is that [0, T ] is supposed to be such that the iterative sequence defined by (168) does not leave B ⊂ H d , thus T depends on the chosen n.
What happens if we hit this v-representability boundary ∂B at some time 0 < t * is most easily seen in the two-site model. The resulting non-linear TDSE is due to (184) given by Obviously, if we leave B the non-linearity becomes infinite since σ x ([Ψ], t * ) = 0. This excludes the application of the Picard-Lindelöf theorem across the boundary to guarantee a unique solution. However, in this model-system we can analyse the behaviour at the v-representability boundary in more detail. By a construction similar to the construction of the explicit density-potential mapping in the continuum of Sec. 3.1, where we expressed the complex wave function in terms of its polar representation, we can deduce an explicit form of the two-site density-potential mapping in terms of σ z (t) by [118] v ± ([n], t) = ± ∂ 2 t n 1 (t) + 2T 2 kin σ z (t) 2 4T 2 kin n 1 (t)n 2 (t) − (∂ t n 1 (t)) where n 1 (t) = (σ z (t) + 1)/2 and n 2 (t) = (1 − σ z (t))/2 are the densities at each site. The ± stands for two different choices of the initial-states phase function. Whenever 4T 2 kin n 1 (t)n 2 (t) − (∂ t n 1 (t)) 2 → 0 (which is equivalent to the maximality condition of (182)) we are at the boundary between the v + and v − realization of the densitypotential mapping [118]. It would therefore be highly desirable to somehow extend the mapping across this boundary uniquely.
At this point we can profit from the fixed-point approach introduced in the previous section. Firstly we point out, that a lattice fixed-point approach amounts to consider an iterative mapping based on which is the variable-transformed version of (185). Obviously we need to impose the same assumption about invertibility ofK(t) as before. We therefore can set up an iterative sequence which is similar to the iterative sequence defined by (168), where every new Ψ k (t) defines by construction a new v k . Thus, while the above non-linear TDSE approach investigates the existence and uniqueness of a fixed point in terms of Ψ k (t) the lattice version of the fixed point approach of Sec. 3.5 considers convergence in terms of v k . If we choose the starting points of these two iterations such that Ψ 0 (t) leads to v[n, Ψ 0 ] = v 0 then both iterations are exactly the same. Therefore, both iterations lead to the same unique solution v[Ψ 0 , n] for some (appropriately short) time interval [0, T ]. But how do we guarantee that the limiting potential really reproduces the prescribed n? In particular, the Picard-Lindelöf theorem implies convergence also for densities which are unphysical, e.g., for an n with wrong initial conditions. To guarantee that the solution Ψ[Ψ 0 , n] and therefore also v[Ψ 0 , n] actually reproduces the prescribed density n we need an equation similar to (125) in the continuum case. This equation can be found in the lattice situation provided that v( z, t) is continuously-differentiable such that we can subtract (183) from (190) and end up with where ρ( z, t) = n([Ψ 0 , v], z, t) − n( z, t). If we then assume that the given n and the propagated n[Ψ 0 , v] have the same initial conditions, i.e., ρ( z, 0) = ∂ t ρ( z, t)| t=0 = 0, than necessarily n = n[Ψ 0 , v].
As a consequence we can conclude, that if the unique solution Ψ[Ψ 0 , n] produces a continuously-differentiable v[Ψ 0 , n] by (168), it reproduces the prescribed n. This implies, that if a potential v[Ψ 0 , n] crosses a v-representability boundary ∂B, then only a continuously-differentiable extension across this boundary is a sensible choice. Thus, in the situation of the two-site problem one needs to choose the extension of v ± ([n], t) according to this condition whenever σ x → 0. Whether this can always be done is not clear. Nevertheless, for any initial state that has an invertibleK(0) (for instance the ground-states of connected lattices [118]) and an appropriately short time interval we have a well-defined lattice density-potential mapping n → v.

Numerical realization of the density-potential mapping
The main ideas and concepts of density-potential mappings were developed alongside those of DFT and TDDFT. However, while DFT and TDDFT without the basic density-potential mappings would not be possible, these mappings, on the other hand, do not rely on DFT methods. Actually, they can be put to use also in other areas of physics and chemistry. For instance, in the context of quantum control theory they augment already existing techniques [119,120,121,122] to steer the dynamics of quantum systems [123,124]. In this section we discuss the density-potential mappings from the point of view of quantum control theory and consider their numerical construction.
We usually employ quantum mechanics to predict the behaviour of a microscopic system. Such a system is modelled by the initial state Ψ 0 and the external potential v which acts on it (see for instance (4)). Then we can, in principle, determine Ψ[v] from the resulting TDSE and calculate all physical observables O(t) = Ψ(t)|ÔΨ(t) . However, we can also use quantum mechanics to control the behaviour of a microscopic system [125], i.e., we can try to determine a v that forces the wave function to show a previously specified behaviour. This can be done in two different ways: The first approach optimizes a chosen observableÔ in time by varying over all possible wave functions starting from a given Ψ 0 in the functional [122] For instance, we could start with the ground state of a quantum system and then try to maximize the occupation of the first excited state, i.e.,Ô = |Ψ 1 Ψ 1 |. This approach is called quantum optimal-control theory [122,125]. In our search for an optimal wave function we cannot allow all Ψ ∈ C 0 ([0, T ], H), since not all of them are connected via the solution of the TDSE to the initial state. To enforce that we only vary over wave functions that are solutions to the TDSE we can either work with a Langrange multiplier on the wave functions, which in this case is a wave function itselfu and obeys a time-reversed TDSE [122], or one can employ that all wave functions are labelled uniquely by their respective external potentials and vary with respect to the potentials [126]. This formulation of optimal control theory [126] employs the mapping v → Ψ introduced in Sec. 2.4 directly. Either way, the resulting control equations to determine an optimal external potential v are numerically extremely demanding, since they usually imply hundreds if not thousands of global iterative solutions of the full TDSE [122]. The second approach avoids these numerically expensive global iterations by (instead of optimizing) prescribing the physical observable O(t) at every time and then try to find a Ψ that reproduces this as a solution of the TDSE with some potential v. However, not every path O(t) can be reproduced by a TDSE with a local potential only. For instance, the control of the non-local observableÔ = |Ψ 1 Ψ 1 | would need a non-local potential that can project Ψ(t) directly onto the first initial state. On the other hand, if we would like to find an external potential v that generates a prescribed dipole moment, we will find multiple solutions. In abstract terms, the mapping v → O is usually not invertible and we cannot guarantee that a prescribed path O(t) is vrepresentablev. If we, however, restrict ourselves to controlling the density, these issues can be avoided. This has to do with the fact, that the Runge-Gross theorem discussed in Sec. 3.4 guarantees the invertibility of the mapping v → n, and thus almost all densities which are consistent with the initial state are v-representable. In principle we can then apply the basic procedure of this so-called local control theorya [119,120,121], where we discretize time and determine an appropriateĤ(0) from solving where n and ∂ 2 t n are prescribed and q[Ψ 0 ] is determined from the initial state. This Hamiltonian is then used to make an Euler time step Ψ(t 1 ) = Ψ 0 − i∆tĤ(0)Ψ 0 b which gives a state that has the prescribed density (without multiple global iterations as was the case in optimal control theory). While the Euler method works well in practice for simple control objectives and small enough time steps ∆t, it is numerically very inefficient and can fail in practice due to round-off errors, which it does in the density case. This is due to the fact that the density (at a given point) may change by orders of magnitude, so that we have to be very precise to stay correct. If we do not, an u To be precise, the Lagragian multiplier will be part of the dual space of the Banach space of wave functions. Therefore in the case of optimal control theory with a Lagrangian multiplier it might be more convenient that one considers the wave functions as part of the self-dual Hilbert space of time and space L 2 ([0, T ], H), such that both, the TDSE wave function and its Lagragian multiplier, are within the same function space. v For v-representability both, the control objective and the control field need to have the same degrees of freedom, i.e, the size of their respective sets need to be the same. a We point out that controllability of the prescribed observable is a major challenge in standard localcontrol schemes. However, local-control schemes that only enforce a monotonic increase/decrease of an observableÔ can overcome most of these problems. b Remembering all the intricacies of the TDSE discussed in Sec. 2.1 one should be a little suspicious about approximating the evolution of the wave function by repeatedly applying the Hamiltonian. We will discuss this issue a little later. extremely strong artificial potential is needed to compensate for the error in the next time step, which makes the resulting algorithm unstable. How we can stabilize this local-control algorithm by employing the iterative scheme introduced in Sec. 3.2, is discussed below.
Finally we can combine local and optimal control theory, provided that we have a way to efficiently calculate Ψ for a given n (besides the local-control approach, in certain physical situations one can also use other schemes [127,128,99,129]). Due to the Runge-Gross theorem we can label the wave functions in terms of their respective densities and thus we can vary with respect to the density in the optimal-control functional of (193), i.e., While the numerical cost to optimize the functional is still the same as in the standard optimal-control approach, the restriction of the search space becomes simpler than in the previous cases. Based on one's physical intuition one can set up a basis of possible densities and then optimize with respect to these (finitely many) degrees of freedom [124]. Further, in the rare cases that observables can be expressed (approximately) in terms of the density (see also Sec. 3.1), one can determine first an optimal density from (193) and afterwards use the density-potential mapping to calculate the respective v.
Now, to stabilize the above local-control algorithm, we start by defining again the iterative procedure We choose a density n (strictly positive on Ω) that satisfies the initial conditions of (162) and (163) on some time interval [0, T ]. Then we propagate the initial state with v k and then calculating q[v k ] as defined in Sec. 2.5. A first important detail is the fact that the iterative procedure can be performed not only on the whole time interval [0, T ], on which we prescribe the density, but also successively in every subinterval of length ∆t (where we use the converged potential of the previous subinterval to determine the new initial state via propagation of the previous initial state). This partitioning of the time interval is also needed to numerically perform the time propagation [7]. If we take the time intervals small enough we can approximate the exact evolution by a time-stepping procedure with time-constant Hamiltonians (see Sec. 2.4). In principle it is then possible to determine the eigenfunctions of the Hamiltonian and calculate the propagator exp(−iĤ∆t) for the time step ∆t. Since the Hamiltonian (and hence the eigenfunctions) change in time, this procedure has to be repeated for all the successive time steps. In practice such a procedure is impossible and hence one usually adopts the approximation for some arbitrary K ∈ N. While this approximation is well-defined if we have discretized our Hamiltonian (to represent the problem on our computer), analytically the Taylor approximation (197) is usually not well-defined as has been discussed in detail in Sec. 2.1. Therefore we have to be specifically careful that this approximation does not violate any analytical constraints. For instance, we need to make sure that the wave function obeys the boundary conditions at all times, i.e., it stays within the domain of the Hamiltonian. Hence for periodic boundary conditions on [0, L] (the multi-dimensional case is straightforward) the wave function always has to stay periodic as dictated by the eigenfunctions of the self-adjoint domain exp(i2πkx/L), and for the zero-boundary case the wave function has to stay odd across the boundaries and periodic on the double domain [−L, L] due to sin(πkx/L) [124]. Thus any external potential (as well as interaction) that is applied via the Taylor approximation (197) to the wave function needs to keep this symmetry. Consequently we restrict in the following to strictly periodic potentials in the case of a periodic quantum system and in the case of zero boundary conditions we restrict to potentials that are periodic on the double domain and even across the boundaries 0 and L [124]. While these conditions come from the propagation of the wave functions and not from the iterative procedure, they are also necessary to make the iterations welldefined. This is the case, since they make sure that we can uniquely (up to a gauge) invert the Sturm-Liouville operator −∇ · (n∇) and find a new v k+1 that again allows the above time-stepping strategy (see also Sec. 3.5). To make this more precise, we first consider the case of a periodic system. For strictly positive densities we can impose periodic boundary conditions on the Sturm-Liouville operator and invert it uniquely (see Sec. 3.3 for details). Therefore we can propagate with the new v k+1 without violating the boundary conditions and perform the next iteration step. In the case of zero boundary conditions on the wave functions, we can invert the Sturm-Liouville operator provided the density does not go faster than x 2 to zero (see Sec. 3.3 for details). If we have made sure that the wave function stays odd across x = 0 then Ψ(x) ∼ x near the boundary and hence we can invert the problem. Since q and n are even across x = 0 the iterated potential v k+1 is even about the boundaries too. The resulting potential is therefore periodic on the double domain and even about the inner and the outer boundaries as required.
In a next step we get rid of the somewhat complicated term q[v k ] by employing (97), which leads to Since we only make one time step, the target density n and the iterated density n[v k ] are usually (even in the first iteration step, where we just take the converged potential of the previous subinterval) very close, such that we can approximately write (use (198) in (196)) Here a further important detail to stabilize the numerical procedure has to be taken into account. The above update formula changes the potential according to how strongly the iterated density differs from the target density. In terms of the iterated wave function this means that we directly control the modulus of the wave function, however its phase we only control indirectly. This indirect control only holds in the exact case, where the phase of a wave function corresponds to its current j (see (91) for the definition), and is determined by the modulus via the continuity equation (90). However, since we will employ a discretization of Ω and will have numerical errors in our algorithm, the continuity equation will not be fulfilled exactly. To avoid that the density is almost exact in the iteration (and numerically we would interpret the potential as converged) while the current is still far off we make the control of the phase explicit by where µ is a non-zero parameter at our disposal. The last term on the right-hand side of (200) therefore measures how well we obey the continuity equation. Now we make the time-stepping explicit. We use instead of the simple on-point Hamiltonian in the original local-control algorithm a mid-point Hamiltonian and thus a mid-point potentialv k (r, t i+1 ) = v k (r, (t i+1 +t i )/2) to make a time-step ∆t = t i+1 −t i [123,124]. For the finite-difference approximation to the time-derivatives in (200) we only employ times prior to the current time. Since for all prior times by assumptions we have converged to the exact density and current we end up with wheren is the mid-point density and A and B are constants depending on the discretization scheme of the time-derivatives and the µ of (200), which effectively leaves the choice of their values at our disposal. Usual choices are A and B between 0.5 and 1 [123,124]. Finally we (equidistantly) discretize Ω and use a (usually seven-point) finitedifference approximation for the spatial derivatives in the Hamiltonian as well as in the above update formula. Due to the fact that we can treat the zero-boundary case in the same manner as a periodic system with double the period, we have the same accuracy in the derivatives at every point of our grid. This allows us to determine the wave functions and the respective iterated potentials to a high accuracy everywhere (also at the boundaries). Especially when the density changes by orders of magnitude (at a point) we need to be very precise, since errors are compensated by the iterative algorithm in the next time step with a large and unphysical potential, which can lead to instabilities [123,124]. By smoothing the iterated potentials, unphysical and extreme differences can be suppressed. For the actual numerical propagation of the wave function the Lanzcos method is employed, since it is the most versatile and numerically cheap approach. The inversion of the discretized Sturm-Liouville operator can be performed with relaxation methods or more efficiently with multi-grid methods [130].
The above local-control algorithm is stable and can treat rapidly changing (by orders of magnitude) densities. It is independent of the dimension of the problem, the number of particles and the initial state as well as the interaction. In practice, the main obstacle to perform this local-control scheme is to store (and then propagate) the interacting many-body wave function, since one quickly runs out of computer memory. The individual update-cycles (201) are usually converged (with respect to the change of the potential differencev k+1 −v k ) within a few iterations [123,124]. However, for non-interacting problems, where the propagation can be performed efficiently, a combination with TDDFT approximations to the exchange-correlation potential allows to find also approximate potentials for large interacting systems [124].
In the following we present a few illustrative examples of the density-potential mapping constructed via the above algorithm. We first consider the system of the two interacting particles already introduced in Sec. 3.1. While before we were interested in the rigid charge transfer, we now want to do a little more and force that the density of the interacting two-particle system changes from the ground-state density n 0 (of the potential v 0 given in (119)) over time to the density n 1 of the first excited state (displayed in Fig. 7). The first excited state is a charge-transfer state, where roughly half of the density is at the left site and the other half is on the right site. We first split the charge in a similar manner as has been done in the example of Sec. 3.1 and then slowly change the density to the one of the charge-transfer state (see Fig. 8).
x (a.u.) t (a.u.) The resulting external time-dependent potential (where we subtracted the static potential, i.e., v ext = v[Ψ 0 , n] − v 0 ) that does this in the interacting system is shown in Fig. 9. The peaks at the boundaries are not numerical artefacts. This complex structure is the same for different spatial and temporal grids and indeed is needed to enforce that the rapidly moving wave function does not change its form too fast. If we instead look at a non-interacting system with the same density profile (starting from the KS ground state with a single orbital of the form ϕ 0 [0, n 0 ] as given in Sec. 3.1.) we find that the external potential (again we have subtracted the static potential v 0 ) that enforces this charge-transfer behaviour does not have these extreme features (see Fig. 10).
These examples demonstrate the capability of this numerical realization of the n → v mapping and also show that the basic ideas of the density-potential mapping can be used in practice also beyond TDDFT and the KS construction.

Extensions to vector potentials and photons
The ideas of the density-potential mappings based on the Runge-Gross approach and its extension by van Leeuwen (see Sec. 3.4) have been applied to a lot of different physical situations beyond the ones described by the standard Hamiltonian of (2), e.g., to superconducting systems [89] and to open quantum systems [93,94,95,96] (for a list of references see Sec. 3.4). In this section, we consider the application of these ideas to quantum systems driven by an external vector potential as an important example, which gives rise to (vector-)potential-current mappings. These mappings form the basis of time-dependent current-density-functional theory (TDCDFT) [84,91]. This density-functional approach can be easily extended to also include the interaction with photons [90,97,100,104].
In all our previous considerations we have neglected two important physical facts: relativity and photons. In principle we should use a kinetic-energy operator that is consistent with special relativity and the charged particles should interact via photons. The standard approach that takes these two requirements into account (and implies spin as well as the existence of positrons) is quantum electrodynamics (QED) [131,12]. While the predictions based on QED are extremely accurate, the theory has severe mathematical problems which express themselves, for instance, in divergent perturbative expressions [132]. Despite these issues, the QED Hamiltonian (or equivalently its Lagrangian) is usally employed as a starting point to derive different approximate quantum theories which describe the properties of charged particles and photons in certain limits. For instance, if we assume that the energies of the charged particles are small compared to mc 2 , then a non-relativistic treatment of the charged particles based on the Pauli-Fierz Hamiltonian [133] is justified. If we further assume magnetic fields to be negligible and take the photons in Coulomb gauge (the polarization is restricted to the two transversal degrees of freedom and thus ∇ ·Â = 0) [132] we end upc with a Hamiltonian [5,104] − dr Ĵ (r, t) · a ext (r, t) + j ext (r, t) ·Â(r) that describes electrons subject to an external vector and scalar potential, i.e, a ext and v respectively, and photons subject to an external charge current j ext and charge density n ext . The interaction between photons and electrons is described with the termsŴ + Ĵ ·Â, wherê is the charge current, the total vector potential isÂ tot =Â + a ext , and the vectorpotential operator is given bŷ where (k, λ) are the two transversal polarization vectors [132,104] and the creation and annihilation operators obey [â k λ ,â † kλ ] = δ(k − k )δ λλ . The total scalar potential is given by v tot = v + dr n ext (r , t)/4π|r − r | and the free-photon energy operator isĤ EM = dk k 2 2 λ=1â † kλâ kλ . The initial state Ψ 0 in this case is a combined initial state of electronic and photonic degrees of freedom. Now, we first consider situations where the coupling term Ĵ ·Â between the photons and the electrons is negligible. This is the case, if the initial state is separable into a purely electronic and photonic state and the transversal part of the internal charge currentĴ can be discarded. In this situation the main contribution comes from the longitudinal internal charge current for which by partial integration the coupling term is zero due to the Coulomb-gauge condition (the interaction between the electrons by the longitudinal charge current is taken into account fully by the Coulomb termŴ ). We therefore can decouple the electronic and the photonic degrees of freedom and have approximately [5,91] c Without further restrictions it cannot be guaranteed that the Pauli-Fierz Hamiltonian is welldefined [134,135]. Nevertheless, if we restrict our considerations to a box with periodic boundary conditions and introduce a highest allowed photon frequency, then the resulting Hamiltonian is selfadjoint. The following arguments, however, do not depend on this procedure and we therefore neglect these subtleties [100,104].
where now also the internal charge currentĴ =ĵ −n a ext does no longer depend on the photon field. If we only allow for scalar external potentials v (and thus a ext = 0) we rederive our original Hamiltonian given in (2). On the other hand, if we keep the external vector potentials explicit we see that with respect to our previous considerations we have more freedom in the choice of our external fields to describe and control a quantum system. Consequently instead of a mapping from potentials to densities in this case we find a mapping of the form where we assume Ψ 0 and the allowed external fields regular enough (for instance infinitely-often differentiable). The continuity equation in this case becomes then and the charge current obeys a local-force equation of the form [5,104] (suppressing all dependencies) where Q k is defined as in Sec. 2.5 and summation over multiple indices is implied.
In the case of the density-potential mapping based on the Hamiltonian of (2) the densities n and the potentials v are functions with the same degrees of freedom. Now, however, we have much more freedom since we can choose (v, a ext ). A reasonable choice to set up a similar mapping would be (n, J) since this pair would have the same degrees of freedom. If we consider the continuity equation (208) together with a prescribed initial state which implies an initial density n 0 , then we see that the charge current determines the density uniquely. Therefore, we need to have a similar reduction of freedom in the external fields if we want to have an invertible mapping for the charge current J. To find this restriction we first take a look at the case of (2), where we can add a time-dependent yet spatially constant function to v and still have the same density. For the inversion of v → n we have to restrict this freedom by fixing a gauge. Now we find that for any differentiable Λ both the pair (v, a ext ) and lead to the same charge current J and physical observables [91,5,104]. By fixing this gauge freedom we find the desired restriction and thus only take into account physically inequivalent external fields. In our case we choose the radiation gauge which fixes v = 0 and which leaves Ψ 0 unchanged by taking the initial condition Λ(r, 0) = 0 [91]. This slightly simplifies the local-force equation (209). Following now the same steps as for the fundamental equation (97)  and to provide the construction of the time-analytic vector potential a ext for a given time-analytic charge current J [91]. We point out that one could start the very same construction directly with the definition of J = j − na ext since there the external vector potential appears already explicitly. That means that we suppress the timederivatives of the complicated terms in (209) and collect them in time-derivatives of j and n [116,104] where again the superindex (k) referes to the k-th time-derivative at t = 0. And consequently we have for vector potentials a ext = a ext (by more than a gauge) which first differ in the k-th order of their respective Taylor expansions that provided n 0 is non-zero everywhere (except maybe at the boundaries) and thus J[a ext ] = J[a ext ]. This forms the basis of TDCDFT and allows us perform a selfconsistent calculation in terms of the charge current J instead of considering the full TDSE with the Hamiltonian of (207) [22,21]. If we discretize the Hamiltonian (207) [117] a rigorous iterative approach to the current-potential mapping similar to the one presented in Sec. 4 has been established by Tokatly [116].
If we do not assume that the photons are negligible, we have further degrees of freedom in the external variables since we can also choose different external charge densities and currents, i.e., n ext and j ext respectively. From the purely photonic limit of (202), we see that j ext couples to A, and thus an invertible mapping j ext → A seems possible. However, we have restricted the freedom of A by the Coulomb-gauge condition to only transversal degrees of freedom and thus we need to find a similar restriction also for the external current j ext . If we determine the equation of motion for the vector-potential operator governed by the full Hamiltonian (202) we find [104] A + ∇ dr ∂ t n ext (r , t) + ∂ t n(r , t) 4π|r − r | = j ext + J.
which only guarantees the Coulomb gauge for A if we impose a continuity equation on the external charge current and density, i.e., ∂ t n ext = −∇ · j ext .
Then the equation of above can be rewritten as A − ∇ dr ∇ j ext (r , t) + ∇ J(r , t) 4π|r − r | = j ext + J, where the second term on the left-hand side cancels explicitly any longitudinal component of A which would arise due to longitudinal components of j ext + J [104]. Consequently, any two inhomogeneities that differ only by a longitudinal function lead to the same internal vector potential A, and hence A and j ext have the same degrees of freedom. Therefore, similar to the radiation-gauge condition on the external vector potentials, we only take into account physically inequivalent external charge currents and restrict to those j ext that differ by more than a longitudinal current. As a consequence of this condition, for a given internal pair (J, A) there exists a unique external charge current j ext determined via Maxwell's Eq.(217). And since we can by a Taylor expansion similar to (213) uniquely determine all the Taylorcoefficients of the (radiation-gauged) external vector potential from (J, A) we have an invertible mapping (a ext , j ext ) from the set of Taylor-expandable external currents and potentials to the corresponding internal currents and potentials. This allows us to perform a densityfunctional treatment of a system of charged particles coupled to photons. In a similar manner also other Hamiltonians describing particle-photon systems give rise to an invertible mapping from external fields to internal fields [90,97,100,104]. If we discretize the Hamiltonian of the charged particles and keep finitely many photonic modes then a rigorous iterative formulation similar to Sec. 4 can be provided [136].

Summary, open questions and outlook
In this topical review we discussed in detail various aspects of the existence, uniqueness, and construction of the density-potential mapping in TDDFT. This problem splits into a number of important subproblems. The most basic of these subproblems is to determine the class of potentials and initial states for which the TDSE has a unique solution. We identified an important class of potentials for which a solution can be guaranteed for any normalizable initial state Ψ 0 . These are the potentials in the Kato class for which also the time-derivative is in the Kato class. We denoted this set of potentials by V. If the initial state is in the domain D(T ) of the kinetic energy operator such thatT Ψ 0 is normalizable then we can also guarantee that the solution stays in this domain and thus has finite total energy. Therefore it is natural to consider this class V as the class of physical potentials. For any potential in V and any initial state in the domain ofT we can find a time-dependent wave-function Ψ(t) and subsequently a density n(t). This defines a mapping from potentials to densities. Our next problem was to know whether this map is invertible or whether it is possible that two potentials map to the same density. An important role in this discussion was played by the local-force equation which gives a direct relation between densities and potentials. We found that for this equation to be well-defined we need to put extra conditions on the potentials. The resulting potentials are similar to those that are generated by charge distributions as calculated from the Poisson equation. Although point charges are not allowed this includes the important physical case of finite atomic nuclei (softened Coulomb potentials). Then we discussed how the local-force equation played an important role in several of the proofs of the density-potential mapping. The original proof by Runge and Gross had to require that the potential was a real-analytic function in time and that the initial state Ψ 0 was infinitely differentiable with respect to spatial coordinates. To remove the condition of temporal analyticity we considered an iterative solution of the local-force equation. This lead us to consider a number of other issues, such as the invertibility of a certain Sturm-Liouville equation and the linear response of the q-operator. We showed how we could use the iterative scheme to prove the existence and uniqueness of the density-potential mapping under certain conditions on the densities and initial states. We further showed how we could define the fixed-point procedure on a lattice and presented a numerical implementation of the fixed-point scheme to construct the density-potential mapping and gave several examples. We finally discussed a TDDFT extension to vector potentials and photons. As is clear from this summary the question whether a time-dependent density can be obtained from some TDSE has many aspects and therefore it is natural that there are still several open issues. The main issue for a rigorous approach to the density-potential mapping in terms of a fixed-point procedure, or equivalently in terms of a nonlinear TDSE, are the properties of the q-operator. We do not yet know sufficient conditions for the differentiability and boundedness of its response functions. However, recent results for the wave function [109] indicate that such conditions should be possible. Closely connected with these issues is also the question of determining the most general set of external potentials for which the fundamental equation of TDDFT (97) is welldefined. This implies the problem of ensuring that an initial state which is fourtimes differentiable keeps this property through time. For now we only know this to be rigorously true in the case of periodic systems with infinitely differentiable potentials and interactions [58]. The next open problem is to guarantee that the iterative procedure to determine the potential for a fixed initial state and density really does reproduce this density as a solution of the TDSE. This holds true if the solution to (125) for zero initial conditions and a general potential indeed is the zero function. Up to now we only know this to hold for analytic potentials. And finally there is the question whether one can extend the invertibility of the Sturm-Liouville equation to all of R 3 . This would be desirable since then one can treat the standard setting of quantum mechanics.
While we face a lot of mathematical challenges if we want to treat the densitypotential mapping in the most general setting, in practice we are usually safe. First of all, the restrictions we had to impose to make the density-potential mapping rigorous, e.g., to only allow for infinitely differentiable initial states in the Runge-Gross theorem, do not really matter when actually solving the TDSE. The time-propagation of a smooth approximation to an initial state with a cusp (as for ground states of Coulomb systems) can be made arbitrarily close to the exact propagation (below any numerical accuracy). Further, since we need to put the TDSE on a grid (or finitely many basis functions) to perform a calculation we are considering indeed an approximation to the original problem in terms of a lattice. In this case we can rely on the results of [118] to guarantee a well-defined density-potential mapping. Therefore, in practice the only real obstacle is to find better and more reliable approximation to the xc potential. Besides going beyond the usual time-local approximations and also include previous times [22], new and promising routes are currently being developed by also employing different initial states and new auxiliary systems such as strictly-correlated electrons [137,138,139].
However, answering the important open questions in the density-potential mapping will lead to new insights into the fundamentals of TDDFT and the KS construction and hopefully will also inspire more accurate approximations to the xc potential and other time-dependent density functionals. of probabilities, but become relevant in this work also as the domains of the densitypotential mapping. They consist of all Lebesgue-measurable functions f : Ω → R or C with finite L p -norm, i.e., The special case p = ∞ represents all functions that are bounded up to a set of measure zero and the associated norm f ∞ is given by the smallest such bound. If we take f roughly with amplitude A and non-zero on a volume V then the L p -norm measures the quantity AV 1/p . This means that lower L p -spaces allow more singularity while higher ones are more forgiving towards spreading, also expressed by the (continuous) embedding L p (Ω) ⊂ L q (Ω) if p > q on bounded domains Ω where the spread is not an issue. As all those spaces are normed vector spaces with always converging Cauchy sequences (completeness) they form proper Banach spaces. In the case p = 2 the norm is directly linked to the usual inner product by f |f = f 2 (note that we typically omitted the index 2 in the norm in this important case) and L 2 (Ω) has all the structure of a Hilbert space that has risen to eminent prominence within quantum theory.
The related class of Sobolev spaces includes the weak derivatives of several orders into its definition. Thus not only amplitude A and volume V are measured but also frequency N with the sensitivity controlled by a parameter m ≥ 0 defining up to what order derivatives get included into the Sobolev W m,p -norm that will be concerned with the quantity AV 1/p N m .d This is already a strong indication towards the important Sobolev embedding theorems relating such spaces. For the definition of the W m,pnorm we use a multi-index notation for the weak α-th partial derivative and note that other equivalent definitions are possible.
f m,p = |α|≤m D α f p Again the case p = 2 yields Hilbert spaces denoted as H m (Ω) = W m,2 (Ω). For bounded intervals I ⊆ R of the real line the relation between absolute continuity and the Lebesgue integration leads to the identification W 1,1 (I) = AC(I). The possibility of unique continuous continuation to the boundary points means we can give meaningful boundary conditions. For m > 1 it holds that H m (I) ⊂ H 1 (I) ⊂ W 1,1 (I) and we used the notation H m 0 (I) for functions in H m (I) with zero-boundary conditions up to the (m − 1)-th derivative, deviating here from standard notation for closed intervals I. In a multi-dimensional setting one defines W m,p 0 (Ω) as the closure of the test functions under the W m,p -norm. The definitive resource on almost all topics relating to Sobolev spaces is [82]. and try to get rid of the unbounded kinetic partT by unitary transformation. Here f (t) = f (r 1 , . . . , r N , t) is a scalar potential acting as a multiplication operator in spatial representation, including all interactions (possibly also of more than two particles) as well as external potentials. We define the unitary free evolution operator U 0 (t) = exp(−iT t) that solves the corresponding Cauchy problem i∂ t Ψ(t) =T Ψ(t) for any initial state Ψ 0 ∈ H (even if the initial state has infinite energy). With the help ofÛ 0 (t) we then perform a unitary transformation Ψ(t) =Û 0 (t)Ψ(t) of the original problem (75) to i∂ tΨ (t) =f (t)Ψ(t) withf (t) =Û † 0 (t)f (t)Û 0 (t) (this is one possible form of the so-called interaction picture). Integrating this problem over time and transforming it back we find the mild form of the time-dependent Schrödinger equation (mild TDSE) The idea is now to guarantee unique solvability of (B.2) by recursively putting Ψ into the integral (thus generating all possible paths of interaction) and showing that this mapping is a contraction and therefore has a unique fixed point. In practice it is enough to show boundedness of the mapping with an estimate involving T then taking the time interval [0, T ] short enough and finally extending to arbitrary time intervals with a continuation procedure like in [54]. To do so, this demands for a purpose-built Banach space of trajectories, i.e., a wave function for all times in [0, T ], that is a subspace of C 0 ([0, T ], H). An important stepping-stone towards such spaces is the Strichartz estimate for solutions to the free Schrödinger equation using the spacetime norm of L θ ([0, T ], L q (R 3N , C)) with (θ, q) fulfilling a certain relation called Schrödinger-admissible [54].
U 0 Ψ 0 q,θ ≤ const · Ψ 0 2 This estimate exhibits a certain smoothing property of the free evolution. An equivalent result for non-free evolution is readily achieved by the fixed-point procedure described above. With the trajectory confined by the given inequality we have C 0 ([0, T ], H) ∩ L θ ([0, T ], L q (R 3N , C)) as the Banach space of quantum trajectories. Note however that the spatial domain here is R 3N , Strichartz estimates for bounded domains are available, although not in this general form.
The set of allowed potentials for the mild TDSE to hold is then a complementary Banach space chosen in a way that f · Ψ is in the topological dual of the trajectory space. A physical consequence to that is Ψ(t)|f (t)Ψ(t) < ∞, i.e., finite energy from the potential.
One drawback of this approach is that it does not include singular Coulombic potentials (as an interaction term or external potential) if more than two electrons in R 3 configuration space are involved [109]. Yet it is general enough on the other hand to include sudden switch-on processes.
e We point out that an even more general definition of a solution to the TDSE would be possible if we defined the time-derivative in a weak sense. However, since such weak solutions are not defined at every instance in time (and thus violate the usual notion of a physical wave function), they are commonly disregarded in physics literature. [48,140] The mild TDSE (B.2) is also the starting point for the study of functional variations of trajectories. These are formed by varying the potential f within its Banach space mentioned above. To fix notation Ψ[f ] is the solution of (B.2) for a fixed initial state and potential (internal and external) f . We form the directional derivative at f in direction g by This variational derivative can be shown to be continuous in f as a linear and bounded mapping from its second argument g to variations of trajectories in the trajectory space and is thus a proper Fréchet derivative [109]. Application of this formalism can be carried over to observables and quantities such as the one-particle density and leads to the well-known non-equilibrium version of Kubo's formula. It further gives important justifications for non-equilibrium density-response theory (for an introduction see [5]) and facilitates the apparatus of variational calculus in the TDDFT context.