Investigation of optimal control problems governed by a time-dependent Kohn-Sham model

TL;DR

This paper explores optimal control problems within the time-dependent Kohn-Sham model of TDDFT, proposing a numerical framework validated by experiments to control multi-electron quantum systems efficiently.

Contribution

It develops a theoretical and numerical framework for optimal control of Kohn-Sham TDDFT models, including existence proofs and a computational scheme.

Findings

Numerical experiments demonstrate the effectiveness of the control approach.

The proposed discretization and conjugate gradient method successfully solve control problems.

The framework can handle various control objectives in quantum systems.

Abstract

Many application models in quantum physics and chemistry require to control multi-electron systems to achieve a desired target configuration. This challenging task appears possible in the framework of time-dependent density functional theory (TDDFT) that allows to describe these systems while avoiding the high dimensionality resulting from the multi-particle Schr\"{o}dinger equation. For this purpose, the theory and numerical solution of optimal control problems governed by a Kohn-Sham TDDFT model are investigated, considering different objectives and a bilinear control mechanism. Existence of optimal control solutions and their characterization as solutions to Kohn-Sham TDDFT optimality systems are discussed. To validate this control framework, a time-splitting discretization of the optimality systems and a nonlinear conjugate gradient scheme are implemented. Results of numerical experiments demonstrate the computational capability of the proposed control approach.

Equations

$$i\frac{\partial }{\partial t}\Psi (x_1,\dots x_N,t)=H\Psi (x_1,\dots x_N,t),$$

$$H=\sum _{i=1}^N\left(-{\nabla }_i^2+V_{ext}(x_i,t)\right)+\sum _{i<j}\frac{1}{|x_i-x_j|},$$

$$\Psi (x_1,x_2,\dots ,x_N,t)$$

$$\rho (x,t)=N\int |\Psi (x,\dots x_N,t){|}^2\mathrm{d}{x}_2\cdots \mathrm{d}{x}_N.$$

$$\rho (x,t)=\sum _{i=1}^N|{\psi }_i(x,t){|}^2,$$

$$i\frac{\partial {\psi }_j}{\partial t}(x,t)$$

$${\psi }_j(x,t)$$

$$V_{Hxc}(x,\rho (x,t))=V_H(x,\rho (x,t))+V_x(\rho (x,t))+V_c(\rho (x,t)).$$

$$V_H(x,\rho )={\int }_{\Omega }\frac{\rho (y)}{|x-y|}\mathrm{d}y.$$

$$V_x^{2D}(\rho )=-\sqrt{\frac{8}{\pi }}\sqrt{\rho },$$

$$V_x:[0,\infty )\to [0,p(2R)],$$

$$p(\rho )=\frac{(n+1)R^{\frac{1}{n}-4}}{4n^2}{\rho }^4-\frac{(4n+5)R^{\frac{1}{n}-3}}{3n^2}{\rho }^3+\frac{2(n+2)R^{\frac{1}{n}-2}}{n^2}{\rho }^2-\frac{4R^{\frac{1}{n}-1}}{n^2}\rho +\frac{(12n^2-11n+17)R^{\frac{1}{n}}}{12n^2},$$

$$V_{ext}(x,t,u)=u(t)V_u(x)+V_0(x).$$

$${\int }_{\Omega }|V_x(\Psi (x,t))\Psi (x,t)-V_x({\rm Y}(x,t)){\rm Y}(x,t){|}^2\mathrm{d}x\le {\int }_{\Omega }{L}^2|\Psi (x,t)-{\rm Y}(x,t){|}^2\mathrm{d}x.$$

$$\Vert V_x(\Psi )\Psi -V_x({\rm Y}){\rm Y}{\Vert }_{X^{\ast }}^2\le \Vert V_x(\Psi )\Psi -V_x({\rm Y}){\rm Y}{\Vert }_Y^2$$

$$={\int }_0^T\Vert V_x(\Psi (x,t))\Psi (x,t)-V_x({\rm Y}(x,t)){\rm Y}(x,t){\Vert }_{L^2(\Omega ;{\mathbb{C}}^N)}^2\mathrm{d}t$$

$$\le {\int }_0^T{L}^2\Vert \Psi (x,t)-{\rm Y}(x,t){\Vert }_{L^2(\Omega ;{\mathbb{C}}^N)}^2\mathrm{d}t=L^2\Vert \Psi -{\rm Y}{\Vert }_Y^2\le L^2\Vert \Psi -{\rm Y}{\Vert }_X^2,$$

$$\begin{array}{rl}i{\left({\partial }_t\Psi (t),\Phi \right)}_{L^2(\Omega ;{\mathbb{C}}^N)}=& {\left(\nabla \Psi ,\nabla \Phi \right)}_{L^2(\Omega ;{\mathbb{C}}^N)}+{\left(V_{ext}(\cdot ,t,u)\Psi ,\Phi \right)}_{L^2(\Omega ;{\mathbb{C}}^N)}\\ & +{\left(V_{Hxc}(\Psi (t))\Psi (t),\Phi \right)}_{L^2(\Omega ;{\mathbb{C}}^N)},\end{array}$$

$$\Psi \in L^2(0,T;H^2(\Omega ;{\mathbb{C}}^N))\cap L^{\infty }(0,T;H_0^1(\Omega ;{\mathbb{C}}^N));$$

$$\begin{array}{rl}J(\Psi ,u)& =\underset{J_{\beta }}{\underbrace{\frac{\beta }{2}{\int }_0^T{\int }_{\Omega }(\rho (x,t)-{\rho }_d(x,t){)}^2\mathrm{d}x\mathrm{d}t}}+\underset{J_{\eta }}{\underbrace{\frac{\eta }{2}{\int }_{\Omega }{\chi }_A(x)\rho (x,T)\mathrm{d}x}}+\underset{J_{\nu }}{\underbrace{\frac{\nu }{2}\Vert u{\Vert }_{H^1(0,T;\mathbb{R})}^2}}.\end{array}$$

$$\underset{(\Psi ,u)\in (W,H^1(0,T;\mathbb{R}))}{\min }J(\Psi ,u)\text{ subject to }\text{eqref}eq:KSweak.$$

$$L(\Psi ,u,\Lambda )=J(\Psi ,u)+L_1(\Psi ,u,\Lambda ),\text{mbox}where$$

$$L_1=\mathrm{Re}\left(\sum _{j=1}^N{\int }_0^T{\int }_{\Omega }\left(i\frac{\partial {\psi }_j(x,t)}{\partial t}-\left(-{\nabla }^2+V_{ext}(x,t,u)+V_{Hxc}(x,t,\rho )\right){\psi }_j(x,t)\right)\overline{{\lambda }_j(x,t)}\mathrm{d}x\mathrm{d}t\right).$$

$$i\frac{\partial {\psi }_m(x,t)}{\partial t}$$

$${\psi }_m(x,0)$$

$$i\frac{\partial {\lambda }_m}{\partial t}$$

$${\lambda }_m(x,T)$$

$$\nu u(t)+\mu (t)=0,$$

$${\nabla }_{\psi }{L}_H:$$

$${\nabla }_{\psi }{L}_{xc}:$$

$$\left(-\frac{{\mathrm{d}}^2}{\mathrm{d}{t}^2}+1\right)\mu =-\mathrm{Re}{\left(\Lambda ,V_u\Psi \right)}_{L^2(\Omega ;{\mathbb{C}}^N)},\mu (0)=0,\mu (T)=0,$$

$$\hat{J}:H^1(0,T;\mathbb{R})\to \mathbb{R},\hat{J}(u):=J(\Psi (u),u),$$

$$\hat{\Psi }={\left(\begin{array}{c}\mathrm{Re}{\psi }_1,\mathrm{Im}{\psi }_1,\cdots ,\mathrm{Re}{\psi }_N,\mathrm{Im}{\psi }_N\end{array}\right)}^T$$

$$c(\Psi ,u):=\tilde{c}(\Psi ,u)-V_{Hxc}(\Psi )\Psi ,\text{ where }\tilde{c}(\Psi ,u):=i\frac{\partial \Psi }{\partial t}-\left(-{\nabla }^2+V_0+V_{u}u(t)\right)\Psi ,$$

$$\Vert B(\Psi ){\Vert }_{\mathcal{L}(W,X^{\ast })}$$

Full text

1 Introduction

Many models of interest in quantum physics and chemistry consist of multi-particle systems, that can be modelled by the multi-particle Schrödinger equation (SE). However, the space dimensionality of this equation increases linearly with the number of particles involved and the corresponding computational cost increases exponentially, thus making the use of the multi-particle SE prohibitive. This fact has motivated a great research effort towards an alternative formulation to the SE description that allows to compute the observables of a quantum multi-particle system using particle-density functions. This development, which starts with the Thomas-Fermi theory in 1927, reaches a decisive point with the works of Hohenberg, Kohn, and Sham [HK64, KS65] that propose an appropriate way to replace a system of $N$ interacting particles in an external potential $V_{ext}$ by another system of non-interacting particles with an external potential $V_{ext}+V_{Hxc}$, such that the two models have the same electronic density. These works and many following ones focused on the computation of stationary (ground) states and obtained successful results that motivated the extension of the density function theory (DFT) theory to include time-dependent phenomena. This extension was first proposed by Runge and Gross in [RG84], and further investigated from a mathematical point of view in the work of van Leeuwen [vL99]. We refer to [ED11] for a modern introduction to DFT and to [MUN*+*06] for an introduction to time-dependent DFT (TDDFT).

Similar to the stationary case, the Runge–Gross theorem proves that, given an initial wavefunction configuration, there exists a one-to-one mapping between the potential in which the system evolves and the density of the system. Therefore, under appropriate assumptions, given a SE for a system of interacting particles in an external potential, there exists another SE model, unique up to a purely time-dependent function in the potential [vL99], of a non-interacting system with an augmented potential whose solution provides the same density as the solution to the original SE problem. We refer to this TDDFT model as the time-dependent Kohn-Sham (TDKS) equation. Notice that the external potential $V_{ext}$ modelling the interaction of the particles (in particular, electrons) with an external (electric) field enters without modification in both the multi-particle SE model and the TDKS model.

This latter fact is important in the design of control strategies for multi-particle quantum systems because control functions usually enter in the SE model as external time-varying potentials. Therefore control mechanisms can be determined in the TDDFT framework that are valid for the original multi-particle SE system. Recently, various quantum mechanical optimal control problems governed by the SE have been studied in the literature, see for example [vWB08], [vWBV10] and [MST06]. Moreover, quantum control problems governed by TDDFT models have already been investigated; see, e.g., [CWG12] and have been implemented in TDDFT codes as the well-known Octopus [CAO*+*06]. However, the available optimization schemes are mainly based on less competitive Krotov’s method and consider only finite-dimensional parameterized controls. Furthermore, much less is known on the theory of the TDDFT optimal control framework and on the use and analysis of more efficient optimization schemes that allow to compute control functions belonging to a much larger function space.

We remark that the functional analysis of optimization problems governed by the TDKS equations and the investigation of optimization schemes requires the mathematical foundation of the governing model. At the best of our knowledge, only few contributions addressing this issue are available; we refer to [RPvL15, Jer15, SCB17] for results concerning the existence and uniqueness of solutions to the TDKS equations.

This work contributes to the field of optimal control theory for multi-particle quantum systems presenting a theoretical analysis of optimal control problems governed by the TDKS equation. To validate the proposed framework, we implement an efficient approximation and optimization scheme for these problems.

This paper is organized as follows. In Section 2, we illustrate multi-particle SE models and the Kohn-Sham (KS) approach to TDDFT. Since these models are less known in the PDE optimal control community, and the literature on these problems is sparse, we make a special effort to provide a detailed presentation and to present results that are instrumental for the discussion that follows. In Section 3, we state a class of optimal control problems and discuss the related first-order optimality systems. The details of the derivation of the optimality system are elaborated in the Appendix A. The analysis of the optimal control problems is presented in Section 4. We show existence of optimal solutions to the control problems and prove necessary optimality conditions. Section 5 is dedicated to suitable approximation and numerical optimization schemes are discussed. We consider time-splitting schemes [BJM02, FOS15] and discuss their accuracy properties. To solve the optimality systems, we implement a nonlinear conjugate gradient scheme. In Section 6, results of numerical experiments are presented that demonstrate the effectiveness of the proposed control framework. A section of conclusions completes this work.

2 The TDKS model

In the Schrödinger quantum mechanics framework, the state of a $N$ electrons system is described by a wave function $\Psi$, whose time evolution is governed by the following Schrödinger equation (SE)

[TABLE]

where $x_{1},\dotsc,x_{N}\in\mathbb{R}^{n}$ are the position vectors of the $N$ particles. We use atomic units, i.e. $\hbar=4\pi\epsilon_{0}=1$, $m_{e}=\frac{1}{2}$.

The Hamiltonian $H$ consists of a kinetic term, the Coulomb interaction between the charged particles (electrons), $\sum_{i<j}\frac{1}{|x_{i}-x_{j}|}$, and an external potential, $V_{ext}$. We have

[TABLE]

where $|\cdot|$ is the Euclidean norm and $\nabla_{i}$ is the $n$-dimensional vector gradient with respect to $x_{i}$.

The Pauli principle states that the wave function of a system of electrons has to be antisymmetric with respect to the exchange of two coordinates. For this purpose, given $N$ orthogonal single particle wave functions (orbitals), $\psi_{j}(x,t)$, that correspond to the Hamiltonian of the $j$th particle, $H_{j}=\left(-\nabla_{j}^{2}+V_{ext}(x_{j},t)\right)$ one can build the following antisymmetric wave function

[TABLE]

This is called Slater determinant. This wave function solves (2.1) if the particles do not interact. However, in the presence of an interaction potential, the solution to (2.1) will be given by an infinite sum of Slater determinants.

This fact shows that the effort of solving a multi-particle SE increases exponentially with the number of particles. To avoid this curse of dimensionality, the approach of DFT is to consider, instead of the wave function on a $nN$-dimensional space, the corresponding electronic density defined of the physical space of $n$-dimensions given by

[TABLE]

For a Slater determinant, one finds that the corresponding density is as follows

[TABLE]

where $\psi_{i}$ represents the wave function of the $i$th particle.

The DFT approach of Kohn and Sham [KS65] to model multi-particle problems was to replace the system of $N$ interacting particles, subject to an external potential $V_{ext}$, by another system of non-interacting particles subject to an augmented potential $V_{ext}+V_{Hxc}$, such that the two models provide the same density. Van Leeuwen [vL99] proved that such a system exists under appropriate conditions on the potentials and on the resulting densities. In particular, for this proof it is required that the wave function is twice continuously differentiable in space and analytic in time and the potential has to have finite expectation values and be differentiable in space and analytic in time.

Based on this development, we consider the time-dependent Kohn-Sham system given by

[TABLE]

It is not trivial to get an initial condition that appropriately represents the interacting system because it may not have a Slater determinant as starting wavefunction. A common choice is to solve the ground state DFT problem and take the eigenstates corresponding to the lowest $N$ eigenvalues as $\psi_{j}^{0}$.

Notice that an TDKS system is formulated in $n$ spatial dimensions and consists of $N$ coupled Schrödinger equations. With the appropriate choice of $V_{Hxc}$, which contains the coupling through the dependence on $\rho$, the solution to (2.5)–(2.6) provides the correct density of the original system, so that all observables, which can be formulated in terms of the density, can be determined by this method.

The main challenge of the DFT framework is to construct KS potentials that encapsulate all the multi-body physics. One class of approximations is called the local density approximation (LDA) [ED11], because in this approach, the KS potential at some point $x$ only depends on the value of the density $\rho(x)$ at this specific point. We use the adiabatic LDA, which means that LDA is applied at every time separately such that $V_{Hxc}(x,t)=V_{Hxc}(\rho(x,t))$. Notice that, if one allows $V_{Hxc}(t)$ to depend on the whole history $V_{Hxc}(\tau)$, $0\leq\tau\leq t$, the resulting adjoint equation would be an integro-differential equation, which would be much more involved to solve.

We remark that the Kohn-Sham potential $V_{Hxc}$ is usually split into three terms, that is, the Hartree potential, and the exchange and correlation potentials as follows

[TABLE]

The Hartree potential $V_{H}$ models electrons interaction due to the Coulomb force. This term dominates $V_{Hxc}$, while the other two parts represent quantum mechanical corrections. Later, we refer to the exchange-correlation part of the potential also as $V_{xc}=V_{x}+V_{c}$. Notice that the classical electric field is given by $\phi(x)=\int_{\Omega}\frac{\rho_{c}(y)}{|x-y|}\mathrm{d}y$, where $x$ and $y$ are positions in space and $\rho_{c}$ is the charge density. Since we are using natural units, charge and particle densities become the same and we have the following

[TABLE]

In quantum mechanics, the Pauli exclusion principle states that two electrons cannot share the same quantum state. This results in a repulsive force between the electrons. In DFT this feature is modeled by introducing two terms: the exchange and the correlation potentials. The exchange potential $V_{x}$ contains the Pauli principle for a homogeneous electron gas. In the LDA framework, it can be calculated explicitly as follows; see, e.g., [Con08, PY89]

[TABLE]

However, quantum mechanics is not applicable for very short distances such that relativistic effects start to play a role. Therefore, we introduce a cut-off of the potential at unphysically large densities, while preserving all required properties in the range of validity of the DFT framework.

We define the exchange potential as follows

[TABLE]

where

[TABLE]

with $\alpha_{2}=-\sqrt{\frac{8}{\pi}}$ and $\alpha_{3}=-\sqrt[3]{\frac{3}{\pi}}$ and $R$ sufficiently large; e.g., $R\approx(10^{57}\frac{1}{\textrm{m}})^{3}N$. This potential is twice continuously differentiable and globally bounded.

The remaining part of the interaction is called the correlation potential $V_{c}$. No analytic expression is known for it. However, it is possible to determine the shape of this potential as a function of the density using Quantum Monte Carlo methods; see, e.g., [AMGGB02]. In Figure 1, we plot the shape of $V_{c}$ used in the numerical experiments in Section 6.

All correlation potentials commonly used are of similar structure; see, e.g., [MOB12]. They are zero for zero density and otherwise negative, while having a convex shape. Furthermore, they are bounded by $V_{x}$ in the sense that $|V_{c}(\rho)|<|V_{x}(\rho)|$ for all $\rho\in\mathbb{R}^{+}$.

It is clear that in applications, confined electron systems subject to external control are of paramount importance. The confinement is obtained considering external potentials such that $\Psi$ is non-zero only on a bounded domain $\Omega$. For this reason, we denote by $V_{0}$ a confining potential that may represent the attracting potential of the nuclei of a molecule or the walls of a quantum dot. A typical model is the harmonic oscillator potential, $V_{0}(x)=u_{0}\,x^{2}$.

A control potential aims at steering the quantum system to change its configuration towards a target state or to optimize the value of a given observable. In most cases, this results in a change of energy that necessarily requires a time-dependent interaction of the electrons with an external electro-magnetic force. For this purpose, we introduce a control potential with the following structure $V_{ext}^{c}(x,t)=u(t)\,V_{u}(x)$, where $u(t)$ has the role of a modulating amplitude. A specific case is the dipole control potential, $V_{ext}^{c}(x,t)=u(t)\,x$.

In our the TDKS system, we consider the following external potential

[TABLE]

In particular, we consider the control of a quantum dot by a changing gate voltage modeled by a variable quadratic potential, $V_{u}(x)=x^{2}$, and a laser control in dipole approximation, $V_{u}(x)=x\cdot p$, with a polarization vector $p$.

With this setting, we have completely specified our TDKS model. Next, we discuss the corresponding functional analytic framework.

We consider our TDKS model defined on a bounded domain $\Omega\subset\mathbb{R}^{n}$, $n=2,3$, with $\partial\Omega\in C^{2}$ and $t\in(0,T)$, together with initial conditions $\psi_{j}(\cdot,0)=\psi_{j}^{0}(\cdot)\in H^{1}_{0}(\Omega;\mathbb{C})$, $j=1,\ldots,N$. For brevity, we denote our KS wavefunction by $\Psi=(\psi_{1},\ldots,\psi_{N})$, and $\Psi^{0}=(\psi_{1}^{0},\ldots,\psi_{N}^{0})$. Therefore, we have $|\Psi(x,t)|^{2}=\rho(x,t)$.

We define the following function spaces. $X=L^{2}(0,T;H^{1}_{0}(\Omega;\mathbb{C}^{N}))$ and $W=\{\Psi\in X|\Psi^{\prime}\in X^{*}\}$ with the norms $\|\Psi\|_{X}^{2}=\int_{0}^{T}\|\Psi\|_{L^{2}}^{2}+\|\nabla\Psi\|_{L^{2}}^{2}\mathrm{d}t$ and $\|\Psi\|_{W}^{2}=\|\Psi\|_{X}^{2}+\|\Psi^{\prime}\|_{X^{*}}^{2}$, where $X^{*}$ is the dual space of $X$; we also need $Y=L^{2}(0,T;L^{2}(\Omega;\mathbb{C}^{N}))$ which is endowed with the usual norm.

To improve readability of the analysis that follows, we write the potentials in (2.7) as functions of $\Psi$ instead of $\rho=|\Psi|^{2}$. We shall also omit the explicit dependence of $V_{H}$ on $x$ if no confusion may arise.

Lemma 1**.**

The exchange potential term $V_{x}(\Psi)\Psi$ is Lipschitz continuous, i.e. $\|V_{x}(\Psi)\Psi-V_{x}(\Upsilon)\Upsilon\|_{L^{2}(\Omega;\mathbb{C}^{N})}\leq L\|\Psi-\Upsilon\|_{L^{2}(\Omega;\mathbb{C}^{N})}$ and $\|V_{x}(\Psi)\Psi-V_{x}(\Upsilon)\Upsilon\|_{X^{*}}\leq L\|\Psi-\Upsilon\|_{X}$.

Proof.

The function $f:\mathbb{C}^{N}\rightarrow\mathbb{C}^{N}$, $f(z)=V_{x}(z)z$ is continuously differentiable with bounded derivative, hence Lipschitz continuous with Lipschitz constant $L$ from $\mathbb{C}^{N}$ to $\mathbb{C}^{N}$. With this preparation, we have Lipschitz continuity from $L^{2}$ to $L^{2}$ as follows

[TABLE]

Similarly, we have Lipschitz continuity from $X$ to $X^{*}$ as follows

[TABLE]

where we use the Gelfand triple $X\hookrightarrow Y\hookrightarrow X^{*}$ and the fact that $V_{x}(\Psi)\Psi\in Y$ as $V_{x}(\Psi)\in L^{\infty}(0,T;L^{\infty}(\Omega;\mathbb{R}))$. ∎

Throughout the paper, we make the following assumptions on the potentials.

Assumption**.**

The correlation potential $V_{c}$ is uniformly bounded in the sense that$|V_{c}(\Psi(x,t))|\leq K,$* $\forall x\in\Omega$, $t\in[0,T]$, $\Psi\in Y$; this corresponds to the correlation potentials from the Libxc library [MOB12] applying a similar approach to (2.9).* 2. 2.

The correlation potential term $V_{c}(\Psi)\Psi$ is continuously real-Fréchet differentiable (see Definition 5) with bounded derivative, c.f. Lemma 9 for the same result on the exchange potential. 3. 3.

The confining potential and the spacial dependence of the control potential are bounded, i.e. $V_{0},V_{u}\in L^{\infty}(\Omega;\mathbb{R})$; as we consider a finite domain, this is equivalent to excluding divergent external potentials. 4. 4.

The control is $u\in H^{1}(0,T;\mathbb{R})$. This is a classical assumption in optimal control, see, e.g. **[vWB08]**.

The following theorem from [SCB17] states existence and uniqueness of (2.5)–(2.6) in a setting well-suited for optimal control.

Theorem 2**.**

The weak formulation of (2.5), with $V_{ext}\in H^{1}(0,T;L^{\infty}(\Omega;\mathbb{R}))$, $\Psi(0)=\Psi^{0}\in L^{2}(\Omega;\mathbb{C}^{N})$, admits a unique solution in $W$, that is, there exists $\Psi\in X$, $\Psi^{\prime}\in X^{*}$, such that

[TABLE]

for all $\Phi\in H_{0}^{1}(\Omega;\mathbb{C}^{N})$ and a.e. in $(0,T)$.

If $\Psi^{0}\in H^{1}_{0}(\Omega;\mathbb{C}^{N})$ and $\partial\Omega\in C^{2}$, then the unique solution to (2.11) is as follows

[TABLE]

if in addition $\Psi^{0}\in H^{2}(\Omega;\mathbb{C}^{N})$, we have $\Psi\in L^{\infty}(0,T;H^{2}(\Omega;\mathbb{C}^{N})\cap H_{0}^{1}(\Omega;\mathbb{C}^{N}))$.

By the continuous embedding $W\hookrightarrow C([0,T];L^{2}(\Omega;\mathbb{C}^{N}))$, see e.g. [Eva10, p. 287], the solution is continuous in time.

Similar problems have been studied in [Jer15]:

Theorem 3**.**

Assuming that $V_{ext}\in C^{1}([0,T];C^{1}(\overline{\Omega};\mathbb{R}))$ and $V_{ext}\geq 0$, and a Lipschitz condition on $V_{xc}$ and a continuity assumption on $V_{xc}$, then (2.11) with $\Psi^{0}\in H_{0}^{1}(\Omega;\mathbb{C}^{N})$ has a unique solution in $C([0,T];H^{1}_{0}(\Omega;\mathbb{C}^{N}))\cap C^{1}([0,T];H^{-1}(\Omega;\mathbb{C}^{N}))$.

Taking into account only the Hartree potential but not the exchange-correlation potential, existence of a unique solution in $C([0,\infty);H^{2}(\mathbb{R}^{3};\mathbb{C}))\cap C^{1}([0,\infty);L^{2}(\mathbb{R}^{3};\mathbb{C}))$ is shown in [CL99].

As the potential in (2.11) is purely real, the norm of the wave function is conserved, see e.g. [SCB17]. We have

Lemma 4**.**

The $L^{2}(\Omega;\mathbb{C})$-norm of the solution to (2.11) is conserved in the sense that$\|\Psi(\cdot,t)\|_{L^{2}(\Omega;\mathbb{C}^{N})}=\|\Psi^{0}\|_{L^{2}(\Omega;\mathbb{C}^{N})}$ for all $t\in[0,T]$.

3 Formulation of TDKS optimal control problems

Optimal control of quantum systems is of fundamental importance in quantum mechanics applications. The objectives of the control may be of different nature ranging from the breaking of a chemical bond in a molecule by an optimally shaped laser pulse to the manipulation of electrons in two-dimensional quantum dots by a gate voltage potential. In this framework, the objective of the control is modeled by a cost functional to be optimized under the differential constraints represented by the quantum model (in our case the TDKS equation) including the control mechanism.

We consider an objective $J$ that includes different target functionals and a control cost as follows

[TABLE]

The first term, $J_{\beta}$ models the requirement that the electron density evolves following as close as possible a given target trajectory, $\rho_{d}$. We remark that $J_{\beta}$ is only well-defined if $\Psi$ is at least in $L^{4}(0,T;L^{4}(\Omega;\mathbb{C}^{N}))$. This is guaranteed by the improved regularity from Theorem 2 for $\Psi^{0}\in H^{2}(\Omega;\mathbb{C}^{N})$. The term $J_{\eta}$ aims at locating the density outside of a certain region $A$. The term $J_{\nu}$ penalizes the cost of the control. We remark that the regularization term $J_{\nu}$ can be any weighed $H^{1}(0,T;\mathbb{R})$ norm, e.g. $\|u^{\prime}\|_{L^{2}(0,T;\mathbb{R})}^{2}+a\|u\|_{L^{2}(0,T;\mathbb{R})}^{2}$ for $a>0$. We assume that the target weights are all non-negative $\beta,\eta\geq 0$, with $\beta+\eta>0$, and the regularization weight $\nu>0$. The characteristic function of $A$ is given by $\chi_{A}(x)=\begin{cases}1&x\in A\\ 0&\text{otherwise}\end{cases}$.

Our purpose is to find an optimal control function $u$, which modulates a dipole or a quadratic potential, such that $J(\Psi,u)$ is minimized subject to the constraint that $\Psi$ satisfies the TDKS equations. This problem is formulated as follows

[TABLE]

Here and in the following, we assume that $\Psi^{0}\in H^{1}(\Omega;\mathbb{C}^{N})$.

The solutions to this PDE-constrained optimization problem are characterized as solutions to the corresponding first-order optimality conditions [BS12, Trö10]. These conditions for (3.2) can be formally obtained by setting to zero the gradient of the following Lagrange function

[TABLE]

The function $\Lambda=(\lambda_{1},\ldots,\lambda_{N})$, where $\lambda_{j}\in L^{2}(0,T;H^{1}_{0}(\Omega;\mathbb{C}))$, $j=1,\ldots,N$, represent the adjoint variables. In the Lagrange formalism, a solution to (3.2) corresponds to a stationary point of $L$, where the derivatives of $L$ with respect to $\psi_{j}$, $\lambda_{j}$ and $u$ must be zero along any directions $\delta\psi$, $\delta\lambda$, and $\delta u$. A detailed calculation of these derivatives can be found in Appendix A. The main difficulty in the derivation is the complex and non-analytic dependence of the Kohn-Sham potential on the wave function, which results in the terms $\nabla_{\psi}L_{H}$, $\nabla_{\psi}L_{xc}$ in (3.4c) below.

The first-order optimality conditions define the following optimality system

[TABLE]

where $m=1,\ldots,N$, and

[TABLE]

Further, $\mu$ is the $H^{1}$-Riesz representative of the continuous linear functional$\left(-\operatorname{Re}\left(\Lambda,\,V_{u}\Psi\right)_{L^{2}(\Omega;\mathbb{C}^{N})},\,\cdot\right)_{L^{2}(0,T;\mathbb{R})}$; see Theorem 18 below. Assuming that $u\in H^{1}_{0}(0,T;\mathbb{R})$, $\mu$ can be computed by solving the equation

[TABLE]

which is understood in a weak sense. For more details, see, e.g., [vWB08].

Theorem 2 guarantees that (3.4a)–(3.4b) is uniquely solvable and hence the reduced cost functional

[TABLE]

is well defined.

To ensure the correct regularity properties of the adjoint variables, we do not use the Lagrange formalism explicitly in this work. Instead, we directly use the existence theorem from [SCB17] for the adjoint equation (3.4c)–(3.4d), an extension of Theorem 2, and derive the gradient of the reduced cost functional $\hat{J}$ from these solutions in Theorem 18 below. Later, we use this gradient to construct a numerical optimization scheme to minimize $\hat{J}$.

4 Theoretical analysis of TDKS optimal control problems

In this section, we present a mathematical analysis of the optimal control problem (3.2). To this end, we first show that both the constraint given by the TDKS equations and the cost functional $J$ are continuously real-Fréchet differentiable. Subsequently, we prove existence of solutions to the optimization problem.

The Kohn-Sham potential depends on the density $\rho$. The density is a real-valued function of the complex wavefunction and can therefore not be holomorphic. As complex differentiability is a stronger property than what we need in the following, we introduce the following weaker notion of real-differentiability.

Definition 5**.**

Let $X,Y$ be complex Banach spaces. A map $f:X\rightarrow Y$ is called real-linear if and only if

$f(x)+f(y)=f(x+y)$* $\forall x,y\in X$ and* 2. 2.

$f(\alpha x)=\alpha f(x)$* $\forall\alpha\in\mathbb{R}$ and $\forall x\in X$.*

The space of real-linear maps from $X$ to $Y$ is a Banach space.

We call a map $f:X\rightarrow Y$ real-Gâteaux (real-Fréchet) differentiable if the standard definition of Gâteaux (Fréchet) differentiability holds for a real-linear derivative operator.

Remark**.**

In complex spaces, the notion of real-Gâteaux (real-Fréchet) differentiability is weaker than Gâteaux (Fréchet) differentiability. However, all theorems for differentiable functions in $\mathbb{R}^{2}$ also hold in $\mathbb{C}$ for functions that are just real-differentiable. This is the case for all theorems that we will make use of, e.g. the chain rule and the implicit function theorem. Therefore, it is enough to show real-Fréchet differentiability of the constraint.

An alternative and fully equivalent approach is to consider a real vector

[TABLE]

and the corresponding matrix Schrödinger equation.

Theorem 6**.**

The map $c:W\times H^{1}(0,T;\mathbb{R})\rightarrow X^{*}$, defined as

[TABLE]

is continuously real-Fréchet differentiable.

Notice that $c(\Psi,u)=0$ represents the TDKS equation.

Remark**.**

We remark that $V_{x}(\Psi)\Psi$, $V_{H}(\Psi)\Psi$, and $V_{u}u\Psi$ are in $Y$. Hence, the operator norm of their derivatives in $\mathcal{L}(W,X^{*})$ can be bounded in $\mathcal{L}(W,Y)$ as follows.

[TABLE]

where $B$ denotes the derivative of $V_{x}(\Psi)\Psi$, $V_{H}(\Psi)\Psi$, or $V_{u}u\Psi$, respectively.

Further, to prove Theorem 6, we need the following lemmas. We begin with studying the nonlinear exchange potential term.

Lemma 7**.**

The map $W\ni\Psi\mapsto V_{x}(\Psi)\Psi\in Y$ is real-Gâteaux differentiable for all $\Psi\in W$ and its real-Gâteaux derivative, denoted with $A(\Psi)$, is specified as follows

[TABLE]

Proof.

First, we need the derivative of the density $\rho$ which is a non-holomorphic function of a complex variable. Differentiation of $\rho$ can be done using the Wirtinger calculus [Rem91] where $\Psi$ and $\overline{\Psi}$ are treated as independent variables. By using these calculus rules, the real-Fréchet derivatives of $\rho=\sum_{m=1}^{M}\psi_{m}\overline{\psi_{m}}=\left(\Psi,\,\Psi\right)_{\mathbb{C}}$ are given by

[TABLE]

The directional derivative of $V_{x}(\Psi)\Psi$ along $\delta\Psi$ is given by $A(\Psi):W\rightarrow X^{*}$,

[TABLE]

Using the definition of $V_{x}$, we have

[TABLE]

and $\delta\Psi\mapsto A(\Psi)(\delta\Psi)$ is obviously linear in $\delta\Psi$ over the real scalars. We are left to show that $A(\Psi)$ is a bounded operator.

[TABLE]

For the first term, we use $V_{x}(\Psi)\in L^{\infty}(0,T;L^{\infty}(\Omega;\mathbb{R}))$ as well as $\|\Psi\|_{Y}\leq\|\Psi\|_{W}$ to obtain

[TABLE]

For the second term, we decompose the domain into $\Omega=\Omega_{1}\cup\Omega_{2}\cup\Omega_{3}$ depending on the size of $\rho$: $\Omega_{1}:=\{x\in\Omega|\,|\Psi|^{2}\leq R\}$, $\Omega_{2}:=\{x\in\Omega|\,R<|\Psi|^{2}<2R\}$, and $\Omega_{3}:=\{x\in\Omega|\,|\Psi|^{2}\geq 2R\}$. Using the fact that $|\Psi|^{2}=\rho$ is bounded by $R$ in $\Omega_{1}$ and by $2R$ in $\Omega_{2}$ as well as that $A_{2}(\Psi)=0$ in $\Omega_{3}$, and that $\frac{\partial p}{\partial\rho}$ is monotonically decreasing between $R$ and $2R$ we obtain

[TABLE]

We have shown that the directional derivative $\delta\Psi\mapsto A(\Psi)(\delta\Psi)$ is a bounded linear map for all $\Psi\in W$, hence $\Psi\mapsto V_{x}(\Psi)\Psi$ is real-Gâteaux differentiable. ∎

Before improving this result to real-Fréchet differentiability, we prove a general result on the $L^{2}$-norm of products of functions. To this end, denote $Y_{1}=L^{2}(0,T;L^{2}(\Omega;\mathbb{C}))$.

Lemma 8**.**

Let $k\in\mathbb{N}$ be arbitrary. Given two functions $f\in L^{\infty}(\Omega;\mathbb{C})$ and $g\in L^{2}(\Omega;\mathbb{C}^{k})$. Then $\|fg\|_{L^{2}(\Omega;\mathbb{C}^{k})}\leq\frac{\sqrt{\mu(\Omega)}}{\sqrt{2\pi}^{n}}\|f\|_{L^{2}(\Omega;\mathbb{C})}\|g\|_{L^{2}(\Omega;\mathbb{C}^{k})}$. Similarly, for $f\in L^{\infty}(0,T;L^{\infty}(\Omega;\mathbb{C}))$ and $g\in Y$, we have$\|fg\|_{Y}\leq\frac{\sqrt{\mu(\Omega)T}}{\sqrt{2\pi}^{n+1}}\|f\|_{Y_{1}}\|g\|_{Y}$.

Proof.

The Fourier transform $\mathcal{F}:L^{2}(\Omega;\mathbb{C})\rightarrow L^{2}(\Omega;\mathbb{C})$ is an isometry by Plancherel’s theorem, where we extend the functions by zero outside of $\Omega$. Furthermore, for $f\in L^{1}(\Omega;\mathbb{C})$, $g\in L^{1}(\Omega;\mathbb{C}^{k})$ the convolution theorem states $\sqrt{2\pi}^{n}\mathcal{F}\left[fg\right]=\mathcal{F}\left[f\right]\star\mathcal{F}\left[g\right]$. Hence, using the embeddings $L^{\infty}(\Omega;\mathbb{C})\hookrightarrow L^{2}(\Omega;\mathbb{C})\hookrightarrow L^{1}(\Omega;\mathbb{C})$, we find

[TABLE]

where Young’s inequality for convolutions was used, see e.g. [Sch05, Theorem 14.6]. The second claim is obtained by applying the same calculations on the domain $\Omega\times(0,T)$. ∎

Lemma 9**.**

The map $W\ni\Psi\mapsto V_{x}(\Psi)\Psi\in Y$ is continuously real-Fréchet differentiable with derivative $\operatorname{D}(V_{x}(\Psi)\Psi)=A(\Psi)$, where $A(\Psi)\in\mathcal{L}(W,Y)$ is given in Lemma 7.

Proof.

We prove that the real-Gâteaux derivative $A(\Psi)$ of $V_{x}(\Psi)\Psi$ at $\Psi$ is continuous from $W$ to $\mathcal{L}(W,Y)$. Then the real-Fréchet differentiability follows immediately from [AH11, Proposition A.3].

Once is proved that $V_{x}(\Psi)\Psi$ is real-Fréchet differentiable, the real-Gâteaux and the real-Fréchet derivatives coincide, and the continuity of the real-Gâteaux derivative carries over to the real-Fréchet derivative. Hence, we have to show the following

[TABLE]

For this purpose, to ease our discussion, we consider $\tilde{A}_{2}(\Psi)\psi_{j}\psi_{m}$, $j,m=1,\dotsc,N$, where

[TABLE]

such that \tilde{A}_{2}(\Psi)\psi_{m}\sum_{j=1}^{N}2\operatorname{Re}(\psi_{j}\overline{\delta\psi_{j}})=\bigl{(}A_{2}(\Psi)(\delta\Psi)\bigr{)}_{m}. For all $\epsilon>0$, $\Psi\mapsto\tilde{A}_{2}(\Psi)$ is continuously differentiable with bounded derivative for $|\Psi|\geq\epsilon$, hence the same holds for $\tilde{A}_{2}(\Psi)\psi_{j}\psi_{m}$, which is therefore Lipschitz continuous. In zero, $\tilde{A}_{2}(\Psi)\psi_{j}\psi_{m}$ is Hölder continuous with exponent $\alpha=\frac{2}{n}$,

[TABLE]

Together, we have Hölder continuity for all $\Psi\in W$, $|\tilde{A}_{2}(\Psi)\psi_{j}\psi_{m}-\tilde{A}_{2}(\Phi)\phi_{j}\phi_{m}|\leq c|\Psi-\Phi|_{\mathbb{C}^{N}}^{\alpha}$ and similarly $|\tilde{A}_{2}(\Psi)\overline{\psi_{j}}\psi_{m}-\tilde{A}_{2}(\Phi)\overline{\phi_{j}}\phi_{m}|\leq c|\Psi-\Phi|_{\mathbb{C}^{N}}^{\alpha}$.

We remark that, if a function $f:\mathbb{C}^{N}\rightarrow\mathbb{C}^{N}$ is Hölder continuous with exponent $\alpha$, i.e. $|f(x)-f(y)|<c|x-y|^{\alpha}$, then $f:L^{2}(\Omega;\mathbb{C}^{N})\rightarrow L^{2}(\Omega;\mathbb{C}^{N})$ is also Hölder continuous with the same exponent as follows

[TABLE]

As $\tilde{A}_{2}(\Psi)\psi_{j}\psi_{m}\in L^{\infty}(0,T;L^{\infty}(\Omega;\mathbb{C}))$ is bounded, we can apply Lemma 8. Now, we turn our attention to $A_{2}(\Psi)$ and use the Hölder continuity of $\tilde{A}_{2}(\Psi)\psi_{j}\psi_{m}$ to obtain the following estimate.

We have

[TABLE]

Furthermore, by Lemma 8 and the Hölder continuity of $V_{x}$, we have

[TABLE]

Now, we have the following

[TABLE]

This completes the proof of the real-Fréchet differentiability of $V_{x}(\Psi)\Psi$. ∎

Lemma 10**.**

The map $\Psi\mapsto V_{H}(\Psi)\Psi$ is continuously real-Fréchet differentiable from $W$ to $Y$ and from $W$ to $X^{*}$ with derivative $\operatorname{D}(V_{H}(\Psi)\Psi)(\delta\Psi)=V_{H}(\Psi)\delta\Psi+\int_{\Omega}\frac{2\operatorname{Re}\left(\Psi,\,\delta\Psi\right)_{\mathbb{C}}}{|x-y|}\mathrm{d}y\Psi$.

Proof.

By [SCB17, Lemma 2], $V_{H}(\Psi)\Psi\in Y$. The following expansion holds

[TABLE]

Hence, we get

[TABLE]

By proof of Lemma [SCB17, Lemma 2], the last three terms are bounded by$C_{N}\|\delta\Psi\|_{L^{2}(\Omega;\mathbb{C}^{N})}\|\Psi\|_{H^{1}(\Omega;\mathbb{C}^{N})}\|\delta\Psi\|_{L^{2}(\Omega;\mathbb{C}^{N})}$, $C_{N}\|\delta\Psi\|_{L^{2}(\Omega;\mathbb{C}^{N})}\|\delta\Psi\|_{H^{1}(\Omega;\mathbb{C}^{N})}\|\Psi\|_{L^{2}(\Omega;\mathbb{C}^{N})}$, respective $C_{N}\|\delta\Psi\|_{H^{1}(\Omega;\mathbb{C}^{N})}\|\delta\Psi\|_{L^{2}(\Omega;\mathbb{C}^{N})}^{2}$. Defining $\operatorname{D}(V_{H}(\Psi)\Psi)(\delta\Psi):=V_{H}(\Psi)\delta\Psi+\int_{\Omega}\frac{2\operatorname{Re}\left(\Psi,\,\delta\Psi\right)_{\mathbb{C}}}{|x-y|}\mathrm{d}y\Psi$, and using the embedding $W\hookrightarrow C([0,T];L^{2}(\Omega;\mathbb{C}^{N}))$, we obtain

[TABLE]

Hence, $V_{H}(\Psi)\Psi$ is real-Fréchet differentiable with derivative $V_{H}(\Psi)\delta\Psi+\int_{\Omega}\frac{2\operatorname{Re}\left(\Psi,\,\delta\Psi\right)_{\mathbb{C}}}{|x-y|}\mathrm{d}y\Psi$ and the derivative is continuous from $W$ to $\mathcal{L}(W,Y)$ and from $W$ to $\mathcal{L}(W,X^{*})$. ∎

Using these results, we can now prove Theorem 6 that states the real-Fréchet differentiability of the map $c$.

Proof of Theorem 6.

We prove that $c$ is a real-Fréchet differentiable function of $\Psi$ and $u$. For this purpose, we first consider the map $\tilde{c}$ defined in (4.1). Simple algebraic manipulation results in the following

[TABLE]

Next, we show that the Fréchet derivative of $\tilde{c}$ is given by

[TABLE]

To bound the reminder $\tilde{c}(\Psi+\delta\Psi,u+\delta u)-\tilde{c}(\Psi,u)-\operatorname{D}\tilde{c}(\Psi,u)(\delta\Psi,\delta u)$, we consider

[TABLE]

where we use the embedding $H^{1}(0,T;\mathbb{R})\hookrightarrow C[0,T]$. Further, using the estimate $\|\delta u\|_{C[0,T]}\leq c\|\delta u\|_{H^{1}(0,T;\mathbb{R})}$ and $\|\delta\Psi\|_{X}\leq\|\delta\Psi\|_{W}$, we can improve the inequality above in the following sense

[TABLE]

With this we can show the Fréchet differentiability as follows

[TABLE]

Hence, $\tilde{c}$ is Fréchet differentiable, and $\operatorname{D}\tilde{c}$ given in (4.3) represents its derivative. The derivative is continuous from $W\times H^{1}(0,T;\mathbb{R})$ to $X^{*}$.

The exchange potential is continuously real-Fréchet differentiable by Lemma 9, the correlation potential by Assumption 2, and the Hartree potential by Lemma 10. To summarize, we have the following real-Fréchet derivative of $c$.

[TABLE]

∎

We have discussed the differentiability properties of the differential constraint $c$ that are required below to prove the existence of a minimizer and for establishing the gradient of the reduced cost functional. Next, we study the cost functional $J$.

Theorem 11**.**

The cost functional $J:W\times H^{1}(0,T;\mathbb{R})\rightarrow\mathbb{R}$ defined in (3.1) is continuously real-Fréchet differentiable.

Proof.

The norm $\|u\|_{H^{1}(0,T;\mathbb{R})}^{2}$ is differentiable by standard results with $\operatorname{D}\left(\frac{\nu}{2}\|u\|_{H^{1}(0,T;\mathbb{R})}^{2}\right)(\delta u)=\nu\left(u,\,\delta u\right)_{H^{1}(0,T;\mathbb{R})}$. The tracking term $J_{\beta}=\frac{\beta}{2}\int_{0}^{T}\int_{\Omega}(\rho(x,t)-\rho_{d}(x,t))^{2}\mathrm{d}x\mathrm{d}t$ is a quadratic functional and hence Féchet differentiable with derivative

[TABLE]

$J_{\eta}=\frac{\eta}{2}\int_{\Omega}\chi_{A}(x)\rho(x,T)\mathrm{d}x$ is well defined, because $\Psi\in C([0,T];L^{2}(\Omega;\mathbb{C}^{N}))$. The directional derivative

[TABLE]

is obviously linear and continuous in $\delta\Psi$, hence it is the real-Gâteaux derivative. Furthermore, we have

[TABLE]

Therefore $J_{\eta}$ is real-Fréchet differentiable from $W$ to $\mathbb{R}$. The Féchet derivative depends linearly on $\Psi$ and is bounded by

[TABLE]

Hence, the derivative is continuous from $W$ to $\mathcal{L}(W,\mathbb{R})$. ∎

4.1 Existence of a minimizer

In this section, we discuss existence of a minimizer of the optimization problem (3.2). Uniqueness cannot be expected because, e.g., the phase of the wave function does not appear in the cost functional $J$.

We start by collecting some known facts. We make use of the following version of the Arzelà-Ascoli theorem which holds also for general Banach spaces, see, e.g., [Cia13, Theorem 3.10-2] and remark thereafter.

Lemma 12** (Arzelà-Ascoli).**

Given a sequence $(f_{n})_{n}$ of functions $f_{n}\in C(K;\mathcal{B})$ where $\mathcal{B}$ is a Banach space and $K=[0,T]$, which is

uniformly bounded, i.e. $\exists M$, such that $\|f_{n}\|_{C(K;\mathcal{B})}\leq M$, $\forall\,n\in\mathbb{N}$, and 2. 2.

equicontinuous, i.e. given any $\epsilon>0$, there exists $\delta(\epsilon)>0$, such that $||f_{n}(t)-f_{n}(s)||_{\mathcal{B}}<\epsilon$ for all $t,s\in K$ with $|t-s|<\delta(\epsilon)$ and all $n\in\mathbb{N}$;

then there exists a subsequence $(f_{n_{l}})_{l}$ and a function $f\in C(K;\mathcal{B})$ such that

[TABLE]

For the purpose of our discussion, notice that the semigroup generated by $H_{0}=-\nabla^{2}+V_{0}$ is denoted by $U(t)=e^{-iH_{0}t}$. For self-adjoint operators $H_{0}$, as in our case, $U(t)$ is strongly continuous by the Stone’s theorem; see [Sto32] and [Yse10, p. 34].

We need the following lemma; see also [RPvL15].

Lemma 13**.**

The Duhamel form of the TDKS equation (2.11) is given by

[TABLE]

Proof.

We follow the approach in [Sal05] and write the following

[TABLE]

which proves the lemma. ∎

Now, we can prove the existence of a solution to (3.2).

Theorem 14**.**

The optimal control problem (3.2) with $\Psi^{0}\in H^{1}_{0}(\Omega;\mathbb{C}^{N})$ and $\beta=0$ admits at least one solution in $(\Psi,u)\in W\times H^{1}(0,T;\mathbb{R})$. In the case $\beta\neq 0$, an optimal solution exists in $W\times H^{1}(0,T;\mathbb{R})$ assuming that $\Psi^{0}\in H^{2}(\Omega;\mathbb{C}^{N})\cap H^{1}_{0}(\Omega;\mathbb{C}^{N})$.

Proof.

For a given control $u\in H^{1}(0,T;\mathbb{R})$, we define $\Psi(u)$ as the unique solution to (2.11). Let $(\Psi_{n},u_{n}):=(\Psi(u_{n}),u_{n})$ be a minimizing sequence of $J$, i.e. $\lim_{n\rightarrow\infty}J(\Psi(u_{n}),u_{n})=\inf_{u\in H^{1}(0,T;\mathbb{R})}J(\Psi(u),u)$. As $J$ is coercive with respect to $u$ in the $H^{1}$ norm, the sequence $(u_{n})_{n}$ is bounded in $H^{1}(0,T;\mathbb{R})$. Hence, we can extract a weakly convergent subsequence again denoted by $(u_{n})_{n}$, $u_{n}\rightharpoonup\hat{u}$ in $H^{1}(0,T;\mathbb{R})$. By the Rellich-Kondrachov theorem, we then have $u_{n}\rightarrow\hat{u}$ in $C[0,T]$.

As the controls $u_{n}$ in the sequence above are globally bounded in $H^{1}(0,T;\mathbb{R})$, by [SCB17, Theorem 3] we have $\|\Psi_{n}\|_{X}\leq K$ and $\|\Psi_{n}^{\prime}\|_{X^{*}}\leq K^{\prime}$, where the constants $K$, $K^{\prime}$ can be chosen to be independent of $u_{n}$. Hence, we can extract weakly convergent subsequences, again denoted by $\Psi_{n}$, $\Psi_{n}^{\prime}$, as follows

[TABLE]

By the Rellich-Kondrachov theorem $H^{1}(\Omega;\mathbb{C}^{N})\Subset L^{2}(\Omega;\mathbb{C}^{N})$ and by [Lio69, 1.5.2] $W\Subset Y$, hence

[TABLE]

By Lemma 9 and 10 and Assumption 2, $\Psi\mapsto V_{Hxc}(\Psi)\Psi$ is real-Fréchet differentiable, hence continuous from $W$ to $Y$. Every $(\Psi_{n},u_{n})$ solves the Schrödinger equation (2.11) and with the strong convergence of $u_{n}$ and $\Psi_{n}$, we can employ [Cia13, p. 291] for the products $u_{n}\Psi_{n}$ and $V_{Hxc}(\Psi_{n})\Psi_{n}$. By standard results, a sequence converging in the $L^{2}(0,T)$-norm contains a subsequence that converges a.e. in $[0,T]$. Hence, we can extract a subsequence, again denoted by $(\Psi_{n})$ such that we can pass to the limit. We have

[TABLE]

holds a.e. in $[0,T]$ and for all $\Phi\in H^{1}_{0}(\Omega;\mathbb{C}^{N})$, where $L:=L^{2}(\Omega;C^{N})$. Hence, $\Psi(\hat{u})=\hat{\Psi}$ and $(\hat{\Psi},\hat{u})$ solves (2.11).

If $\eta\neq 0$, then the evaluation of $\rho$ at the final time $T$ in the cost $J$ is required, such that strong convergence in $C([0,T];L^{2}(\Omega;\mathbb{C}^{N}))$ is needed. So far, we only have weak convergence in $W$, which is continuously embedded into $C([0,T];L^{2}(\Omega;\mathbb{C}^{N}))$, but not compactly embedded. To overcome this problem, we improve our convergence result by using Lemma 12. The required uniform bound is given by Lemma 4, we have $\|\Psi_{n}(t)\|_{L^{2}}=1$, for all $n$ and all $t\in[0,T]$. We are left to show the equicontinuity of $\Psi_{n}$ for the Arzelà-Ascoli theorem, i.e. find an $\delta(\epsilon)$ that does not depend on the $n$ of the sequence. To this end, we take a fixed but arbitrary $\epsilon>0$. With $\|\Psi_{n}\|_{L^{2}}=1$, the Lipschitz continuity of $V_{xc}$, and estimates in [SCB17, Lemma 2, Theorem 3], we have

[TABLE]

where $g$ is defined in Lemma 13. As $U(t)=e^{-iH_{0}t}$ is continuous by Stone’s theorem, there exists a $\delta_{1}$ such that

[TABLE]

Using the Duhamel form from Lemma 13, we find for two different times $t$, $t^{\prime}$, with $|t-t^{\prime}|<\delta$, the following

[TABLE]

The first term is bounded by

[TABLE]

For the second term, we have

[TABLE]

Using the fact that that $e^{iH_{0}s}$ is unitary, we obtain

[TABLE]

where we used the Lipschitz continuity of $V_{xc}=V_{x}+V_{c}$ and [SCB17, Lemma 2] for the estimate on $V_{H}$. Furthermore, we used Lemma 4 $\|\Psi_{n}\|_{L^{2}(\Omega;\mathbb{C}^{N})}=1$, and according to [SCB17, Theorem 3] we have $\|\Psi(t)\|_{H^{1}(\Omega;\mathbb{C}^{N})}\leq C(\Psi^{0})$. Moreover, since the $u_{n}$ are from a bounded sequence, then $\|u_{n}\|_{C[0,T]}$ is bounded by a global constant.

Now, we define $\delta:=\min\{\delta_{1},\frac{\epsilon}{2\sqrt{K}}\}$. Then

[TABLE]

This means, that the sequence $(\Psi_{n})_{n}$ is equicontinous. As it is also uniformly bounded, we can invoke the Arzelà-Ascoli theorem (Lemma 12) to conclude that there exists a subsequence $\Psi_{n_{l}}$ and a function $\hat{\Psi}\in C(0,T;L^{2}(\Omega;\mathbb{C}))$, such that $\lim_{l\rightarrow\infty}\Psi_{n_{l}}=\hat{\Psi}$ in $C(0,T;L^{2}(\Omega;\mathbb{C}))$. With the strong convergence in $C(0,T;L^{2}(\Omega;\mathbb{C}))$, we have

[TABLE]

With this result, we have the convergence of (4.1) also for a target depending only on the final time.

For the target $J_{\beta}$ to be well defined, we need to assume higher regularity. Assuming $\Psi^{0}\in H^{2}(\Omega;\mathbb{C}^{N})$, we get $\Psi\in L^{\infty}(0,T;H^{2}(\Omega;\mathbb{C}^{N}))$ from [SCB17]. Due to the embedding $L^{\infty}(0,T;H^{2}(\Omega;\mathbb{C}^{N}))\hookrightarrow L^{\infty}(0,T;L^{\infty}(\Omega;\mathbb{C}))$ the wavefunction is therefore globally bounded in space and time by a constant depending on the initial condition and $\|u\|_{H^{1}(0,T;\mathbb{R})}$. As $\|u_{n}\|_{H^{1}(0,T;\mathbb{R})}$ is bounded, we have $\|\Psi_{n}\|_{L^{\infty}(0,T;L^{\infty}(\Omega;\mathbb{C}))}<K$, $\|\hat{\Psi}\|_{L^{\infty}(0,T;L^{\infty}(\Omega;\mathbb{C}))}<K$ and the same holds for the squares $\rho_{n},\hat{\rho}$.

As $\Psi_{n}\stackrel{{\scriptstyle Y}}{{\rightarrow}}\hat{\Psi}$, so does $\overline{\Psi}_{n}\stackrel{{\scriptstyle Y}}{{\rightarrow}}\overline{\hat{\Psi}}$. Hence the product converges in the $L^{1}$-norm, $\rho_{n}\stackrel{{\scriptstyle L^{1}(0,T;L^{1}(\Omega;\mathbb{R}))}}{{\longrightarrow}}\hat{\rho}$. Furthermore, $(\rho_{n})_{n}$ converges in the $Y$-norm as follows

[TABLE]

As $J_{\beta}(\Psi)=\frac{\beta}{2}\|\rho-\rho_{d}\|_{Y}^{2}$ and the norm is continuous, we can pass the limit.

The norm $\|u\|_{H^{1}(0,T;\mathbb{R})}$ is weakly lower semicontinuous. Hence, $(\hat{\Psi},\hat{u})$ minimizes (3.2),

[TABLE]

4.2 Necessary optimality conditions

In this section, we discuss the Lagrange multiplier $\Lambda$ and state the first-order optimality condition for a minimum.

We start showing the existence of a Lagrange multiplier in the Lagrange framework.

Theorem 15**.**

Given a control $u$ and a corresponding state $\Psi$, there exists a Lagrange multiplier $\Lambda\in X$ associated with $(\Psi,u)$.

Proof.

The constraint $c$ and the objective $J$ are continuously real-Fréchet differentiable by Lemma 11 and Theorem 6, and the derivative $\operatorname{D}c(\Psi,u):W\times H^{1}(0,T;\mathbb{R})\rightarrow X^{*}$ is surjective by [SCB17], where the results can readily extended to a nonzero right-hand side $F\in X^{*}$. Hence the constraint qualification of Zowe and Kurcyusz is fulfilled [ZK79]. Therefore, we have the existence of a Lagrange multiplier $\Lambda\in X$; see, e.g., [Trö10, Section 6.1]. ∎

As we need higher regularity, namely $\Lambda\in W$, we do not make further use of this result. Instead, we proceed in a different way using the existence of a unique solution of the adjoint equation in $W$ from [SCB17].

As $c(\Psi,u)=0$ is uniquely solvable, we have the following equivalent formulation of the optimization problem.

Lemma 16**.**

The minimization problem (3.2) is equivalent to the unconstrained minimization of the reduced cost functional $\hat{J}(u):=J(\Psi(u),u)$,

[TABLE]

To calculate the gradient of the reduced cost functional, we make use of the implicit function theorem; see, e.g., [Cia13, p. 548]. We apply the chain rule for the real-Fréchet derivative to the cost functional as follows.

Lemma 17**.**

The derivative of the reduced cost functional is given by

[TABLE]

where $\operatorname{D}_{u}\Psi(u)\delta u=\delta\psi$ is the solution to the linearized equation

[TABLE]

Proof.

The chain rule for the real-Fréchet derivative is given by

[TABLE]

To show that the term $\operatorname{D}_{u}\Psi(u)$ is well defined, we use the implicit function theorem, which ensures differentiability of the map $u\mapsto\Psi(u)$, $H^{1}(0,T;\mathbb{R})\rightarrow W$. To this end, we show that the Fréchet derivative $\operatorname{D}_{\Psi}c(\Psi,u):W\rightarrow X^{*}$ is a bijection at any $(\Psi,u)\in W\times H^{1}(0,T;\mathbb{R})$.

We consider the derivative of $c$ with respect to $\delta\Psi$. This is given by

[TABLE]

As shown in [SCB17], the equation $\operatorname{D}_{\Psi}c(\Psi,u)(\delta\Psi)=F$ with the initial condition $\delta\Psi(t=0)=0$ is uniquely solvable for any right-hand side $F\in X^{*}$. Hence $\operatorname{D}_{\Psi}c(\hat{\Psi},\hat{u})$ is a bijection.

Now, from $c(\Psi,u)=0$, we have

[TABLE]

Solving this equation for $\delta\Psi$ results in the following

[TABLE]

This means that the solution $\delta\Psi$ of the linearized equation $\operatorname{D}c=0$ is in fact $\operatorname{D}_{u}\Psi(u)$.

All together, we find

[TABLE]

where $\delta\Psi$ is given by (4.11). ∎

Theorem 18**.**

The $H^{1}$-gradient of the reduced optimization problem is given by

[TABLE]

where $\mu$ is the $H^{1}$-Riesz representative of the continuous linear functional$\left(-\operatorname{Re}\left(\Lambda,\,V_{u}\Psi\right)_{L^{2}(\Omega;\mathbb{C}^{N})},\,\cdot\right)_{L^{2}(0,T;\mathbb{R})}$, $\Psi$ is the unique solution of $c(\Psi,u)=0$, and $\Lambda$ is the unique solution of $a(\Psi,u,\Lambda)=0$, where

[TABLE]

with the terminal condition

[TABLE]

Proof.

We want to calculate the Riesz representative of $\operatorname{D}\hat{J}(u)(\delta u)$. By Lemma 17, we have

[TABLE]

The derivative $\operatorname{D}_{u}J(\Psi(u),u)(\delta u)=\nu\left(u,\,\delta u\right)_{H^{1}(0,T;\mathbb{R})}$, hence the gradient of the second term is given by

[TABLE]

To express the first term as an operator acting on $\delta u$, we need to use the fact from Lemma 17, that $\operatorname{D}_{u}\Psi(u)(\delta u)=\delta\Psi$ is the solution of the linearized equation. Using the directional derivative calulated in the Appendix A.2.1, we obtain for the first term

[TABLE]

To simplify the equation, we focus on the second term in (4.16). By the continuous embedding $W\hookrightarrow C([0,T];L^{2}(\Omega;\mathbb{C}^{N}))$, we can invoke the fundamental theorem of calculus in time.

[TABLE]

Observing that $\delta\Psi(0)=0$ and using equations (4.13), (4.14) for $\Lambda$ and the fact (4.4) that $\delta\Psi$ is the solution of the linearized equation, as well as the previous result, we obtain

[TABLE]

Using integration by parts with the zero boundary condition, this simplifies to

[TABLE]

Using this result in (4.16), we obtain

[TABLE]

as, by definition of $\mu$, $\left(-\operatorname{Re}\left(\Lambda,\,V_{u}\Psi\right)_{L^{2}(\Omega;\mathbb{C}^{N})},\,\delta u\right)_{L^{2}(0,T;\mathbb{R})}=\left(\mu,\,\delta u\right)_{H^{1}(0,T;\mathbb{R})}$.

Adding (4.17) and (4.15) together, we obtain (4.12). ∎

Theorem 19**.**

Given a local solution $(\Psi,u)\in W\times H^{1}(0,T;\mathbb{R})$ of the minimization problem (3.2), i.e. $\hat{J}(u)\leq\hat{J}(\tilde{u})$ for all $\|u-\tilde{u}\|_{H^{1}(0,T;\mathbb{R})}<\epsilon$ for some fixed $\epsilon>0$. Then there exists a unique Lagrange multiplier $\Lambda\in W$, such that the following first order optimality system is fulfilled.

[TABLE]

where $\nabla\hat{J}(u)$ is given by Theorem 18.

Proof.

For all admissible pairs $(\Psi,u)$ that satisfy (4.18a), the adjoint problem (4.18b) has a unique solution in $W$ by [SCB17]. As $\hat{J}$ is differentiable, a local minimum is characterized by a zero gradient; see, e.g., [BS12], hence (4.18c) holds. ∎

5 Numerical approximation and optimization schemes

In this section, we illustrate an approximation scheme for the TDKS equation and its adjoint and discuss an optimization algorithm to solve the optimality system (3.4a)–(3.4e) (or equivalently (4.18a)–(4.18c)).

To solve the TDKS equations (3.4a), we use the Strang splitting [FOS15] given by

[TABLE]

where $V=V_{ext}+V_{Hxc}$ and $j=1,\dotsc,N$. To solve the adjoint TDKS equation (3.4c), we have to include the inhomogeneous right-hand side as follows.

[TABLE]

with a right-hand side $g(t)$, e.g. $g_{j}(t)=-2\beta(\rho(t)-\rho_{d})\psi_{j}(t)$. The Laplacian $-\nabla^{2}$ is evaluated spectrally. With this setting, we obtain a discretization scheme that provides second-order convergence in time and analytic convergence in space. This is proved for constant potentials in [BJM02], and there is numerical evidence that this accuracy performance also holds in the case of variable potentials.

Furthermore, this time-splitting scheme is unconditional stable and norm preserving, as well as time reversible and gauge invariant. The latter means that adding a constant to the potential changes only the phase of the wave function in such a way that discrete quadratic observables are not changed [BJM02].

We are not able to give a proof for the convergence rate of the time-splitting discretization scheme in the presence of nonlinearity as it appears in the TDKS equation. However, we refer to results on the Strang splitting for similar SE. For the inhomogeneous case, the importance of an inhomogeneity that vanishes at the boundary is stressed in [FOS15]. In [Tha12] the Strang splitting for the Gross-Pitaevskii equation is studied, and second-order convergence in time and spectral convergence in space is proved. Further, quantum models with Lipschitz nonlinearities are studied in [BBD02], where second-order convergence in time is obtained.

Even though these results cannot be applied directly to our problem, they suggest that similar accuracy can be expected in our case. Therefore, we study this accuracy issue numerically. To this end, we solve the TDKS equations for two interacting particles in a harmonic trap $V_{ext}=50x^{2}$. The used initial condition is given by the coherent states of two non-interacting particles in the harmonic oscillator. We set $\Omega=(0,L)^{2}$, $L=7$, and the time interval is $[0,0.1]$. Since an analytic solution is not available, we consider a reference solution $\Psi_{\text{reference}}$ obtained solving the problem on a very fine mesh. In Figure 2, we report results of numerical experiments showing second-order convergence in time (slope factor $1.995$) and spectral convergence in space.

Next, we illustrate our implementation of the nonlinear conjugate gradient (NCG) method to solve our optimization problem. We follow the approach in [HZ05]. The minimization algorithm is given in Algorithm 1, where Algorithm 2 is called to compute the reduced gradient.

6 Numerical experiments

In this section, we present results of numerical experiments to numerically validate our optimization framework.

In all experiments, we consider $N=2$ interacting electrons in $n=2$, so $x=(x_{1},x_{2})^{T}\in\mathbb{R}^{2}$. Confining electrons to a two-dimensional surface region is used to model quantum dots, which are nowadays widely used in applications; see, e.g., [HKP*+*07] for a review on quantum dots.

As initial guess for the control, we take $u_{0}=0$. Different results with different choices of the optimization parameters are denoted differently. In the figures below, the results obtained with $\nu=10^{-5}$ are shown with dotted lines; in the case $\nu=10^{-6}$ we plot dash-dotted lines, and with $\nu=10^{-7}$, we use dashed lines.

In our first experiment, our objective is that the density of 2 electrons follows a prescribed trajectory ($\beta=1$, $\eta=0$). Our target trajectory is produced by an oscillating strength of the harmonic confinement $V_{ext}(x,t,u)=50(x_{1}^{2}+x_{2}^{2})+u_{d}(t)(x_{1}^{2}+x_{2}^{2})$ with a forcing $u_{d}$. Therefore, our purpose is to track the density resulting from the prescribed forcing term $u_{d}$ (solid line in Figure 3(a)). The stopping criterion for convergence is $\|\nabla J\|_{H^{1}(0,T;\mathbb{R})}<5\cdot 10^{-7}$.

The results of this experiment are presented in Figure 3. We see that the trajectory is tracked more closely for smaller values of $\nu$, which require larger computational effort, while with larger $\nu$ the stopping criterion is met after fewer iterations.

The second experiment is as in [CWG12]. In this case, we consider two electrons in the following asymmetric double well potential

[TABLE]

At $t=0$, the electrons are in their ground state which is centred around the global minimum at $(x_{1},x_{2})=(-3.6,0)$. Our objective is to spatially shift this ground state, that is, to move the 2 electrons to the right-half space, $x_{1}>0$. For this purpose, we consider the following cost functional

[TABLE]

The stopping criterion for convergence is $\|\nabla J\|_{H^{1}(0,T;\mathbb{R})}<5\cdot 10^{-5}$.

As the results presented in Figure 4 show, the cost functional can be reduced by approximately 4 orders of magnitude and the density is almost completely localized in the desired set, since $\int_{x_{1}<0}\rho(x,T)\mathrm{d}x=2.6\cdot 10^{-4}$. This result is obtained with $\nu=10^{-7}$. For smaller values of $\nu$, the stopping criterion is matched later and the objective can be further improved.

7 Conclusion

In this paper, an optimal control framework for the time-dependent Kohn-Sham model was presented and analyzed. The purpose of this work was to provide a mathematical rigorous proof of the existence of optimal controls and their characterization as solutions to the corresponding optimality system. For this purpose, the differentiability properties of the nonlinear Kohn-Sham potential and of the TDKS equation were investigated. A proof of the existence of a minimizer was given and the first-order optimality system was discussed. To validate the proposed optimization framework, an efficient Strang splitting discretization scheme was implemented and validated. The optimization problem was solved using a NCG scheme. Results of numerical experiments demonstrated the control ability of our scheme with different control problems.

Appendix A Derivation of the optimality system

The derivation of the optimality system is done by calculating the directional derivatives of the Lagrange functional (3.3) with respect to the states $\psi_{j}$, the adjoint variables $\lambda_{j}$, and the control $u$.

A.1 The forward equation (derivative by $\lambda_{j}$)

Since $J$ does not depend on $\Lambda$ and $L_{1}$ depends linearly on $\lambda_{j}$, so the derivative by $\lambda_{j}$ gives the TDKS equation for $\psi_{j}$:

[TABLE]

A.2 The adjoint equation (derivative by $\psi_{j}$)

Both $L_{1}$ and the target term $J$ depend on $\Psi$. We start with $L_{1}$ which is split into two terms: the linear part $L_{2}$ and the nonlinear part from the Kohn-Sham potential $L_{3}$. We have

[TABLE]

We begin with the linear part $L_{2}$. Differentiating $L_{2}$ by $\psi_{j}$ gives the TDKS equation for $\lambda_{j}$ as it is linear in $\psi_{j}$ and the sign from integration by parts is canceled by the complex conjugation. One boundary term $B_{1}$ appears.

[TABLE]

For the calculation details, we fix one particle index $j$:

[TABLE]

Now, we use integration by parts and use the fact that $\lambda_{j}$ and $\psi_{j}$ are zero on the boundary. We obtain

[TABLE]

At time $t=0$, all wave functions have to fulfill the initial condition, hence $\delta\psi_{j}(x,0)=0$. We are left with one boundary term

[TABLE]

which can be removed by prescribing the terminal condition $\lambda_{j}(x,T)=0$ in the case of $\eta=0$, otherwise by the terminal condition derived below. We obtain

[TABLE]

The derivative of the nonlinear part $L_{3}$ is given in the Lemmas 7 and 10.

A.2.1 Derivative of the target functional

Consider the trajectory term $J_{\beta}=\frac{\beta}{2}\int_{0}^{T}\int_{\Omega}(\rho(x,t)-\rho_{d}(x,t))^{2}\mathrm{d}x\mathrm{d}t$, and notice the following

[TABLE]

With this preparation, we habe

[TABLE]

Similarly, we find for the terminal term $J_{\eta}=\int_{\Omega}\rho(x,T)\chi_{A}(x)\mathrm{d}x$ the following

[TABLE]