External constraints on optimal control strategies in molecular orientation and photofragmentation: Role of zero-area fields

# External constraints on optimal control strategies in molecular orientation and photofragmentation: Role of zero-area fields

D. Sugny111Laboratoire Interdisciplinaire Carnot de Bourgogne (ICB), UMR 5209 CNRS-Université de Bourgogne, 9 Av. A. Savary, BP 47 870, F-21078 DIJON Cedex, FRANCE, dominique.sugny@u-bourgogne.fr, S. Vranckx222Service de Chimie Quantique et Photophysique, CP 160/09 Université Libre de Bruxelles, B-1050 Brussels, Belgium, M. Ndong, O. Atabek333Institut des Sciences Moléculaires d’Orsay (ISMO), Bât. 350 UMR 8214 CNRS-Université Paris-Sud, Orsay, France, M. Desouter-Lecomte444Laboratoire de Chimie Physique (LCP), UMR 8000 CNRS-Université Paris-Sud, Orsay, France
###### Abstract

We propose a new formulation of optimal and local control algorithms which enforces the constraint of time-integrated zero-area on the control field. The fulfillment of this requirement, crucial in many physical applications, is mathematically implemented by the introduction of a Lagrange multiplier aiming at penalizing the pulse area. This method allows to design a control field with an area as small as possible, while bringing the dynamical system close to the target state. We test the efficiency of this approach on two control purposes in molecular dynamics, namely, orientation and photodissociation.

## 1 Introduction

In recent years, advances in quantum control have emerged through the introduction of appropriate and powerful tools coming from mathematical control theory [1, 2] and by the use of sophisticated experimental techniques to shape the corresponding control fields [3, 4, 5]. All these efforts have lead nowadays to an unexpected and satisfactory agreement between theory and experiment. On the theoretical side, one major tool to design the control field is Optimal Control Theory (OCT) [1, 2]. Over the past few years, numerical iterative methods have been developed in quantum control to solve the optimization problems. Basically, they can be divided into two families, the gradient ascent algorithms [6] and the Krotov [7, 8] or the monotonic [9, 10] ones. The design of optimal control fields by standard iterative algorithms [5] can require the computation of several hundreds numerical propagations of the dynamical quantum system. While the efficiency of this procedure has been established for low dimensional quantum systems, this approach can numerically be prohibitive for large dimensions. In this latter case, it is possible to use more easily accessible numerical methods, such as the local control approach [11, 12]. Roughly speaking, the optimization procedure is built from a Lyapunov function over the state space, which is minimum (or maximum) for the target state. A control field that ensures the monotonic decrease (or increase) of is constructed with the help of the first derivative of . Note that this approach has largely been explored in quantum control [13, 14]. Thanks to the progresses in numerical optimization techniques, it is now possible to design high quality control along with some experimental imperfections and constraints. Recent studies have shown how to extend the standard optimization procedures in order to take into account some experimental requirements such as the spectral constraints [15, 16], the non-linear interaction between the system and the laser field [17, 18], the robustness against experimental errors [19, 20]. In view of experimental applications in quantum control, it is also desirable to design pulse sequences with a zero global time-integrated area. Several works have pointed out the inherent experimental difficulties associated with the use of non-zero area fields [21, 22, 24, 23, 25], in particular for laser fields in the THz regime. Since the dc component of the field is not a solution of Maxwell’s equation, such pulses are distorted when they propagate in free space as well as through focusing optics. The standard optimization procedures do not take into account this basic requirement, designing thus non-physical control fields. In this framework, a question which naturally arises is whether one can adapt the known optimization algorithms to this additional constraint. This paper aims at taking a step toward the answer of this open question by proposing new formulations of optimal and local control algorithms. The zero-area requirement for the laser fields is mathematically fulfilled by the introduction of a Lagrange multiplier and the derivation of the corresponding optimal equations.

The goal of this paper is to explore the efficiency of the generalized optimal and local control algorithms on two key control problems of molecular dynamics: orientation and photodissociation. The remainder of this paper is organized as follows. The new formulations of optimization control algorithms are presented in Sec. 2. Section 3 is devoted to the application of the optimal control algorithm to the enhancement of molecular orientation of CO by THz laser fields at zero temperature. Local control is used in Sec. 4 to manipulate efficiently the photodissociation of HeH. Conclusion and prospective views are given in Sec. 5.

## 2 Formulation of optimal and local control theory with zero-area constraint

### 2.1 The optimal control algorithm

We consider a quantum system whose dynamics is governed by the following Hamiltonian:

 H(t)=H0+E(t)H1 (1)

where is the control field. The state of the system satisfies the differential equation:

 i∂|ψ⟩∂t=H(t)|ψ⟩, (2)

with the initial () state. The units used throughout the paper are atomic units. As mentioned in the introduction, we add the physical constraint of zero-area on the control field:

 A(tf)=∫tf0E(t)dt=0, (3)

where is the control (or total pulse) duration. The goal of the control problem is to maximize the following standard cost functional at time :

 Joc=R[⟨ψf|ψ(t)⟩]−λ∫tf0E(t)2dt, (4)

where is the target state. Other cost functionals could of course be chosen. More specifically, a monotonic optimal control algorithm can be formulated to satisfy the zero-area constraint. For such a purpose let us consider the following cost functional

 Joc=R[⟨ψf|ψ(tf)⟩]−μ[∫tf0E(t)dt]2−λ∫tf0[E(t)−Eref(t)]2/S(t)dt, (5)

where is a reference pulse and an envelope shape given by [26]. The function ensures that the field is smoothly switched on and off at the beginning and at the end of the control. The positive parameters and , expressed in a.u., weight the different parts of , which penalize the area of the field (- term) and its energy (- term). In this algorithm, we determine the field at step from the field at step , such that the variation . At step , the reference field is taken as and we denote by its time-integrated area. Note that the choice leads to a smooth transition of the control field between two iterations of the algorithm [27, 28].

The variation can be expressed as follows:

 ΔJoc=∫tf0(Ek+1−Ek)[I[⟨χk|H1|ψk+1⟩]−λ(Ek+1−Ek)/S(t)−2μAk(Ek+1−Ek)]dt, (6)

is obtained from backward propagation of the target taken as an intial state for Eq.(2) and we assume that . One deduces that the choice:

 Ek+1−Ek=S(t)I[⟨χk|H1|ψk+1⟩]λ−(Ek+1−Ek)−2μλAkS(t) (7)

ensures a monotonic increase of . This leads to the following correction of the control field:

 Ek+1=Ek+S(t)I[⟨χk|H1|ψk+1⟩]2λ−μλS(t)Ak. (8)

The monotonic algorithm can thus be schematized as follows:

1. Guess an initial control field, which at step is taken as .

2. Starting from , propagate backward with

3. Evaluate the correction of the control field as obtained from Eq. (8), while propagating forward the state of the system from with .

4. With the new control, go to step 2 by incrementing the index by 1.

### 2.2 Local control

The same formalism can be used to extend the local control approach to the zero-area constraint. We refer the reader to [11] for a complete introduction of this method. Using the same notations as in the previous section, we consider as Lyapunov function the cost functional:

 Jlc(t)=⟨ψ(t)|O|ψ(t)⟩−μA(t)2, (9)

where is any operator such that . The cost is maximum at a given time if is maximum and . The computation of the time derivative of gives:

 ˙Jlc(t)=E(t)⟨ψ(t)|[O,H1]i|ψ(t)⟩−2μA(t)E(t). (10)

One can choose:

 E(t)=ε(⟨ψ(t)|[O,H1]i|ψ(t)⟩−2μA(t)), (11)

to ensure the monotonic increase of from the initial state of the system . The parameter is a small positive parameter which is adjusted to limit the amplitude of the control field. In order to illustrate the efficiency of the algorithm with the zero-area constraint, we introduce two measures:

 (12)

where and are the optimized field with and without the zero-area constraint, respectively.

## 3 Molecular orientation dynamics with zero-area fields

### 3.1 Model system

In this section, we present the theoretical model for the control of molecular orientation by means of THz laser pulses [29, 30, 31, 32, 33]. We consider the linear CO molecule described in a rigid rotor approximation and driven by a linearly polarized electric field along the - axis of the laboratory frame. The molecule is assumed to be in its ground vibronic state. The Hamiltonian of the system is given by:

 H(t)=BJ2−dcosθE(t). (13)

The first term is the field-free rotational Hamiltonian, where is the rotational constant of the molecule and the angular momentum operator. The second term represents the interaction of the system with the laser field, the parameter being the permanent molecular dipole moment. The spatial position of the diatomic molecule is given in the laboratory frame by the spherical coordinates . Due to the symmetry of revolution of the problem with respect to the - axis, the Hamiltonian depends only on , the angle between the molecular axis and the polarization vector. Numerical values of the molecular parameters are taken as and . In the case of zero rotational temperature () which is assumed hereafter, the system is described at time by the wave function and its dynamics by the time-dependent Schrödinger equation (Eq. (2)). The dynamics is solved numerically using a third-order split operator algorithm [35]. We finally recall that the degree of orientation of the molecular system is given by the expectation value .

### 3.2 Orientation dynamics

In this paragraph, we test the new formulation of the optimal control algorithm presented in Sec. 2.1 on the specific example of CO orientation dynamics. The initial state is the ground state denoted in the spherical harmonics basis set . Due to symmetry, the projection of the angular momentum on the axis is a good quantum number. This property leads to the fact that only the states , , will be populated during the dynamics. In the numerical computations, we consider a finite Hilbert space of size . Note that this size is sufficient for an accurate description of the dynamics considering the intensity range of the laser field used in this work. To define the cost functional, rather than maximizing the expectation value , we focus on a target state maximizing this value in a finite-dimensional Hilbert space spanned by the states with . The details of the construction of can be found in Refs. [36, 37]. Here, the parameter is taken to be . The guess field is a Gaussian pulse of 288 fs full width at half maximum, centered at . The optimization time, , is taken to be the rotational period of the molecule. The parameter defined in the cost functional is fixed to 100. Both, with and without zero-area constraint, the optimized field achieves a fidelity higher than 99. The results are shown in Fig.1. Note that, since the ratio has the dimension of time, we found convenient to plot the evolution of the measures with respect to (panel (b)). The upper right panel, Fig.1(a), shows the deviation from unity of the fidelity (or final objective)

 Ctf=|⟨ψ(tf)|ψf⟩|2 (14)

as a function of the number of iterations.

As could be expected, optimization with zero-area constraint requires a large number of iterations to reach a high fidelity. For , i.e. without zero-area constraint, the final cost reaches a value of 0.99 after 30 iterations, while in the case for which , this value is reached after 40 iterations. The variations of the and with respect to the parameter are shown in Fig.1(b). For , the area of the optimized field obtained with the zero-area constraint is two orders of magnitude smaller than the one obtained without the zero-area constraint. Note that the values of and displayed in Fig.1(b) are computed with optimized fields which achieve a fidelity higher than 99. Figure 1(c) compares the guess (reference) field and the optimized fields with and without the zero-area constraint. The optimized field with the zero-area constraint shown in Fig. 1(c) is obtained with . The optimized fields are similar to the guess field at the beginning of the control and become slightly different for . Figure 1(d) corresponds to the time evolution of during two rotational periods (one with the field and one without). It compares the dynamics of induced by the three fields shown in the lower left panel. Clearly, the dynamics with the guess field is very different from the one obtained with the optimized fields, which have similar features.

## 4 Photodissociation

### 4.1 Model system

We consider the photodissociation of HeH through the singlet excited states coupled by numerous nonadiabatic interactions. A local control strategy has recently been applied to find a field able to enhance the yield in two specific dissociation channels either leading to He, or to He [38]. The diatomic system is described by its reduced mass and the internuclear distance . Due to the short duration of the pulses (as compared with the rotational period) a frozen rotation approximated is valid. In addition, the molecule is assumed to be pre-aligned along the -direction of the laboratory frame. The initial state is the vibrationless ground electronic adiabatic state. Dynamics is performed in the diabatic representation. The total Hamiltonian is written as:

 H=H0−N∑i,j=1Mdij(R)E(t)=N∑i,j=1|φdi⟩(Tijδij+Vdij(R))⟨φdj|−N∑i,j=1Mdij(R)E(t) (15)

where is the field free Hamiltonian, is the number of diabatic electronic states under consideration and . All the adiabatic potential energy curves and the nonadiabatic couplings have been computed in Ref. [39]. are the diabatized dipole transition matrix elements. We shall consider only parallel transitions among the singlet states induced by the dipole operator by assuming the internuclear axis pointing along the -direction. The adiabatic-to-diabatic transformation matrix has been derived by integrating from the asymptotic region where both representations coincide.

### 4.2 Dissociation control

The local control of a nonadiabatic dissociation requires a careful choice of the operator referred to in Eq. (9) since it has to commute with the field free Hamiltonian [38]. In the nonadiabatic case, the projectors on either adiabatic or diabatic states are thus not relevant as they do not commute with this Hamiltonian due to kinetic or potential couplings respectively. This crucial problem can be overcome by using projectors on eigenstates of , i.e. on scattering states correlating with the controlled exit channels. In this example, the operator takes the form:

 O=∑p∈S∫dk|ξ−p(k)⟩⟨ξ−p(k)| (16)

where represents the two channels leading to the target He fragments, the objective being . The local control field now reads

 E(t)=ηI(∑p∈S∫dk⟨ψ(t)|ξ−p(k)⟩⟨ξ−p(k)|dz|ψ(t)⟩)−2μA(t) (17)

involving two adjustable parameters and .

The ingoing scattering states are estimated using a time-dependent approach based on Møller-operators [40] as derived in Ref. [38]. This control strategy remains local in time but can preemptively account for nonadiabatic transitions that occur later in the dynamics. The photodissociation cross section [41] shows that there is no spectral range where the He dissociation channels dominate. The local control finds a very complicated electric field which begins by an oscillatory pattern followed after 10 fs by a strong complex positive shape whose area is obviously not zero (see Fig. 2(a)). We therefore choose this example to check the efficiency of the zero-area constraint algorithm. We use 4.2 and different values of to assess the effect of the zero-area constraint on the control field. Figure 2(a) shows the dramatic evolution of the pulse profile for two values of (0.05 and 0.10). The main change is the suppression of the erratic positive structure which can be interpreted as a Stark field. However, as shown in Fig. 2(c) the variation of is not monotonic and the best correction is actually obtained for a given value of around , leading to almost zero value for . Close to the optimal , one observes a clear improvement of the zero-area constraint, but at the price of a decrease of the objective (from without constraint to with ) as shown in Fig. 2(b). However, a brute force strategy by filtration of near-zero frequencies already provides a large correction to the non vanishing area while also decreasing the yield. The field with and after filtration of the low frequencies is shown in black line in Fig. 2(d) and the corresponding yield is given in Fig. 2(e). The objective is divided by a factor 2. We also filter the fields generated by the constrained control with equal to 0.05 and 0.10 which improves even further the zero-area constraint. The measures evaluated with the filtered fields (Fig. 2(f)) are smaller than the unfiltered ones by about two orders of magnitude, while the objectives are reduced by about a factor 1.5, as seen by comparing Fig. 2(e) and Fig. 2(b).

## 5 Conclusion

We have presented new formulations of optimization algorithms in order to design control fields with zero area. The procedure is built from the introduction of a Lagrange multiplier aiming at penalizing the area of the field. The value of the corresponding parameter is chosen to express the relative weight between the area and the other terms of the cost functional, i.e. the projection onto the target state and the energy of the electric field. The efficiency of the algorithms is exemplified on two key control problems of molecular dynamics: namely, the molecular orientation of CO and the photodissociation of HeH.

From a theoretical view point, the zero-area constraint is crucial since all electromagnetic fields which are physically realistic must fulfill this requirement. Surprisingly, this physical requirement has been considered, up to date, by only very few works in quantum control. As such, this work can be viewed as a major advance of the state-of-the-art of optimization algorithms that could open the way to promizing future prospects.

## 6 Acknowledgments

S. V. acknowledges financial support from the Fonds de la Recherche Scientifique (FNRS) of Belgium. Financial support from the Conseil Régional de Bourgogne and the QUAINT coordination action (EC FET-Open) is gratefully acknowledged by D. S. and M. N..

## References

• [1] Pontryagin, L. et al. Mathematical theory of optimal processes; Mir, Moscou, 1974.
• [2] Bryson, A. E.; Ho, Y.-C. Applied optimal control; Taylor & Francis, New York London, 1975.
• [3] Rice, S.; Zhao, M. Optimal control of molecular dynamics; Wiley, New York, 2003.
• [4] Shapiro, M.; Brumer, P. Principles of quantum control of molecular processes; Wiley, New York, 2003.
• [5] Brif, R.C.C.; Rabitz, H. New. J. Phys. 2010, 12 075008.
• [6] Skinner, T. E.; Reiss, T. O.; Luy, B.; Khaneja, N.; Glaser, S. J. J. Magn. Reson. 2003, 163 8.
• [7] Somlói, J.; Kazakovski, V. A.; Tannor, D. J. Chem. Phys. 1993, 172 85.
• [8] Reich, D. M.; Ndong, M.; Koch, C. P. J. Chem. Phys. 2012, 136 104103.
• [9] Zhu, W.; Botina, J.; Rabitz, H. J. Chem. Phys. 1998, 108 1953.
• [10] Maday, Y.; Turinici, G. J. Chem. Phys. 2003, 118 8191.
• [11] Engel, V.; Meier, C.; Tannor, D. Adv. Chem. Phys. 2009, 141 29.
• [12] D’Alessandro, D. Introduction to quantum control and dynamics; Boca Raton, FL, CRC Press 2007.
• [13] Sugawara, M. J. Chem. Phys. 2003, 118 6784.
• [14] Wang, X.; Schirmer, S. G.; Phys. Rev. A 2009, 80 042305.
• [15] Gollub, C.; Kowalewski, M.; Vivie-Riedle R. Phys. Rev. Lett. 2008, 101 073002.
• [16] Lapert, M.; Tehini, R.; Turinici, G.; Sugny, D. Phys. Rev. A 2009, 79 063411.
• [17] Ohtsuki, Y.; Nakagami, K. Phys. Rev. A 2008, 77 033414.
• [18] Lapert, M.; Tehini, R.; Turinici, G.; Sugny, D. Phys. Rev. A 2008, 78 023408.
• [19] Kobzar, K.; Skinner, T. E.; Khaneja, N.; Glaser, S. J.; Luy, B. J. Magn. Reson. 2004, 170 236.
• [20] Zhang, Y.; Lapert, M.; Braun, M.; Sugny, D.; Glaser, S. J. J. Chem. Phys. 2011, 134 054103.
• [21] You, D.; Bucksbaum, P. H. J. Opt. Soc. Am. B 1997, 14 1651.
• [22] Liao, S.-L.; Ho, T.-S.; Rabitz, H.; Chu, S.-I. Phys. Rev. A 2013, 87 013429.
• [23] Ortigoso, J. J. Chem. Phys. 2012, 137 044303.
• [24] Sugny, D.; Keller, A.; Atabek, O.; Daems, D.; Guérin, S.; Jauslin, H. R. Phys. Rev. A 2004, 69 043407.
• [25] Lapert, M.; Sugny, D. Phys. Rev. A 2012, 85 063418.
• [26] Sunderman, K.; de Vivie-Riedle, R. J. Chem. Phys. 1999, 110 1896.
• [27] Palao, J. P.; Kosloff, R. Phys. Rev. Lett. 2002, 89 188301.
• [28] Bartana, A.; Kosloff, R.; Tannor, D. J. Chem. Phys. 2001, 267 195.
• [29] Seideman, T.; Hamilton, E. Adv. At. Mol. Opt. Phys. 2006, 52 289.
• [30] Stapelfeldt, H.; Seideman, T. Rev. Mod. Phys. 2003, 75 543.
• [31] Shu, C. C.; Yuan, K. J.; Hu, W. H.; Cong S. L. J. Chem. Phys. 2010, 132, 244311.
• [32] Fleischer, S.; Zhou, Y.; Field, R. W.; Nelson, K. A. Phys. Rev. Lett. 2011, 107 163603.
• [33] Shu, C.-C.; Henriksen, N. E. Phys. Rev. A 2013, 87 013408.
• [34] Viellard, T.; Chaussard, F.; Sugny, D.; Lavorel, B.; Faucher, O. J. Raman Spec. 2008, 39 694.
• [35] Feit, M. D.; Fleck, J. A.; Steiger, A. J. Comput. Phys. 1982, 47 412.
• [36] Sugny, D.; Keller, A.; Atabek, O.; Daems, D.; Dion, C. M.; Guérin, S.; Jauslin, H. R. Phys. Rev. A 2005, 71 063402.
• [37] Sugny, D.; Keller, A.; Atabek, O.; Daems, D.; Dion, C. M.; Guérin, S.; Jauslin, H. R. Phys. Rev. A 2005, 72 032704.
• [38] Bomble, L.; Chenel, A.; Meier, C.; Desouter-Lecomte M. J. Chem. Phys. 2011, 134, 204112.
• [39] Loreau, J.; Palmeri, P.; Quinet, P.; Liévin, J.; Vaeck, N. J. Phys. B: At. Mol. Opt. Phys. 2010, 43 065101.
• [40] Tannor, D. J.; Quantum Mechanics: A Time Dependent Perspective 2007, University Science Book: Sausalito, CA.
• [41] Sodoga, A.; Loreau, J.; Lauvergnat, D.; Justum, Y.; Vaeck, N.; Desouter-Lecomte, M. Phys. Rev. A 2009, 80 033417.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters