Quantum optimal quantum control field design using logarithmic maps

16 November 2001

Chemical Physics Letters 348 (2001) 440±446 www.elsevier.com/locate/cplett

Quantum optimal quantum control ®eld design using logarithmic maps Julie S. Biteen, J.M. Geremia, Herschel Rabitz

*

Department of Chemistry, Princeton University, 121 Frick Laboratory, Washington Road and William St., Princeton, NJ 08544-1009, USA Received 23 May 2001

Abstract A mapping technique is introduced to expand the capabilities of current control ®eld design procedures. The maps relate the control ®eld to the logarithm of the time evolution operator. Over the dynamic range of the maps, which begins at the sudden limit and extends beyond, they are found to be more accurate than ®eld ! observable maps. The maps may be used as part of an iterative computational algorithm for ®eld design. This process is illustrated for the design of a ®eld to meet a population transfer objective. Ó 2001 Elsevier Science B.V. All rights reserved.

1. Introduction The control of quantum dynamics phenomena has many potential applications in chemistry, physics, and engineering [1±3]. Control is commonly achieved using tailored laser ®elds in order to manipulate a system's Hamiltonian [4±6]. Modulated laser ®elds have a rich structure that can induce subtle molecular dynamical behavior to discriminate amongst the di€erent possible control outcomes. Optimal control theory [7] provides a general framework for designing laser ®elds [5,6], however, solving the design equations can entail a large computational e€ort. Closed-loop learning procedures [1,2,8±16], operating directly in the

*

Corresponding author. Fax: +1-609-258-0967. E-mail address: [email protected] (H. Rabitz).

laboratory, provide the means to identify the control ®elds that best achieve a target objective. Such experiments could bene®t from prior designs as an initial start and also as a means for providing physical insight into the control process. The present Letter aims to expand the control ®eld design capabilities through the introduction of a special ®eld ! output mapping technique. Recently, new mapping algorithms have been proposed capable of achieving control designs [17,18]. The ®rst stage of such an algorithm involves constructing a map between a set of laser ®eld variables and their outcome by sampling a representative set of ®elds. The next stage is an optimization performed by searching the map for ®elds that best meet the control objective. In principle, no further solutions of the Schr odinger equation need to be obtained after constructing the map, provided that the map is suciently accurate and includes in its domain one or more

0009-2614/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 0 9 - 2 6 1 4 ( 0 1 ) 0 1 1 4 4 - 7

J.S. Biteen et al. / Chemical Physics Letters 348 (2001) 440±446

441

control ®elds that can adequately achieve the desired goal. As searching a map can be very fast, global techniques (e.g., genetic algorithms (GAs) [19,20]) can be employed to thoroughly explore for viable controls. If the control objective is altered, often no additional solutions of the Schr odinger equation are needed since a search for the control ®eld that reaches this new goal can be conducted using previously constructed maps. Furthermore, the maps can be used to attain control mechanism information, including the relationship between the di€erent control variables [17]. Several mapping procedures have been introduced, and any of these could be used in a mapfacilitated learning algorithm. Procedures have been developed that generate a series of maps to approximate the ®eld ! observable relationship as linear over their domain [21], but each map is only accurate over a small region. Nonlinear mapping techniques provide a more e€ective model for the relationship between a ®eld and its control outcome. Alternative forms of these techniques have been proposed to operate in the frequency [17] and time [18] domains. Here, we present a new mapping technique that could operate in either of the latter domains while taking advantage of the special form of the quantum mechanical time evolution operator, U. Rather than mapping directly to quantum observables or to U itself, the present algorithm maps the electric ®eld to ln U . Treatment of observables for ®eld optimization can then be readily performed with the ln U map. Section 2 describes the Magnus exponential form of U, lays out the nonlinear technique used to construct the maps, and shows how such maps may be used in control problems. Section 3 presents a series of simulations to explore the reliability of the ln U maps and demonstrate their capabilities in control problems. Some brief conclusions are given in Section 4.

®eld ! output relationship is learned by sampling a representative set of ®elds in a speci®ed domain and recording the output dynamical response of the system. During the optimization stage, the map domain is explored by interpolation of these results to identify one or more ®elds that best achieve the control target. Such a search, which requires no further solutions of the Schr odinger equation, can be conducted using any suitable algorithm. GAs [19] are attractive because they have global search capabilities and propagate a large solution family in an attempt to identify multiple control ®elds. The map-facilitated control concept is general and can incorporate any suitable mapping technique. Previous studies [17,18] utilized ®eld ! observable maps, while the present technique is based on the use of field ! ln U maps to exploit the expected faster convergence of ln U with respect to the ®eld variables. The ®eld, E, may be expressed either in time or frequency domain variables for mapping purposes.

2. Field mapping

where T is the time ordering operator. An alternative expression for the solution to Eq. (2) was introduced by Magnus [22], where

The map-facilitated control algorithm proposed in this Letter is accomplished in two stages: (1) map construction and (2) control ®eld optimization. In the map construction stage, the

2.1. The field ! ln U relationship A quantum mechanical system is propagated from its initial state, W…t0 †, at time, t0 , by the time evolution operator, U …t; t0 †, such that, W…t† ˆ U …t; t0 †W…t0 †:

…1†

The operator, U …t; t0 †, satis®es ih

oU …t; t0 † ˆ H …t†U …t; t0 †; ot

…2†

where H …t† ˆ H0

lE…t†:

…3†

Here, H0 is the Hamiltonian of the unperturbed system, and l is the transition dipole operator. U …t; t0 † can be formally expressed as   Z i t H …t0 † dt0 ; U …t; t0 † ˆ T exp …4† h t0

U …t; t0 † ˆ exp ‰ M…t; t0 †Š with M ˆ ln U , and

…5†

442

J.S. Biteen et al. / Chemical Physics Letters 348 (2001) 440±446

i h 

M…t; t0 † ˆ

1 2 ‡

1 4 Z



Z t0

Z Z

t

t t0 t t0

t00 t0

dt0 H …t0 †

dt00 dt000

Z

t00

t0

Z

t0

dt0 ‰H …t0 †; H …t00 †Š

t000

dt00

 0

dt ‰H …t0 †; H …t00 †Š; H …t000 †



Z t Z t000 1 dt000 dt00 12 t0 t0 ) Z t00   0 0 00 000  dt H …t †; ‰H …t †; H …t †Š ‡


: ‡

t0

…6† Previous studies have indicated that the series for M may have rapid convergence [23], especially under conditions of either a weak ®eld or the near sudden regime where the higher commutators diminish rapidly [24]. These observations suggest that E ! ln U maps may have simpler behavior than maps of E ! U or E ! observables. The actual maps will be constructed for ln U …T ; t0 †, where T is the ®nal time at which control is desired. 2.2. The mapping algorithm In order to construct the map, E…t† is described in terms of a set of N control variables, x ˆ fx1 ; . . . ; xN g. The assignment of these variables can be chosen on convenience or physical grounds such that, E…t† ) E…x; t†

…7†

which re-expresses ln U …T ; t0 † as a function, f …x†, of x, ln U …T ; t0 † ' f …x†:

…8†

With this transformation made, the mapping procedure can begin with a choice of an appropriate input variable domain, X, where max Xi ˆ ‰xmin Š speci®es the range of each variable, i ; xi xi . For control ®elds within this domain, f …x† may be considered an operational substitute for ln U …T ; t0 †.

Map construction is done by discretely sampling ®elds in X. For each sample x, the corresponding f …x† is found by numerical integration of Eq. (2) followed by the formation of ln U …T ; t0 †. The ®eld sampling must be done to a resolution that is sucient for eventual map interpolation over X. However, sampling f …x† in the traditional sense, at a rate of S points per variable, xi , for N dimensions, leads to computational complexity growing as O…S N †. The number of ®eld variables, N, is typically large, so it becomes impractical to map the system in such a way. This curse of dimensionality problem can be overcome by using an alternative high dimensional model representation (HDMR) [17,18,25±33] technique. With HDMR, the output, f …x†, is expanded as a hierarchy of functions of increasing dimensions, and HDMR has already been applied to control problems [17,18]. The HDMR expansion is given by f …x† ˆ f0 ‡

N X iˆ1

fi …xi † ‡

N X

fij …xi ; xj † ‡




i
‡ f1;...;N …x1 ; . . . ; xN †:

…9†

The hierarchy of the terms in Eq. (9) expresses the increasing degrees of cooperativity amongst all the variables, x1 ; . . . ; xN , in f …x†. The L ˆ 0 term, f0 , is the nominal value of f …x† over its domain, the L ˆ 1 term, fi …xi †, is a one-dimensional function that describes the e€ect of the single variable, xi , on the output, and the L ˆ 2 term, fij …xi ; xj †, is a twodimensional function that describes the cooperative e€ect of the pair of variables, xi and xj , etc. The functions increase in dimension, L, up to a ®nal L ˆ N term, f1;...;N …x1 ; . . . ; xN †, which accounts for any non-separable overall cooperativity between all the variables. The degree of cooperativity found in many real systems with well-chosen variables has been observed to be low order (L 6 3) [17,18,26,28±33]. With this convergence behavior, the truncation of Eq. (9) at L  N allows for accurate representation by a map whose construction cost grows only polynomically with N. The total number of sample points needed is O…SN † for a ®rst order (L ˆ 1) map, O…S 2 N 2 † for a second order (L ˆ 2) map, etc.

J.S. Biteen et al. / Chemical Physics Letters 348 (2001) 440±446

There are several formulations of HDMR which are formally equivalent [25±27]. The maps in this Letter are constructed using the cut-center approach [17,18,33] to HDMR, with the component functions constructed as: f0 ˆ f …x1 ; . . .; xN †;

…10a†

fi …xi † ˆ f …x1 ; . . .; xi ; . . .; xN † f0 ; fij …xi ; xj † ˆ f …x1 ; . . .; xi . . .; xj ; . . .; xN †

…10b†

fi …xi †

…10c†

fj …xj †

f0 ;

.. .  is a point in the N-dimensional space where x called the cut-center, and . . . implies that the corresponding variables are taken at their cut-center values. Further details on map construction can be found in the literature [17,18,25±33]. 2.3. Relating logarithmic maps to control problems The quantum control observable, U…T †, at the target time, T, is given by, U…T † ˆ hW…T †jOjW…T †i;

…11†

where O is the Hermitian operator associated with the observable being targeted. Equation (11) may be re-expressed in terms of f …x† from Eq. (8) by noting that jW…T †i ' exp‰f …x†ŠjW…t0 †i, producing, U…T † ' g‰f …x†Š ˆ h exp‰f …x†ŠW…t0 †jOj exp‰f …x†ŠW…t0 †i:

…12†

Optimal ®elds can now be sought by searching the map over X to ®nd the point, x 2 X, for which g‰f …x †Š is as close as possible to the target goal, U…target† . This task is accomplished by minimizing the cost function,

2 …13† J…x† ˆ U…target† g‰f …x†Š over x 2 X. J…x† can also be extended to optimize other features (e.g., smoothness of the ®eld, robustness of the solution, etc.) beyond simply achieving the target observable [34]. For the purposes of this Letter, we use only the simple cost function in Eq. (13). The overall learning algorithm involves closedloop iterations drawing on the map output, f …x†,

443

to minimize the cost function, J…x†. Once an HDMR map is constructed, the GA seeks, over a series of generations, to discover an optimal ®eld or family of ®elds for a given target by choosing and assessing a set of pulses at each generation. When the optimization routine converges, the best ®eld, x , from the map is tested against the truth (i.e., a new solution of the Schr odinger equation), and if its output is within an acceptable range of the target, then E…x ; t† is optimal. Otherwise, the algorithm recommences. A sequence of quasi-local maps may be the most ecient way of performing design with the algorithm [17,18]. 3. Illustrations In order to demonstrate the utility of the procedure discussed in Section 2, E ! ln U maps were compared in terms of their accuracy to the previously proposed E ! U…T † maps [17], and assessed for their utility in control problems. In both cases, E was expressed in the frequency domain variables described below. The maps were constructed for a series of simple systems with four nondegenerate energy levels and ®ve possible resonance frequencies. These systems evolved under the in¯uence of an electric ®eld, E…x; t†, de®ned in terms of a modulated Gaussian pulse of the form,   5  t2 X E…x; t† ˆ exp Aj sin…xj t† 2 2r jˆ1  ‡ Bj cos…xj t† ; …14† where Aj and Bj are the sine and cosine Fourier amplitudes corresponding to each of the ®ve resonance frequencies, xj , of H0 . The N ˆ 10 control variables are de®ned as, x ˆ A [ B ˆ fA1 ; . . . ; A5 ; B1 ; . . . ; B5 g. The pulse had a total length of 1.0 ps, and a width, r ˆ 0:15 ps. 3.1. Map accuracy A series of nondegenerate four-state systems, related by an energy level scaling factor, a, was considered. The resonance frequencies and transition dipole matrix elements for the a ˆ 1 system

444

J.S. Biteen et al. / Chemical Physics Letters 348 (2001) 440±446

are given by: x12 ˆ 2:3703  1012 , x13 ˆ 3:7925 1012 , x23 ˆ 1:4222  1012 , x24 ˆ 2:3703  1012 , and x34 ˆ 0:9481  1012 rad/s; l12 ˆ 10:3709 10 32 , l13 ˆ 2:1581  10 32 , l14 ˆ 0, l23 ˆ 14:1584  10 32 , l24 ˆ 3:6457  10 32 , and l34 ˆ 16:7058  10 32 C m. The other systems in the series had the same transition dipoles, but the system energy levels were scaled by a over the range 0.25±3.00 in increments of 0.25. For each system, a series of ®rst order HDMR maps with domain windows, Dxi ˆ xmax xmin i i , of  were constructed about a cut-center ®eld, 0.5 V/A x1 ˆ x2 ˆ


ˆ x10 ˆ E0 , where the value of E0 for  in increeach map varied from 0.5 to 12.5 V/A  ments of 0.5 V/A. For each of the 300 combinations of a and E0 , two maps were constructed: E ! ln U and E ! U…T †. For the latter maps, fU…T † …x† ' U…T † was chosen to be the fourth state ®nal population.

(a)

The output of the former maps, fln U …x† ' ln U …T ; t0 †, was converted by the transformation in Eq. (12) such that g‰fln U …x†Š gave an equivalent population measurement. Each map was tested by 1000 random ®elds in its domain. For each test ®eld, the fourth state ®nal population according to the map was compared against the true ®nal population of that state found from another solution of the Schr odinger equation. The measure of the E ! ln U ‰or E ! U…T †Š map quality, Eln U (or EU…T † ), was de®ned as the fraction of the tests which reproduced the true population within a 10% relative error. The results of the tests on the E ! ln U maps are illustrated in a plot of Eln U vs. a and E0 in Fig. 1a. Each square in this ®gure corresponds Eln U for a single map, with darker shading denoting higher values of Eln U (i.e., more accurate maps). The map accuracy was observed to improve as a decreased.

(b)

Fig. 1. Measure of the accuracy of the E ! ln U maps and comparison to E ! U…T † map accuracy. The squares in (a) indicate Eln U , the fraction of 1000 test ®elds in each E ! ln U map domain that had  10% relative error compared to the truth, as a function of the energy scaling factor, a, related to the magnitude of the resonance frequencies of the system, and the amplitude of the Fourier components of the cut-center ®eld, E0 . The squares in (b) show DE, the di€erence in accuracy between the E ! ln U and E ! U…T † maps. Darker shading represents higher values (i.e., better accuracy for the E ! ln U maps), lighter shading represents lower values, and hashed squares indicate negative values.

J.S. Biteen et al. / Chemical Physics Letters 348 (2001) 440±446

The di€erence in the map quality for the two techniques, DE ˆ Eln U EU…T † , was then calculated for each combination of a and E0 . A plot of DE vs. a and E0 is shown in Fig. 1b. The shaded region in Fig. 1b indicates the …a; E0 † combinations for which the E ! ln U maps are more accurate than the E ! U…T † maps. Darker shading corresponds to maps with Eln U  EU…T † and lighter shading shows maps for which DE was still positive, but smaller. The hashed squares in the ®gure indicate …a; E0 † combinations for which Eln U < EU…T † . The DE > 0 case, showing better accuracy for the E ! ln U maps, was observed over a sizeable domain. In particular, the behavior of Eln U , which increases at low a, is not mirrored by EU…T † . As a result, DE is large for small values of a. The increased accuracy of the E ! ln U maps at low a occurs because the Magnus relation outlined in Eq. (6) dictates that the HDMR maps should be exact to ®rst order at the sudden (a ˆ 0) limit. In order to ascertain that these maps were useful beyond the sudden limit, the fourth state population of the cut-center ®eld, as found from the sudden approximation, was compared to the fourth state population found by the map for each map in Fig. 1a,b. The range over which the E ! ln U maps were more accurate than the E ! U…T † maps (i.e., the range over which DE > 0) extends well beyond the sudden region, which is a K 0:25 in the present case. 3.2. Control design with field ! ln U maps The control algorithm [17,18] discussed in Section 2.3 was applied to a population transfer problem using the E ! ln U maps. The goal was to ®nd an optimal ®eld that drove the population from its initial distribution, p0 ˆ ‰1; 0; 0; 0Š, to a target distribution, pT ˆ [0, 0, 0, 1], i.e., to move all population from j1i to j4i. This was done by seeking the ®eld that minimized the cost function, J…x† ˆ kpT

g‰f …x†Šk2 :

…15†

Applying this algorithm to the a ˆ 1:0 system and beginning with a cut-center ®eld at E0 ˆ 1:0 V/  the algorithm successfully found ®elds that A,

445

transferred > 80% of the population into the target state. This control procedure was achieved utilizing a series of smaller maps for improved accuracy and eciency [17,18]. The map domain was readjusted to be centered about the previous domain's best ®eld, and a ®nal population of 0.804 was achieved after the construction of only six maps. 4. Discussion This Letter introduced a new technique for mapping the quantum dynamical e€ects of a laser ®eld acting on a system. The E ! ln U maps were shown to predict the observable more accurately than the E ! U…T † counterpart over a signi®cant regime that extended well beyond the sudden limit. The E ! ln U mapping technique was used to create an e€ective substitute for further calculations on a system. The maps could be thoroughly searched by a global optimization algorithm to discover ®elds that approached the desired objective. As discussed in Section 2.3, the control algorithm used in this Letter is an iterative process. At the end of each optimization, the output of the best ®eld is tested against the truth. If the true output is not within a tolerable range of the target, a new map is constructed and the next iteration begins. This iterative process is mainly intended for the purpose of ®nding an optimal ®eld for any target regardless of whether a good ®eld exists within the initial domain, but it carries an additional bene®t: the iterations will tend to drive the maps toward areas of greater accuracy. In the present mapping procedure, some error is introduced by the truncation of the HDMR expansion in Eq. (9) and by the map interpolation which may be sensitive to phase changes in ln U . However, the algorithm has broad tolerance for such errors, as its iterative nature overcomes them by operating in the reliable domain of a series of maps. The optimal ®elds designed using the algorithm in this Letter would be best suited for attaining an estimate of the optimal ®eld to serve as a starting point for laboratory learning control procedures [2,8±16].

446

J.S. Biteen et al. / Chemical Physics Letters 348 (2001) 440±446

Acknowledgements The authors acknowledge support from the NSF and the DOD. References [1] H. Rabitz, R. de Vivie-Riedle, M. Motzkus, K. Kompa, Science 288 (2000) 824. [2] R. Judson, H. Rabitz, Phys. Rev. Lett. 68 (1992) 1500. [3] D. Tannor, S. Rice, J. Chem. Phys. 83 (1985) 5013. [4] S. Shi, A. Woody, H. Rabitz, J. Chem. Phys. 11 (1988) 6870. [5] A. Pierce, M. Dahleh, H. Rabitz, Phys. Rev. A 37 (1988) 4950. [6] R. Koslo€, S. Rice, P. Gaspard, S. Tersigni, D. Tannor, J. Chem. Phys. 139 (1989) 201. [7] H. Rabitz, W. Zhu, Accounts Chem. Res. 33 (2000) 572. [8] A. Assion, M. Baumbert, T. Bergt, B. Brixner, V. Kiefer, M. Seyfried, M. Strehle, G. Gerber, Science 282 (1998) 919. [9] C. Bardeen, V. Yakovlev, K. Wilson, S. Carpenter, P. Weber, W. Warren, Chem. Phys. Lett. 280 (1997) 151. [10] T. Feurer, A. Glass, T. Rozgonyi, R. Sauerbrey, G. Szabo, Chem. Phys. 267 (2001) 223. [11] T. Weinacht, J. White, P. Bucksbaum, J. Phys. Chem. A 103 (1999) 10166. [12] R. Bartels, B. Backus, E. Zeek, L. Misoguti, G. Vdovin, I. Christov, M. Murnane, H. Kapteyn, Nature 406 (2000) 164. [13] G. Toth, A. Lorincz, H. Rabitz, J. Chem. Phys. 101 (1994) 3715. [14] P. Gross, D. Neuhauser, H. Rabitz, J. Chem. Phys. 98 (1993) 4557.

[15] S. Vajda, A. Bartelt, E.-C. Kaposta, T. Leisner, C. Lupulescu, S. Minemoto, P. Rosendo-Francisco, L. W oste, Chem. Phys. 267 (2001) 231. [16] R. Levis, G. Menkir, H. Rabitz, Science 292 (2001) 709. [17] J. Geremia, E. Weiss, H. Rabitz, Chem. Phys. 267 (2001) 209. [18] J. Biteen, J. Geremia, H. Rabitz, Chem. Phys., in press (2001). [19] D. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wesley, Reading, MA, 1989. [20] S. Forrest (Ed.), Proceedings of the Fifth International Conference on Genetic Algorithms, Morgan Kaufmann, San Mateo, 1993. [21] M. Phan, H. Rabitz, J. Chem. Phys. 217 (1997) 389. [22] W. Magnus, Commun. Pure Appl. Math. 7 (1954) 649. [23] P. Pechukas, J. Light, J. Chem. Phys. 44 (1966) 3897. [24] L. Smith, S. Augustin, H. Rabitz, J. Comp. Phys. 45 (1982) 417. [25] O. Alis, H. Rabitz, J. Math. Chem. 25 (1999) 197. [26] H. Rabitz, O. Alis, J. Shorter, K. Shim, Comp. Phys. Comm. 117 (1999) 11. [27] H. Rabitz, O. Alis, in: Sensitivity Analysis, Wiley, Chichester, 2000, p. 199. [28] J. Shorter, P. Ip, H. Rabitz, Geophys. Res. Lett. 27 (2000) 3485. [29] J. Shorter, P. Ip, H. Rabitz, J. Phys. Chem. A 36 (1999) 7192. [30] H. Rabitz, K. Shim, J. Chem. Phys. 111 (1998) 10640. [31] K. Shim, H. Rabitz, Phys. Rev. B 58 (1998) 12874. [32] H. Rabitz, K. Shim, J. Chem. Phys. 111 (1999) 10640. [33] J. Geremia, H. Rabitz, C. Rosenthal, J. Chem. Phys. 114 (2001) 9325. [34] J. Geremia, W. Zhu, H. Rabitz, J. Chem. Phys. 113 (2000) 10841.