Direct Optimal Control Approach to Laser-Driven Quantum Particle Dynamics

Optimal control theory is usually formulated as an indirect method requiring the solution of a two-point boundary value problem. Practically, the solution is obtained by iterative forward and backward propagation of quantum wavepackets. Here, we propose direct optimal control as a robust and flexible alternative. It is based on a discretization of the dynamical equations resulting in a nonlinear optimization problem. The method is illustrated for the case of laser-driven wavepacket dynamics in a bistable potential. The wavepacket is parameterized in terms of a single Gaussian function and field optimization is performed for a wide range of particle masses and lengths of the control interval. Using the optimized field in a full quantum propagation still yields reasonable control yields for most of the considered cases. Analysis of the deviations leads to conditions which have to be fulfilled to make the semiclassical single Gaussian approximation meaningful for field optimization.

In quantum optimal control theory the goal of optimizing the expectation value of a target operator such as a projector onto a certain state, is formulated as a variational problem for a cost functional subject to certain constraints. The latter includes, for instance, some penalty for high field intensities or that the wavepacket should fulfill the Schrödinger equation. This control problem is usually solved using an indirect approach, i.e., the cost functional is not minimized directly. Instead, the stationarity condition for the cost functional is converted to a two-point boundary problem for two coupled Schrödinger equations. A numerical solution is obtained by iterative forward and backward propagation of the actual wavepacket and an auxiliary wavepacket, respectively (e.g., [18]). This procedure is sometimes referred to as the optimize and then discretize paradigm [19]. Indirect methods for optimal control are in use in other areas of physics, e.g., stochastic control [20], but also in engineering and biology [21].
Direct optimal control, in contrast, follows the discretize and then optimize paradigm, i.e. the cost functional is minimized directly using methods from nonlinear optimization. Although being popular, for instance, in applied mathematics [22], engineering [23], and biology [21], there have been no applications to quantum molecular dynamics so far. The present paper is devoted to fill this gap.
Indirect optimal control requires to solve iteratively two timedependent Schrödinger equations where the numerical effort scales exponentially with the number of degrees of freedom. To cope with this situation the Multi-Configurational Time-Dependent Hartree (MCTDH) approach is most suited [24,25]. An OCT implementation has been reported in Ref. [26], for an application see also Ref. [27]. The solution of the timedependent Schrödinger equation requires a priori knowledge of the potential energy surface. But, when driving the wavepacket into a particular region of configuration space using laser control, a global potential might not be needed. Thus on-the-fly approaches, e.g., in the context of MCTDH [28,29] could be of advantage. On the other hand, semiclassical approximations in terms of Gaussian wavepackets play a prominent role in molecular quantum dynamics [30] and indeed there has been a semiclassical formulation of indirect OCT reported in Refs. [31,32] (for related work using Wigner space sampling, see Ref. [33]).
In this paper we explore direct OCT using a representation of the wavepacket dynamics in terms of a single Gaussian function. Although this choice has been made for numerical convenience, it also facilitates exploration of its limitations by comparison with solutions of the time-dependent Schrödinger equation. Specifically, for the considered problem of quantum particle motion in a bistable potential we are able to identify conditions for which the single Gaussian approximation is adequate.

Equations of Motion
The equations for the time evolution of a quantum mechanical state can be obtained from the time-dependent variational principle starting with the stationarity condition for the action S, i.e. [34].
where the quantum Lagrangian is given by (Note that atomic units are used throughout) In the following we will focus on one-dimensional systems (coordinate x and momentum p) coupled to a radiation field, E(t), in dipole approximation (dipole operator μ(x)). Thus the Hamiltonian operator in the coordinate representation is given by Equation 1 yields the condition [34].
Assuming that the time-dependence of the wavepacket is implicitly parameterized by the set of time-dependent real parameters a(t) {a 1 (t), . . . , a n (t)}, this yields δΨ n j 1 zΨ za j δa j .
Inserting Eq. 5 into Eq. 4 gives the equations of motion for the general set of parameters used to describe the wavepacket with K ij being the elements of the inverse of the matrix formed by Im〈zΨ/za i zΨ/za j 〉.
In order to connect to on-the-fly approaches and to reduce the number of differential equations of motion (and thus the computational cost) we assume that the wavepacket has the following Gaussian form [30] at all times where α and β describe the width and tilt of the phase space Gaussian. Further, x 0 and p 0 are the average position and momentum, respectively. Hence, we identify a(t) {α(t), β(t), x 0 (t), p 0 (t)} and using Eq. 6 gives the following set of coupled differential equations subject to some initial conditions at time t 0 . Here, we defined the time-dependent expectation value of the potential In the next section we will focus on the control problem assuming that these equations of motion can be solved, which implies that the expectation value of the potential and its derivatives are available.

Statement of the Control Problem
Let us start with a brief summary of optimal control theory [9,10,35]. Given a functional of the form (13) where T and R are the terminal and running cost, respectively, the task is to find the state trajectory a(t), external control u(t) (where the time t ∈ [t 0 , t f ]) and the set of static parameters k that minimize the functional J [a, u, k]. The minimization is performed subject to the following differential constraints Further, there can be path constraints and event constraints such as Here, the subscript L and and U denotes the lower and upper boundary, respectively, defining the constraints. Notice that in contrast to path constraints, event constraints are not timedependent, but could include a functional, F, of, e.g., the state trajectory or the external control (see below).
Next, we specify this general control problem to the model introduced in Section 2.1. The state is characterized by the set a(t) {α(t), β(t), x 0 (t), p 0 (t)} and the external control is given by the laser field u(t) E(t). Additional time-independent parameters, k, will not be used. The differential constraints (14) are given by Eqs. [8][9][10][11].
The goal of the optimization can be stated as follows. Given some initial quantum state |Ψ(t 0 )〉, parameterized by such that the overlap is maximized between the time-evolved final state at t t f , Ψ(t f )〉, and some target state Φ t 〉. Thus, the terminal cost in Eq. 13 is given by (notice the minus sign because the terminal cost will be minimized and we want to maximize the overlap) Here, for simplicity we will use the parametrization of Eq. 7 for the target state as well, labeling the target parameters as a t {α t , β t , x t 0 , p t 0 }. The running cost will be chosen as follows Besides the field intensity we have included a factor κ scaling the penalty for high field strengths as well as a shape function s(t), which ensures that the field increases(decreases) slowly when turned on(off) [36]. Note that ϵ is a small parameter introduced to avoid division by zero and numerical problems at times t 0 and t t f . Throughout the text we have used ϵ 0.005.
For the application presented below we don't use any path constraints, but event constraints. Given the event upper and lower bounds will be chosen equal as follows Hence, the parameters of the initial state are fixed and not subject to optimization. Further, we enforce the zero-net-force condition by demanding that The optimization problem will be solved using a direct method, i.e. by means of discretization of the differential equations. Details will be specified in the next section.

Model System and Computational Details
The direct optimal control approach will be applied to the problem of particle dynamics in a bistable potential. This could represent, for instance, proton or hydrogen atom transfer in a tautomerization reaction [38,39]. The following potential will be used Here, x B is the distance between the minimum of the potential and the top of the barrier, and V B is the barrier height.
The system-field interaction is treated in semiclassical approximation, taking the polarization of the field in the same direction as the dipole, and assuming a linear model for the latter (q is the charge) Specific parameters for the numerical simulations have been chosen to mimic typical situations in proton transfer reactions [38,39], i.e. x B 2a 0 ( ≈ 1.06Å), V B 0.01E h ( ≈ 6.3 kcal/mol), and q 1 ( 1e). The particle's mass, m, will be used to tune the 'quantumness' of the dynamics. Exemplary, we show potential and eigenstates for two choices of the masses in Figure 1.
Comparing the two cases we note that in particular the number of eigenstates below the barrier is 8 and 16 for masses of 1 m H and 5 m H respectively (where m H is the hydrogen mass).
Using Eqs. 21, 22 together with Eq. 7 one can calculate the time-dependent expectation value of the potential, Eq. 12, and its derivatives with respect to α and x 0 required for the equations of motion (Eqs. 9, 11). Although in the present case the required expectation value could have been calculated analytically, we have used a more general prescription. To this end the potential is globally approximated by a sum of Gaussians of the form Using Eq. 23 one obtains for Eq. 12 where and For the solution of the control problem the software package PSOPT has been used [40]. This package employs an approximation for the state trajectory of the form where t k are the Gauss-Lobatto quadrature nodes (a(t k ) a N (t k )) and L k are the Lagrange basis polynomials. This approximation allows to transform the performance functional (Eq. 13) into the performance function and the differential constraints into a set of holonomic constraints for the decision vector y (u(t 0 ), . . . , u(t N ), a(t 0 ), . . . , a(t N ), k, t 0 , t f ); w k are the Gauss-Lobatto weights. For more details see Ref. [40]. The performance function (30) is optimized using nonlinear programming (NLP) algorithms, such as the ones implemented in IPOPT [41]. PSOPT provides different discretization schemes. The global pseudospectral Legendre and Chebyshev discretization yield very slow convergence for non-smooth functions [19], as it is the case for the solutions found for α(t) and β(t) (see first and second row, (b) and (d) columns of Figure 2 below). Increasing the number of nodes is not an option for these discretization schemes because of the nonsparsity of the Jacobian matrices which cannot be handled properly by the implemented IPOPT NLP solver. This issue translates into a disproportional increase of computational time. The local methods available are trapezoidal and Hermite-Simpson discretization. In order to check their performance we simulated the case of a particle of mass of 1 m H and a final time of t f 20, 000 au. In doing so the number of time discretization nodes has been scanned from 200 to 6,000. To evaluate the discretization error we use the maximum relative local error, ε disc , defined in Ref. [40]. The results are shown in Figure 3. If the number of nodes is below 1,000 the trapezoidal method has a smaller error ε disc compared to Hermite-Simpson for the same number of nodes. Beyond 1,000 nodes, Hermite-Simpson outperforms the trapezoidal discretization. However, this comes at the expense of an increased computational time as can be seen in the lower panel of Figure 3. For the simulations reported below we have used Hermite-Simpson discretization with 2,000 nodes, which offers a good balance between accuracy and speed.
In order to quantify the importance of quantum effects beyond the simple Gaussian ansatz for the wavepacket, Eq. 7, MCTDH simulations have been performed using the optimized field. For this purpose the Heidelberg MCTDH package has been used [42].

Laser-Controlled Proton Transfer
In the following we present a proof-of-principle application of direct OCT using the example of proton transfer in a bistable potential. Specifically, the two cases (particle masses) given in Figure 1 will be considered. For the initial state we choose the parameters of a Gaussian in the left well, and as the target state we choose a symmetrically located Gaussian in the right side well. The Gaussian parameters have been optimized to the ground state using a local harmonic approximation. Although direct The optimal solutions for the two particle masses are given in Figures 2B,D. Apparently, the optimal field is able to drive the center of the wavepacket across the barrier into the right minimum at t t f . In this respect one should note that the optimal fields have a relatively simple shape and little resemblance with the initial guess. This is one of the major advantages of the direct approach to optimal control problems, i.e. the convergence region of the initial guess is very broad. The dynamics is rather similar, i.e., in both cases the trajectory passes the barrier coming from the turning point at the left hand side. Just before and after the barrier the wavepacket gets localized in coordinate and delocalized in momentum space, whereas the position-momentum correlation (β) vanishes. The wavepacket passes the top of the barrier with large momentum.
The question now arises if the optimum field found for a single Gaussian wavepacket is able to trigger the same particle dynamics in the full quantum case. To this end the optimal field is used within a quantum dynamics simulation. The results are compared in Figure 4 in terms of coordinate and momentum expectation values and the respective standard deviation. Until after the barrier crossing, Gaussian and full quantum results are rather similar. Indeed, if the goal would have been to trigger the localization of the wavepacket somewhere in the region of the right well at a particular time, the optimal field would still perform this task also in the quantum case. Of course, the agreement between classical and quantum propagation is better in case of the heavier mass even though there is considerable larger spread of the wavepacket in the quantum case after reflection at the right turning point. For the lighter mass the agreement after barrier crossing is less favorable due to the larger spread and the structured character of the quantum wavepacket which cannot be captured by a single Gaussian. Frontiers in Physics | www.frontiersin.org April 2021 | Volume 9 | Article 615168

Region of Validity of the Gaussian Wavepacket Approximation
Single Gaussians cannot capture the dynamics of structured wavepackets. Nevertheless, the agreement between Gaussian and full quantum results is at least qualitative, even for the lighter particle. This provides the motivation for the investigation of the validity of the Gaussian approximation over a wider range of parameters. Again the optimum field is obtained following the procedure described in Section 3.1, but now for different final times (ranging from 5, 000 au to 20, 000 au in steps of 1, 000 au) and masses (ranging from 1 to 10 m H in steps of 1 m H ). To evaluate the performance of the optimum field to drive the wavepacket to the right well in the full quantum case we choose the following error: whereΨ(t f ) is the exact quantum wavefunction at the final time. This error will be between 0 and 1 if the expectation value of the quantum wavepacket crossed the barrier and greater than 1 if it did not. Results are shown in Figure 5.   In general, we can see from Figure 5 that the Gaussian optimal control fields are able to drive the particle reaction on a broad range of masses and final times. As expected the performance deteriorates for the lighter masses. There are some features which deserve closer attention. For example, there are regions where the Gaussian wavepacket approach works exceptionally well (characterized by stripes of intense blue color). In these regions the final time is matching a total integer number of well oscillations plus the barrier crossing time. Assuming that these oscillations are harmonic with period T and taking the barrier crossing time as being half of the harmonic period, these final times can be estimated. The middle green line in Figure 5 corresponds to a final time of 5T/2. It nicely matches with the dark blue region where the approach works well. Thus, in general one would expect regions with (2n + 1)T/2 and nT where the approximation works well and not so well, respectively. This is roughly seen in Figure 5, although the deviation from the harmonic approximation causes some quantitative disagreement. This analysis points to the importance of the final time t f for the effect of the quantumness of the dynamics on the overlap with the target. In passing we note that in principle direct optimal control offers the possibility to optimize the final time as well, e.g., to fulfill some constraints with respect to the spread of the wavepacket.
Another interesting feature apparent from Figure 5 are the isolated "islands" of poor performance, e.g. at t f 14, 000 au and m 7 m H . To rationalize this behavior Figure 6 shows various expectation values for t f 14, 000 au and m 6 and 7 m H . The first and second row compares Gaussian and quantum results and we can notice that the corresponding trajectories diverge considerably more for 7 m H (b) than for 6 m H (a), even though a naive consideration would suggest that the performance of the single Gaussian approximation is better for the more massive particle. In general we observe that while in the good performing cases the wavepacket essentially stays localized, the opposite is true for the poor performing cases, which stands out as a likely reason for the discrepancy between Gaussian and quantum propagation in the later case. This holds irrespective of the actual mass of the particle. From the second and fourth rows of Figure 6 we notice that the cases m 6 and 7 m H differ in the momentum and thus kinetic energy when crossing the barrier. While in the former case the momentum is maximum at the barrier top, in the latter the particle is slowed down when reaching the barrier. As a consequence it becomes rather delocalized in position space and thus the single Gaussian approximation fails.
In principle one could expect that decreasing the penalty factor κ would alleviate this problem, i.e. stronger fields would imply higher momentum. However, after inspecting Figure 6, it is apparent that for a given final time it depends on the initial direction of momentum whether the wavepacket will pass the barrier with high or low momentum. This idea supports the conclusion that not only the mass of the particle, but also the specific optimal path, are important for the validity of the single Gaussian approximation. Controlling the initial direction in a way which works in a black-box fashion for all cases covered in Figure 5 has not been successfull. However, in contrast to indirect  control, where one would have to compute running cost derivatives with respect to state variables to get coupling terms between forward and backward Schrödinger equation, including additional running costs is straightforward in direct control. To demonstrate this we have added a second term to the running costs of Eq. 18, which serves to maximize the kinetic energy, i.e.
Here, η is a penalty scaling factor and the minus sign ensures that this term gets maximized. It is expected that this will lead to barrier crossing with high momentum and thus a reduced error, Eq. 31.
The results shown in Figure 7 clearly support our hypothesis, i.e. adding the running cost functional Eq. 32 leads to the elimination of the poor-performing islands. Hence, using the flexibility of the direct optimal control approach the region of validity of the single Gaussian approximation could be extended.

CONCLUSION
In this paper we have introduced a new tool for quantum optimal control. In contrast to indirect methods, which require the solution of a two-point boundary value problem, the present direct method builds on the first discretize and then optimize paradigm. Thus, by construction there is no need for explicit propagation of a wavepacket. So far direct methods have found application mostly in engineering [23,40]. The performance and capabilities of the direct method have been demonstrated for the case of one-dimensional particle transfer in a bistable potential. For simplicity the wavepacket has been approximated by a single Gaussian function, but in principle other forms are possible, e.g., superposition of Gaussians [28] or even expansions in terms of an eigenstate basis. Of course, Gaussians have the potential advantage of being suited for on-the-fly simulations, which brings OCT into the realm of the dynamics of complex molecular systems, at least in principle. At this point it will be required to explore the scaling of the numerical effort associated with the direct method more thouroughly. Here, we merely explored the dependence on the number of nodes. But the number of parameters will be another limiting factor. Preliminary calculations performed on regular hardware showed that about 50 parameters and 500 nodes are feasible.
For a simple test system the question has been addressed whether the quantumness of the dynamics influences the final control yield, given a field which has been optimized for the single Gaussian approximation. Interestingly, it turned out that nearly complete particle transfer can be achieved for a wide range of masses and final times. Here, the important point is whether the wavepacket crosses the barrier with high or low momentum, which for the given model is decided by the sign of the momentum during the initial dynamics. As a consequence, even the optimization based on a simple Gaussian wavepacket, possibly using on-the-fly dynamics, may provide reasonable control fields.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.