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

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

M. Sprengel Institut für Mathematik, Universität Würzburg, Emil-Fischer-Strasse 30, 97074 Würzburg, Germany (martin.sprengel@mathematik.uni-wuerzburg.de).    G. Ciaramella Section de mathématiques, Université de Genève, 2-4 rue du Lièvre 1211 Genève 4, Switzerland (gabriele.ciaramella@unige.ch).    A. Borzì Institut für Mathematik, Universität Würzburg, Emil-Fischer-Strasse 30, 97074 Würzburg, Germany (alfio.borzi@mathematik.uni-wuerzburg.de).
July 3, 2019
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ö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.

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 interacting particles in an external potential by another system of non-interacting particles with an external potential , 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 [MUN06] 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 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 [CAO06]. 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 electrons system is described by a wave function , whose time evolution is governed by the following Schrödinger equation (SE)

(2.1)

where are the position vectors of the particles. We use atomic units, i.e. , .

The Hamiltonian consists of a kinetic term, the Coulomb interaction between the charged particles (electrons), , and an external potential, . We have

(2.2)

where is the Euclidean norm and is the -dimensional vector gradient with respect to .

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 orthogonal single particle wave functions (orbitals), , that correspond to the Hamiltonian of the th particle, one can build the following antisymmetric wave function

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 -dimensional space, the corresponding electronic density defined of the physical space of -dimensions given by

(2.3)

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

(2.4)

where represents the wave function of the th particle.

The DFT approach of Kohn and Sham [KS65] to model multi-particle problems was to replace the system of interacting particles, subject to an external potential , by another system of non-interacting particles subject to an augmented potential , 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

(2.5)
(2.6)

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 eigenvalues as .

Notice that an TDKS system is formulated in spatial dimensions and consists of coupled Schrödinger equations. With the appropriate choice of , which contains the coupling through the dependence on , 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 only depends on the value of the density at this specific point. We use the adiabatic LDA, which means that LDA is applied at every time separately such that . Notice that, if one allows to depend on the whole history , , 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 is usually split into three terms, that is, the Hartree potential, and the exchange and correlation potentials as follows

(2.7)

The Hartree potential models electrons interaction due to the Coulomb force. This term dominates , while the other two parts represent quantum mechanical corrections. Later, we refer to the exchange-correlation part of the potential also as . Notice that the classical electric field is given by , where and are positions in space and is the charge density. Since we are using natural units, charge and particle densities become the same and we have the following

(2.8)

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 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]

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

(2.9)

where

(2.10)

with and and sufficiently large; e.g., . This potential is twice continuously differentiable and globally bounded.

The remaining part of the interaction is called the correlation potential . 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 used in the numerical experiments in Section 6.

Figure 1: The Quantum Monte Carlo fit used in the numerical experiments of this work from [AMGGB02]. The values for the limits are and .

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 in the sense that for all .

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 is non-zero only on a bounded domain . For this reason, we denote by 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, .

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 , where has the role of a modulating amplitude. A specific case is the dipole control potential, .

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

In particular, we consider the control of a quantum dot by a changing gate voltage modeled by a variable quadratic potential, , and a laser control in dipole approximation, , with a polarization vector .

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 , , with and , together with initial conditions , . For brevity, we denote our KS wavefunction by , and . Therefore, we have .

We define the following function spaces. and with the norms and , where is the dual space of ; we also need 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 instead of . We shall also omit the explicit dependence of on if no confusion may arise.

Lemma 1.

The exchange potential term is Lipschitz continuous, i.e. and .

Proof.

The function , is continuously differentiable with bounded derivative, hence Lipschitz continuous with Lipschitz constant from to . With this preparation, we have Lipschitz continuity from to as follows

Similarly, we have Lipschitz continuity from to as follows

where we use the Gelfand triple and the fact that as . ∎

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

Assumption.
  1. The correlation potential is uniformly bounded in the sense that , , ; this corresponds to the correlation potentials from the Libxc library [MOB12] applying a similar approach to (2.9).

  2. The correlation potential term 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. The confining potential and the spacial dependence of the control potential are bounded, i.e. ; as we consider a finite domain, this is equivalent to excluding divergent external potentials.

  4. The control is . 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 , , admits a unique solution in , that is, there exists , , such that

(2.11)

for all and a.e. in .

If and , then the unique solution to (2.11) is as follows

(2.12)

if in addition , we have .

By the continuous embedding , see e.g. [Eva10, p. 287], the solution is continuous in time.

Similar problems have been studied in [Jer15]:

Theorem 3.

Assuming that and , and a Lipschitz condition on and a continuity assumption on , then (2.11) with has a unique solution in .

Taking into account only the Hartree potential but not the exchange-correlation potential, existence of a unique solution in 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 -norm of the solution to (2.11) is conserved in the sense that for all .

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 that includes different target functionals and a control cost as follows

(3.1)

The first term, models the requirement that the electron density evolves following as close as possible a given target trajectory, . We remark that is only well-defined if is at least in . This is guaranteed by the improved regularity from Theorem 2 for . The term aims at locating the density outside of a certain region . The term penalizes the cost of the control. We remark that the regularization term can be any weighed norm, e.g. for . We assume that the target weights are all non-negative , with , and the regularization weight . The characteristic function of is given by .

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

(3.2)

Here and in the following, we assume that .

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

(3.3)

The function , where , , represent the adjoint variables. In the Lagrange formalism, a solution to (3.2) corresponds to a stationary point of , where the derivatives of with respect to , and must be zero along any directions , , and . 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 , in (3.4c) below.

The first-order optimality conditions define the following optimality system \cref@addtoresetequationparentequation

(3.4a)
(3.4b)
(3.4c)
(3.4d)
(3.4e)

where , and

Further, is the -Riesz representative of the continuous linear functional; see Theorem 18 below. Assuming that , can be computed by solving the equation

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

(3.5)

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 from these solutions in Theorem 18 below. Later, we use this gradient to construct a numerical optimization scheme to minimize .

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 are continuously real-Fréchet differentiable. Subsequently, we prove existence of solutions to the optimization problem.

The Kohn-Sham potential depends on the density . 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 be complex Banach spaces. A map is called real-linear if and only if

  1. and

  2. and .

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

We call a map 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 also hold in 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

and the corresponding matrix Schrödinger equation.

Theorem 6.

The map , defined as

(4.1)

is continuously real-Fréchet differentiable.

Notice that represents the TDKS equation.

Remark.

We remark that , , and are in . Hence, the operator norm of their derivatives in can be bounded in as follows.

where denotes the derivative of , , or , respectively.

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

Lemma 7.

The map is real-Gâteaux differentiable for all and its real-Gâteaux derivative, denoted with , is specified as follows

Proof.

First, we need the derivative of the density which is a non-holomorphic function of a complex variable. Differentiation of can be done using the Wirtinger calculus [Rem91] where and are treated as independent variables. By using these calculus rules, the real-Fréchet derivatives of are given by

The directional derivative of along is given by ,

Using the definition of , we have

and is obviously linear in over the real scalars. We are left to show that is a bounded operator.

For the first term, we use as well as to obtain

For the second term, we decompose the domain into depending on the size of : , , and . Using the fact that is bounded by in and by in as well as that in , and that is monotonically decreasing between and we obtain

We have shown that the directional derivative is a bounded linear map for all , hence is real-Gâteaux differentiable. ∎

Before improving this result to real-Fréchet differentiability, we prove a general result on the -norm of products of functions. To this end, denote .

Lemma 8.

Let be arbitrary. Given two functions and . Then . Similarly, for and , we have.

Proof.

The Fourier transform is an isometry by Plancherel’s theorem, where we extend the functions by zero outside of . Furthermore, for , the convolution theorem states . Hence, using the embeddings , we find

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 . ∎

Lemma 9.

The map is continuously real-Fréchet differentiable with derivative , where is given in Lemma 7.

Proof.

We prove that the real-Gâteaux derivative of at is continuous from to . Then the real-Fréchet differentiability follows immediately from [AH11, Proposition A.3].

Once is proved that 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

For this purpose, to ease our discussion, we consider , , where

such that . For all , is continuously differentiable with bounded derivative for , hence the same holds for , which is therefore Lipschitz continuous. In zero, is Hölder continuous with exponent ,

Together, we have Hölder continuity for all