Local Control and Representability of Correlated Quantum Dynamics
Abstract
We present a local control scheme to construct the external potential that, for a given initial state, produces a prescribed timedependent density in an interacting quantum manybody system. The numerical method is efficient and stable even for large and rapid density variations irrespective of the initial state and the interactions. The method can at the same time be used to answer fundamental representability questions in densityfunctional theory. In particular, in the absence of interactions, it allows us to construct the exact timedependent KohnSham potential for arbitrary initial states. We illustrate the method in a correlated onedimensional twoelectron system with different interactions, initial states and densities. For a KohnSham system with a correlated initial state we demonstrate the interplay between memory and initialstate dependence as well as the failure of any adiabatic approximation.
pacs:
31.15.ee, 32.80.Qk, 71.15.MbThe representability question is one of the outstanding problems in density functional theory (DFT) DFT (); Chayes1985 (); Lammert2010 (). In timedependent DFT (TDDFT) RvL1999 (); Baer2008 (); Li2008 (); Verdozzi2008 (); KurthS2011 (); CarstenBook (); GFPP (); Farzanehpour2012 (), the question is whether, for a given initial state, there exists a local external potential that yields a prescribed density by solution of the timedependent Schrödinger equation (TDSE). In the case of a noninteracting system, representability amounts to the existence of a KohnSham (KS) system DFT (); CarstenBook (). The KS system has played a major role in the study of correlated manybody systems as it allows for the treatment of interacting systems in an effective oneparticle framework. This feature greatly reduces computational costs Andrade2012 () and (TD)DFT has hence been one of the leading methods in electronic structure theory Burke2012 (). In practice the accuracy of the method is limited by the approximate nature of the density functionals that are used. In TDDFT the most commonly used density functionals are based on the adiabatic approximation in which the KS potential only depends on the instantaneous density. These functionals can, however, fail in important cases CarstenBook () and therefore there is a great need for better functionals. To develop and benchmark such new functionals the availability of exact timedependent KS potentials is highly desirable. Although such potentials can be constructed in special cases Lein2005 (); Godby2012 () no general practical scheme has been available so far. In this Letter we provide such a scheme based on a recently introduced fixedpoint formulation of TDDFT GFPP (); GFPP 1D (). It is at the same time an efficient local control scheme based on the density which augments other important control methods OCT1 (); OCT2 (); LCT1 (); LCT2 () already of extensive use in laser physics, quantum optics Mancini2005 () and the physics of ultracold gases CHU02 (). The scheme is closely related to existing methods LCT1 (); LCT2 () but targets a spatially extended quantity instead. It is applicable to general interactions and initial states and can deal with fast and large density changes. We demonstrate the approach for an interacting onedimensional twoelectron system with different interactions, initial states and densities. For a KS system with a nonseparable initial state we illustrate the connection between memory and initial state dependence as well as the failure of any adiabatic approximation.
The global fixed point method. We consider a electron system with a timedependent Hamiltonian , where is the kinetic energy, the timedependent external potential and the manybody interaction (which may even be timedependent). The expectation values and of the density and current operators (atomic units are used throughout)
satisfy equations of motion given by
(1)  
(2) 
Here the internal local force is defined by
where is the timedependent manybody state obtained from the TDSE with potential and given initial state. Eqs.(1) and (2) imply
where is regarded as a functional of through the state . For a fixed density and initial state this is an implicit equation for the potential. To solve this implicit equation we define an iterative sequence of potentials by the iterative solution of
(3) 
In previous works GFPP (); GFPP 1D () we proved, for general initial states and interactions,
that under mild restrictions on the density the sequence converges in Banach norm sense to a potential
which is both fixedpoint of the equation and produces the prescribed density .
Although the fixed point method itself is welldefined it is highly nontrivial
to develop a stable numerical algorithm. To do this we found it advantageous
to make explicit use of also the
current (still being a functional of the density).
This is most easily done for onedimensional systems since the continuity Eq.(1) can be integrated analytically.
We find where
(4) 
and where is an arbitrary point. For this reason, and for simplicity of presentation, we restrict ourselves to the onedimensional case in this Letter. We first show how we can eliminate the quantity from our equations. By integrating Eq.(3) and using Eq.(4) we obtain
where is an integration constant. From Eq.(2) for a system with potential we then find
(5)  
To obtain an equation that is only dependent on densities we can use Eq.(1) and Eq.(4) to find
(6)  
where is a new constant. While mathematically equivalent Eqs.(5) and Eq.(6) are not numerically equivalent as their discretizations on a spacetime grid generally differ. In practice, it is therefore advantageous to use
(7) 
as follows immediately by multiplying Eq.(5) by and Eq.(6) by and
adding the results. Here is a parameter at our disposal and is a new constant.
This equation defines an iterative procedure to determine from .
The constant in this equation is uniquely determined by the spatial boundary conditions on
(and hence depends on ).
When then and and Eq.(7)
implies
that . Since we also obtain after convergence
we can calculate any observable, and in particular the current .
This is an explicit realization of the RungeGross result RungeGross () that any observable
is a functional of the density and the initial state.
Numerical procedure.
The iterative method based on Eq. (7) should be implemented stepwise in time
for high efficiency.
We use a midpoint based timestepping method
which uses the midpoint potentials
to propagate the wave function on a timegrid with timepoints .
For this we implemented the Split Operator and Lanczos method Leforestier1991 ().
Let us now suppose that we have obtained ,
and hence the giving the required density for .
Then to determine , and hence ,
we define an iterative procedure
in which we guess an initial potential and loop over potentials until we converge to the desired :

Use to calculate from by timestepping.

From calculate and .

Calculate from
(8) where .
Eq.(8) is obtained from Eq.(7) by a discretization w.r.t. time using only times with for the derivatives. We further used the fact that we have already converged up to time and replaced by on the left hand side of the equation since we found that this does not affect the convergence. The constants and depend on the discretization scheme and the of Eq.(7), which effectively leaves the choice of their values at our disposal. The constant in Eq.(8) depends on the boundary conditions and hence the geometry of the system. Below we will present examples for a periodic system. In that case the constant is determined by the periodicity condition on the spatial interval . This yields the condition
(9) 
During timepropagation numerical errors can build up in the modulus and phase of the wave function. The density is determined by the modulus only while the current also depends on the phase. This implies, for example, that errors in the phase can lead to inaccurate currents while still producing an almost correct density. This tends to happen in the case that since in that case the procedure enforces the correct density without constraints on the current as can be seen from Eq.(8). The opposite happens in the case that . By taking nonzero values for and we control the accuracy of both the density and the current and hence the modulus and phase of the wave function. This suffices to stabilize the algorithm in most cases. Perfect stability is obtained by spatially smoothening the potentials as they are obtained. Without smoothening the best algorithm is obtained when dominates for nonzero while these values are not so important with smoothening (we used and ). We find that iterations generally suffice to converge and that the precision of the potential is limited mainly by the timestepping method for the wavefunction (assuming a sufficient spatial resolution). By increasing the precision thereof almost arbitrary precision can be achieved even when the density changes by orders of magnitude.
Translating and splitting a given density. To illustrate the algorithm we consider two electrons on a quantum ring of length over a time period of length . We start by calculating the singlet ground state (which has a spatially symmetric wave function) of a (properly periodic) Hamiltonian with external potential and interaction given by
where is the interaction strength. The ground state density is denoted by . We then construct the (spatially periodic) timedependent densities and by:
The density describes a situation where the initial density is rigidly translated around the ring exactly once
whereas the density describes a situation where the initial density is split in equal halves
that are rigidly translated in opposite directions to rejoin at times and .
We have used our algorithm to calculate the potentials that
produce these prescribed densities and via timepropagation of the initial state by the TDSE.
This was done for the interaction strengths and .
In Fig. 1 we present the corresponding potentials and densities (insets).
We see large differences in the potentials for the interacting case (panels (c) and (d))
as compared to the noninteracting case (panels (a) and (b)).
The convergence of our algorithm shows that the prescribed densities are indeed representable
and that the algorithm can be used for density changes of orders of magnitude.
Exact KS potential for a nonseparable initial state. As a second example we construct an exact KS system, i.e. a noninteracting system having the same timedependent density as that of an interacting reference system. For the KS system we also need to specify an initial state with the correct initial density . This state does not need to be the KS ground state (the ground state of a noninteracting system with density ) as the RungeGross theorem RungeGross () allows for general initial states (see for further discussion ElliotMaitra ()). Here we take the KS initial state to be identical to the true correlated ground state of the interacting system. As the interacting reference system we consider a system forever kept in the ground state of the previous example for . The density is therefore stationary and equal to . Since is not an eigenstate of a noninteracting system the KS state and potential will in general be timedependent, but in such a way that they still produce the static density . We denote the KS potential by and the KS Hamiltonian is thus given by
We have determined the timedependent potential with our algorithm and displayed it in Fig.2. The square of the corresponding KS wave function is displayed in Fig.3 at four times corresponding to extreme values of the KS potential. We see strong internal motions in the wave function as it passes through states in which the electrons are wellseparated and states where they are confined to the same region in space , although the corresponding density is completely static. The wave function at these times as well as the intermediate times and are in correspondence with the extreme values of the potential in Fig.2. We also see the exact potential cannot be an adiabatic functional of the density, and hence must have memory, as an adiabatic functional produces a static potential when we insert the exact density, in conflict with Fig.2. This can be illustrated further by choosing as (the spatial part of) the initial KS state a separable state of the form
(10) 
with . In this case the KSpotential is static and given by
(11) 
up to an arbitrary constant.
In this case the exact KS potential is static as would also have been predicted by any
adiabatic approximation.
This explicitly demonstrates the interplay between memory and initial statesMaitra2002 (); MaitraBook ().
Outlook.
We presented a stable and fast algorithm to construct the external potential that, for a given initial state,
produces a prescribed timedependent density in an interacting manybody system.
The method will be valuable for further development of density functionals and local control theory.
Especially exciting is the possibility to use more advanced (multiconfigurational) initial states in DFT
in combination with existing and new approximate functionals and to test them using our benchmarking algorithm.
This can open up new possibilities for the study of strongly correlated systems within a DFT framework.
Acknowledgement. S.E.B.N. acknowledges support from the Lundbeck Foundation. M.R. acknowledges support by the Erwin Schrödinger Fellowship J 3016N16 of the FWF (Austrian Science Fonds). We further thank Prof. J. Olsen for valuable discussions.
References
 (1) R.M. Dreizler and E.K.U. Gross, Density Functional Theory  An Approach to the Quantum ManyBody Problem (SpringerVerlag, 1990).
 (2) J.T. Chayes, L. Chayes, and M.B. Ruskai, J. Stat. Phys. 38, 497 (1985).
 (3) P.E. Lammert, Phys. Rev. A 82, 012109 (2010).
 (4) R. van Leeuwen, Phys. Rev. Lett. 82, 3863 (1999).
 (5) R. Baer, J. Chem. Phys. 128, 044103 (2008)
 (6) Y. Li and C.A. Ullrich, J. Chem. Phys. 129, 044105 (2008).
 (7) C. Verdozzi, Phys. Rev. Lett. 101, 166401 (2008)
 (8) S. Kurth and G. Stefanucci, Chem. Phys. 391, 164 (2011)
 (9) C. A. Ullrich, TimeDependent DensityFunctional Theory (Oxford University Press, 2012).
 (10) M. Ruggenthaler and R. van Leeuwen, Euro Phys. Lett. 95, 13001 (2011).
 (11) M. Farzanehpour and I.V. Tokatly, arXiv:1206.6267v1 (2012).
 (12) X. Andrade et al., J. Phys.: Condens. Matter 24, 233202 (2012).
 (13) K. Burke, J. Chem. Phys. 136, 150901 (2012).
 (14) M. Lein and S. Kümmel, Phys. Rev. Lett. 94, 143003 (2005).
 (15) J. D. Ramsden and R.W. Godby, Phys. Rev. Lett. 109, 036402 (2012).
 (16) M. Ruggenthaler, K.J.H. Giesbertz, M. Penz, and R. van Leeuwen, Phys. Rev. A 85, 052504 (2012).
 (17) I. Serban, J. Werschnik, and E. K. U. Gross, Phys. Rev. A 71, 053810 (2005).
 (18) W. Zhu, J. Botina, and H. Rabitz, J. Chem. Phys. 108, 1953 (1998);
 (19) P. Gross, H. Singh, and H. Rabitz, Phys. Rev. A 47, 4593 (1993).
 (20) W. Zhu and H. Rabitz, J. Chem. Phys. 119, 3619 (2003).
 (21) S. Mancini, V.I. Manko, and H.M. Wiseman, J. Opt. B: Quantum Semiclass. Opt. 7, S177 (2005).
 (22) S. Chu, Nature 416, 206 (2002)
 (23) E. Runge and E.K.U. Gross, Phys. Rev. Lett. 52, 997 (1984).
 (24) C. Leforestier et al., J. Comp. Phys. 94, 59 (1991).
 (25) P. Elliott and N. Maitra, Phys. Rev. A 85, 052510 (2012)
 (26) N.T. Maitra, K. Burke, and C. Woodward, Phys. Rev. Lett. 89, 023002 (2002).
 (27) N.T. Maitra, Lect. Notes Phys. 837, 167 (Springer, Berlin, Heidelberg, 2012).