ORIGINAL RESEARCH article
Direct Optimal Control Approach to Laser-Driven Quantum Particle Dynamics
- Institut für Physik, Universität Rostock, Rostock, Germany
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.
“Teaching lasers to control molecules” has been a long-standing goal in molecular physics . Among the various methods of the early days [1–5], optical control theory (OCT) emerged as a versatile tool. Originally developed by Rabitz et al. [6, 7] and Kosloff et al. , numerous methodological extensions have been developed over the years (for reviews, see e.g., [9–12]). In terms of practical realizations of chemical reaction control, the feedback strategy [1, 13, 14] as well as straightforward resonant excitation schemes [15–17] have been most successful.
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., ). This procedure is sometimes referred to as the optimize and then discretize paradigm . Indirect methods for optimal control are in use in other areas of physics, e.g., stochastic control , but also in engineering and biology .
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 , engineering , and biology , 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 time-dependent 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. , for an application see also Ref. . The solution of the time-dependent 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  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. ).
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.
2 Theoretical Methods
2.1 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. .
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,
Assuming that the time-dependence of the wavepacket is implicitly parameterized by the set of time-dependent real parameters
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  at all times
where α and β describe the width and tilt of the phase space Gaussian. Further,
subject to some initial conditions at time
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.
2.2 Statement of the Control Problem
Further, there can be path constraints
and event constraints such as
Here, the subscript
Next, we specify this general control problem to the model introduced in Section 2.1. The state is characterized by the set
The goal of the optimization can be stated as follows. Given some initial quantum state
Here, for simplicity we will use the parametrization of Eq. 7 for the target state as well, labeling the target parameters as
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
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.
2.3 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
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.
FIGURE 1. Eigenstates for a particle of mass (A)
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
We have used
For the solution of the control problem the software package PSOPT has been used . This package employs an approximation for the state trajectory of the form
and the differential constraints into a set of holonomic constraints for the decision vector
FIGURE 2. Initial guess (A) and (C) and optimal solution (B) and (D) for state,
FIGURE 3. Maximum relative local error (upper panel) and timing (lower panel) for trapezoidal (blue) and Hermite-Simpson (orange) as a function of the number of nodes.
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 .
3.1 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 control in principle allows to vary the final time, in the present application the final time has been fixed to
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
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.
FIGURE 4. Comparison of the coordinate (top row) and momentum (bottom row) expectation values and their respective standard deviation (shaded area), using the Gaussian approximation (blue) and the full quantum propagation (orange). Both trajectories are propagated under the influence of the optimal control field as obtained for the Gaussian (A)
3.2 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
FIGURE 5. Error according to Eq. 31 as a function of different final times and masses. Green lines represent an odd number of half harmonic oscillation periods for the corresponding mass
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
Another interesting feature apparent from Figure 5 are the isolated “islands” of poor performance, e.g. at
FIGURE 6. Expectation values of coordinate and momentum (shaded areas indicate the standard deviation), optimal field, as well as total (
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.
FIGURE 7. Error according to Eq. 31 as a function of different final times and masses. Running cost according to Eq. 32 has been used together with Eq. 18. The penalty scaling factor was
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  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.
ARR has performed the work and analyzed the results. OK has designed the project and supervised the scientific work. All authors have discussed and interpreted the results and contributed to writing the article.
This work has been funded by the grant Ku952/10–1 from the Deutsche Forschungsgemeinschaft (DFG).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
15. Stensitzki T, Yang Y, Kozich V, Ahmed AA, Kössl F, Kühn O, et al. Acceleration of a ground-state reaction by selective femtosecond-infrared-laser-pulse excitation. Nat Chem 10 (2018) 126–31. doi:10.1038/nchem.2909
16. Nunes CM, Pereira NAM, Reva I, Amado PSM, Cristiano MLS, Fausto R. Bond-breaking/bond-forming reactions by vibrational excitation: infrared-induced bidirectional tautomerization of matrix-isolated thiotropolone. J Phys Chem Lett (2020) 11:8034–9. doi:10.1021/acs.jpclett.0c02272
18. Zhu W, Rabitz H. A rapid monotonically convergent iteration algorithm for quantum optimal control over the expectation value of a positive definite operator. J Chem Phys (1998) 109:385–91. doi:10.1063/1.476575
21. Chen-Charpentier BM, Jackson M. Direct and indirect optimal control applied to plant virus propagation with seasonality and delays. J Comput Appl Math (2020) 380:112983. doi:10.1016/j.cam.2020.112983
23. Pardo D, Moller L, Neunert M, Winkler AW, Buchli J. Evaluating direct transcription and nonlinear optimization methods for robot motion planning. IEEE Robot Autom Lett (2016) 1:946–53. doi:10.1109/lra.2016.2527062
25. Beck M, Jäckle A, Worth GA, Meyer HD. The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets. Phys Rep (2000) 324:1–105. doi:10.1016/s0370-1573(99)00047-2
28. Richings GW, Polyak I, Spinlove KE, Worth GA, Burghardt I, Lasorne B. Quantum dynamics simulations using Gaussian wavepackets: the vMCG method. Int Rev Phys Chem (2015) 34:269–308. doi:10.1080/0144235x.2015.1051354
36. Sundermann K, de Vivie-Riedle R. Extensions to quantum optimal control algorithms and applications to special problems in state selective molecular dynamics. J Chem Phys (2000) 110:1896. doi:10.1063/1.477856
38. Došlić N, Kühn O, Manz J. Infrared laser pulse controlled ultrafast H-atom switching in two-dimensional asymmetric double well potentials. Ber Bunsen Ges Phys Chem (1998) 102:292–7. doi:10.1002/bbpc.19981020303
40. Becerra VM. Solving complex optimal control problems at no cost with PSOPT. In: IEEE international symposium on computer-aided control system design; 2010 Sept 8–10; Yokohama, Japan. Piscataway, NJ: IEEE (2010). p. 1391–6.
Keywords: optimal control, quantum dynamics, semiclassical dynamics, Gaussian wavepackets, proton transfer
Citation: Ramos Ramos AR and Kühn O (2021) Direct Optimal Control Approach to Laser-Driven Quantum Particle Dynamics. Front. Phys. 9:615168. doi: 10.3389/fphy.2021.615168
Received: 08 October 2020; Accepted: 09 February 2021;
Published: 23 April 2021.
Edited by:Tamar Seideman, Northwestern University, United States
Reviewed by:Ilya Averbukh, Weizmann Institute of Science, Israel
Ilia Tutunnikov, Weizmann Institute of Science, Rehovot, Israel, in collaboration with reviewer [IA]
Regina De Vivie-Riedle, Ludwig Maximilian University of Munich, Germany
Copyright © 2021 Ramos Ramos and Kühn. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: O. Kühn, email@example.com