Gaining Insights into Laser Pulse Shaping by Evolution Strategies Ofer M. Shir1 , Joost N. Kok1 , Thomas B¨ack1? and Marc J.J. Vrakking2 1
Natural Computing Group, Leiden University Niels Bohrweg 1, 2333 CA Leiden, The Netherlands {oshir, joost, baeck}@liacs.nl, http://natcomp.liacs.nl 2
Amolf-FOM, Institute for Atomic and Molecular Physics Kruislaan 407, 1098 SJ Amsterdam, The Netherlands
[email protected]
Abstract. We consider the numerical evolutionary optimization of dynamic molecular alignment by shaped femtosecond laser pulses. We study a simplified model of this quantum control problem, which allows the full physical investigation of the optimal solutions. By using specific variants of Derandomized Evolution Strategies, subject to parameterizations which are known to be superior for this problem, the numerical results reveal different conceptual physics structures for the different optimization procedures. These results are strong both from the algorithmic as well as from the physics perspectives. This shows that Natural Computing techniques can be used to derive new insights into Physics.
1
Introduction
To investigate and, more importantly, to control the motion of atoms or molecules by irradiating them with laser light, one has to provide laser pulses with durations on the same time scale as the motion of the particles. Recent technological development has made lasers with pulse lengths on the order of femtoseconds (1 fs=10−15 s) routinely available. Moreover, the time profile of these laser pulses can be shaped to a great extent. By applying a self-learning loop using an evolutionary mechanism, the interaction between the system under study and the laser field can be steered, and optimal pulse shapes for a given optimization target can be found. To this end, the role of the experimental feedback in the self-learning loop is played by a numerical simulation [1]. We plan to transfer, later on, these techniques to the laboratory for experimental optimization [2]. Control of molecular motion with shaped laser pulses is subject of intense current theoretical and experimental efforts. The success of femtosecond laser pulse shaping can considerably contribute even to the field of computation, as its application to Molecular Quantum Computing has been suggested (see, e.g., [3]). More specifically, molecular alignment is of considerable interest in this context ?
NuTech Solutions, Martin-Schmeisser-Weg 15, 44227 Dortmund, Germany.
because of its many practical consequences: a multitude of chemical and physical processes ranging from bimolecular reactions [4] to high harmonic generation [5] are influenced by the angular distribution of the molecular sample. Furthermore, in many fundamental molecular dissociation or ionization experiments the interpretation of the collected data becomes much easier when the molecules are known to be aligned with respect to a certain axis. Hence, techniques to generate molecular alignment are much needed. For a review, see [6]. A recent study presented a survey of modern evolutionary approaches to the problem, and showed that it payed off to use more elaborated optimization schemes, and in particular Derandomized Evolution Strategies (DES), for such a high-dimensional optimization task [7]. We rely on that study in our choice of two DES variants. In this study we focus in the simplified variant of the original problem, at zero temperature and with only a single rotational level in the initial distribution. The motivation for this simplification is to allow studying the physics nature of the optimal solutions, which would not have been possible for the general case, for reasons that will be explained in section 4.4 of this paper. The optimization procedure is subject here to two different parameterizations, which have been proposed for the laser problem [8]. The remainder of this paper is organized as follows. In section 2 we provide the reader with the details concerning the optimization routines in use. This is followed in section 3 with the introduction of the dynamic molecular alignment problem. In section 4 we present the experimental procedure, and section 5 outlines the conclusions which were drawn from this study.
2
Algorithms
Based on previous experience with this laser pulse shaping problem, and due to experimental results that showed that certain variants of Derandomized Evolution Strategies perform well with respect to other Evolutionary Algorithms on those problems [7], we restrict our study to these state-of-the-art algorithms. Our goal here is not to compare performance of algorithms, but rather to show that natural computing techniques can be used to derive new insights in physics. In this section we provide a short background of the specific variants in use.
2.1
Evolution Strategies
Evolution Strategies [9] are a canonical evolutionary algorithm for real-valued function optimization, due to their straightforward real-valued encoding, their specific variation operators, as well as to their high performance in this domain in comparison with other methods on benchmark problems. The higher the dimensionality of the search space, the more suitable a task becomes for an ES.
2.2
Derandomized Evolution Strategies
Mutative step-size control tends to work well for the adaptation of a global step-size, but tends to fail when it comes to the individual step size. This is due to several disruptive effects [10] as well as to the fact that the selection of the strategy parameters setting is indirect. The so-called derandomized mutative step-size control aims to tackle those disruptive effects. It is important to note that the different variants of derandomized-ES hold different numbers of strategy parameters to be adapted, and this is a factor in the speed of the optimization course: it is either a linear or quadratic order in terms of the dimensionality of the search problem n, and there seems to be a trade-off between the number of strategy parameters and the time needed for the adaptation/learning process of the step sizes. The (1, λ)-DR2 Algorithm The DR2 Algorithm [11] is considered to be the second generation of the derandomized Evolution Strategies. This variant uses a linear number in n of strategy parameters, and it aims to accumulate information about the correlation or anticorrelation of past mutation vectors in order to adapt the step size: xg+1 = xg + δ g δ gscal Z k
δ g+1
Z k = N (0, 1)
Z g = cZ sel + (1 − c) Z g−1 β g 1 |Z | −1+ = δ g · exp √ q c 5n n 2−c
g
βscal
|Z | g q δ g+1 + 0.35 scal = δ scal · +
c 2−c
(1) (2) (3)
(4)
n
where ξ scal = N (0, 1) , Z ∈ {−1, +1} , and β, βscal , b and ξ k are constants. The (µW , λ) Covariance Matrix Adaptation ES The (µW , λ)-CMA-ES algorithm [10] is known as the state-of-the-art among of the derandomized ES variants (could also be considered as DR4). It has been successful for treating correlations among object variables, where it applies principal component analysis (PCA) to the selected mutations during the evolution, also referred to as “the evolution path”, for the adaptation of the covariance matrix of the distribution. The concept of weighted recombination is introduced: applying intermediate multi-recombination on the best µ out of µ λ with given weights {wi }i=1 . The result is denoted with hxiW . Furthermore, (g) (g) pσ ∈ Rn is the so-called evolution path, pc ∈ Rn , sum of weighted differences
of points hxiW , C(g) ∈ Rn×n , the covariance matrix of the mutation distribution ¡ ¢T (C(g) = B(g) D(g) B(g) D(g) ): xg+1 = hxiW + σg Bg Dg z g+1 k
(5)
g+1
pg+1 = (1 − cc ) · pgc + cuc · cW Bg Dg hziW c ¡ g+1 ¢T pc Cg+1 = (1 − ccov ) · Cg + ccov · pg+1 c pg+1 σ
= (1 − cσ ) ·
σ g+1 = σ g · exp
pgσ
Ã
+
cuσ
1 · dσ
g
g+1 hziW
· cW B °! ° g+1 °pσ − χ ˆn ° χ ˆn
(6) (7) (8) (9)
where χ ˆn is the expected length of pσ . cc , ccov , cσ and dσ are p learning/adaptation µ u cc (2 − cc ), cW := rates, {w } are the recombination weights, and c := i i=1 c P µ p wi i=1 u √P 2 µ and cσ := cσ (2 − cσ ) are derived respectively. wi i=1
All weighting variables and learning rates were applied as suggested in the given citations, and particular in [10].
3
Dynamic Molecular Alignment
In this section we describe the optimization problem under investigation. The reader who wants to abstract from the physics details can view the problem as a single-criterion 80-dimensional optimization task, subject to maximization, with a punishment term for handling a physics constraint. 3.1
Quantum Control: Physics Background
The interaction of a generic linear molecule with a laser field is described within the framework introduced in [1]. We calculate the time evolution of a thermal ensemble of molecules quantum mechanically, by considering a single initial rotational level, characterized by the rotational quantum number Jinitial = 0 (and the projection of the angular momentum on the laser polarization axis Minitial , respectively). We take the molecule to be a rigid rotor, which allows a description of its wavefunction solely in terms of the rotational wave functions |JKM i (where K = 0 for a diatomic molecule). Two electronic states are taken into account, the ground state denoted by X and an off-resonant excited state denoted by A. Hence the wavefunction for a given M is expanded as ΨM (t) =
JX max
αXJM (t)ψXJM + αAJM (t)ψAJM
(10)
J=M
The time dependence of the molecular wave function is given by i
∂ΨM = (H0 + V )ΨM (t) ∂t
(11)
The Hamiltonian consists of a molecular part H0 and the interaction with the laser field, given by V = µ · E(t) cos(ωt)
(12)
The Eigenenergies of H0 are given by E(J) = hcBJ(J + 1)
(13)
where B is the rotational constant of the molecule. The laser field induces transitions between the rotational states which, in the off-resonant case, occur via subsequent Raman processes. The transitions between X and A were assumed to proceed via the selection rules ∆J = ±1, ∆M = 0. The envelope of the laser field (which completely determines the dynamics after the transition to the rotating frame has been performed) is described by Z E(t) = A(ω) exp(iφ(ω)t) exp(iωt) dω. (14) The spectral function A(ω) is taken to be a fixed Gaussian. The control function is the phase function φ(ω), which defines the phase at a set of n frequencies that are equally distributed across the spectrum of the pulse. These parameters are taken to be the decision parameters of the evolutionary search; the search space is therefore an n-dimensional hypercube spanning a length of 2π in each dimension. 3.2
Optimization
We consider the goal of optimizing the alignment of a sample of generic diatomic molecules undergoing irradiation by a shaped femtosecond laser. We have used ® the maximum cos2 (θ) that occurs under field free conditions after the laser pulse, where θ is the angle between the molecular and the laser polarization axis, as a measure of the alignment. The temperature of the ensemble was T = 0K and the rotational constant was chosen to be B = 5cm−1 . The peak Rabi frequency between the two electronic states X and A, that determines the interaction strength, was ΩXA = 1.6 · 1014 s−1 . Since we want to achieve a high degree of alignment with a peak intensity as low as possible, an additional constraint was introduced as a punishment for pulses that are too intense. We have used Z Ip = E 2 (t)Θ(E 2 (t) − Ithr ) dt (15)
(where Θ(x) is the Heaviside step function) for this purpose, so that the fitness function assigned to a pulse shape was ® F = maxE(t)=0 cos2 (θ) − βIp . (16)
By choosing β large enough, Ithr was shown [7] to effectively operate the evolutionary algorithms only on a subset of pulses whose maximum peak intensity approaches the threshold intensity from below. We have typically used β = 1; unless otherwise specified, Ithr was 0.36 · IF T L .
4
Experimental Procedure
Given the optimization routines which were introduced in section 2, we describe here the different parameterizations in use, present the numerical results, and finally apply an analytical tool for the investigation of the optimal solutions. 4.1
Parameterization
As introduced earlier, the phase function φ(ω) is the target function to be calibrated. Here, two parameterizations of φ(ω) are considered. The traditional approach was to interpolate φ(ω) at n frequencies {ωi }ni=1 ; the n values {φ(ωi )}ni=1 are the decision parameters to be optimized. In order to achieve a good trade-off between high resolution and optimization efficiency, the value of n = 80 turned to be a good compromise. We define this calibration of φ(ω), i.e. learning n = 80 function values and interpolating, as the so-called ’plain-parameterization’ optimization: φP (ω) = (φ(ω1 ), φ(ω2 ), ..., φ(ωn ))
(17)
An alternative parameterization for the phase function, which has been presented at [8], considers the Hermite polynomials, © ª¢ © ª dk ¡ k exp −x2 , k = 0, 1, ... Hk (x) = (−1) exp x2 k dx
(18)
as building blocks for the phase function, and aims to learn their coefficients in order to form the phase function: φH (ω) =
KX max k=0
ck · Hk (ω)
(19)
Note that the Hermite polynomials form a complete set of functions ©over the ª infinite interval −∞ < x < ∞ with respect to the weight function exp − 21 x2 . We define this as the ’Hermite-parameterization’ optimization. Kmax = 40 turned out to perform best. 4.2
Setup
Some technical details concerning the experimental setup and our modus operandi are outlined: – Every function evaluation has the duration of 7 seconds on a P4 2.6GHz. – Based on past experience, we choose the (1, 10) strategy for the DR2; (7, 15)CMA for n = 40 and (8, 17)-CMA for n = 80. – Each run is limited to 20, 000 function evaluations. – Implementation was done in Fortran for the numerical simulation (in order to stay close to the systems used by the physics researchers), and in MATLAB 7.0 for the optimization routines.
4.3
The Optimization: Numerical Results
Table 1 presents the mean values and the standard deviations of the cosinesquared alignment, obtained after 20 runs of 20, 000 function evaluations. As can be observed, the DR2 routine clearly outperforms the CMA in the plain parameterization, in consistency with previous results on the general problem initial Jmax = 7 (see, e.g., [7]). For the Hermite parameterization, however, the picture is different. The DR2 does not seem to deliver, and fails to obtain highquality solutions. The CMA does succeed in this task, with highly-satisfying results. Essentially, this suggests that the Hermite parameterization introduces strongly correlated decision parameters, whereas the plain parameterization can be tackled successfully by a strategy which does not consider the correlations between the decision parameters. It is nevertheless surprising to observe the low performance of the CMA on the latter. We may conclude that the Hermite pa-
maxE(t)=0 < cos2 (θ) > DR2 CMA Plain Parameterization 0.9559 ± 0.0071 (0.9622) 0.9413 ± 0.0058 (0.9508) Hermite Parameterization 0.9501 ± 0.0043 (0.9570) 0.9583 ± 0.0026 (0.9618) Table 1. Maximizing the cosine-squared alignment over 20 runs with 20, 000 function evaluations per run; mean and std ; the maximal value obtained is in brackets.
rameterization is slightly better, but not dramatically superior on this simplified variant, especially due to the fact that it requires the CMA routine in order to obtain optimal solutions, in comparison to the DR2 with the plain parameterization. 4.4
Investigation of Optimal Solutions
An optimal solution is represented by its phase function φ (ω), and the electric field respectively, but one can also examine the revival structure. Only due max to our simplified variant, i.e. Jinitial = 0 at initialization, it is possi3 ble to study the population of the rotational levels as a function of time. Otherwise, in the general case, all levels are initially populated, and a thermal averaging is applied. Explicitly, the wavefunction can be expressed as a superposition of those levels, X (t) Ej t (20) ψ= aj · |ji · e−i ~ j
the expectation of the cosine-squared alignment (the objective function) is calcu(t) lated directly from these complex amplitudes aj , whereas the population of the 3
The careful reader should note that ’population’ is used here exclusively in the context of quantum mechanics, e.g., populating quantum levels.
¯ ¯ ¯ (t) ¯2 rotational levels is ¯aj ¯ . This population of rotational levels can be analyzed in a fairly simple technique, known as the Sliding Window Fourier Transform (SWFT), which provides us with a powerful visual tool. Given the revival structure of an optimal solution, a sliding time window is Fourier transformed, to produce the frequency picture through the alignment process. This windowing creates a transformation which is localized in time. Due to the quantization of the rotational levels, only certain frequencies (or energy levels, respectively) are expected to appear.
Analysis Results We applied the SWFT routine to the optimal solutions which were found in the various runs. Figures 3, 4, 5 and 6 visualize the typical population process of the rotational levels for four typical solutions of the different optimization procedures (2 parameterizations times 2 DES variants). The quantum energy levels are indeed observed as expected from theory. The results reveal two different conceptual physical structures, which correspond to optimal and sub-optimal solutions in terms of the success-rate, i.e., the cosine-squared alignment. The Plain-DR2 as well as the Hermite-CMA procedures obtain the best solutions, which share the same structure - they are characterized by the dominant population of the 4th quantum energy level in the SWFT picture. On the other hand, the Plain-CMA and Hermite-DR2 procedures obtain inferior solutions, which are characterized by a gradually increasing population of the energy levels. The original revival structures for two optimal solutions, representing the two conceptual structures, are given in Figures 1 and 2. The optimal family of solutions (Fig. 1) presents a dramatic revival structure, with a typical strong pulse in the train which lies on the boundary of the punished regime (I ≈ 0.36). This pulse seems to be essential in giving the molecules the right ’kick’, and most likely responsible for the dominant population of the 4th quantum energy level in the SWFT picture. The sub-optimal family of solutions (Fig. 2) yields a revival structure with a smooth exponential envelope, and thus has a gradual building-up of the quantum energy levels in the SWFT picture, respectively. It typically contains a train of medium pulses and lacks a dominant one. We would like to emphasize the fact that we obtained the same family of optimal solutions, with the same structure, from two different optimization approaches - the one learning a vector of successful variations of interpolated points of the target function, whereas the other learning a covariance matrix of coefficients of Hermite polynomials that span the target function. Based on our experience with the alignment problem, and due to the new results, we claim that this might suggest that the regime of the global optimum has been reached. This could have some strong physics consequences, which needs to be investigated carefully by the physicists, but as far as the algorithmic perspective is concerned, this seems to be the case.
Alignment and Revival Structure of two obtained solutions (Fig. 1 and 2). Thin red line: alignment; thick black line: intensity of the laser pulse. 5
10 1
−10 1
−8
−6
−4
−2
0
2
4
6
8
10 1
0.9
0.9
0.8
0.8
0.8
0.8
0.7
0.7
0.7
0.7
0.6
0.6
0.6
0.6
0.5
0.5
0.5
0.5
0.4
0.4
0.4
0.4
0.3
0.3
0.3
0.3
0.2
0.2
0.2
0.2
0.1
0.1
0.1
0.1
0 −5
0
0.9
0 10
5
0 −10
−8
−6
−4
−2
t [ps]
0
2
4
6
Field Intensity
0
Field Intensity
−5 1 0.9
0 10
8
t [ps]
Fig. 1. Optimal solution type: no governing structure for the revival structure is observed; < cos2 (θ) >= 0.9622.
Fig. 2. Suboptimal solution type: a smooth exponential envelope on top of the revival structure; < cos2 (θ) >= 0.9505.
Each of the SWFT figures (Fig. 3 - 6) represents a Fourier transform applied to the revival structures of the optimal solutions (the thin-red alignment curves of Fig. 1-2). The values are log-scaled, and represent how high the rotational levels of the molecules are populated as a function of time. Thus, an exponential envelope (Fig. 2) is represented by a gradual building-up of frequencies (Fig. 4). Note that the quality of the laser pulse cannot be measured in those plots. 2
2
Plain−DR2: =0.9622
−11
1
x 10
Plain−CMA: =0.94695
−11
1
x 10
−6
−6
0.5
0.5 −8
−8
−10
−12 −0.5
−10
0
Time
Time
0
−12 −0.5
−14
−14
−16
−16
−1
−1 −18
−1.5
0
5
10
15
Energy [B
rot
20
Units]
25
30
−18 −1.5
35
0
5
20
Units]
25
30
35
Fig. 4. All five first rotational levels are populated gradually after the interaction.
2
−11
x 10
15
Energy [B
rot
Fig. 3. 4th rotational level is mostly populated after the interaction.
1
10
2
Hermite−DR2: =0.94942
−11
1
x 10
Hermite−CMA: =0.96185
−6
−6
0.5
0.5 −8
−8
−10
−12 −0.5
−14
−10
0
Time
Time
0
−12 −0.5
−14
−16 −1
−16 −1
−18 −1.5
0
5
10
15
Energy [B
rot
20
Units]
25
30
35
Fig. 5. The four first rotational levels are populated gradually after the interaction.
−18 −1.5
0
5
10
15
Energy [B
rot
20
Units]
25
30
35
Fig. 6. 4th rotational level is mostly populated after the interaction.
5
Discussion and Outlook
We have applied derandomized Evolution Strategies, subject to two parameterizations, to the numerical optimization of dynamic molecular alignment. Two different approaches obtained optimal solutions with numerical results in the same rank of quality. Furthermore, the investigation of the optimal solutions, which was performed here for the first time, revealed two typical physics structures, through the time-evolution of the populated quantum energy levels. The first structure typically contains a strong pulse, which is followed by the population of a specific quantum rotational level for the molecules, whereas the second structure is a train of medium pulses with an outcome of a gradual population of the rotational levels. This confirms the multi-modality of the search space, and provides us with strong physics intuition with respect to optimal versus sub-optimal solutions, e.g., the shape of optimal pulses as well as the optimal population of rotational levels. Our new observation suggests that the regime of optimal solutions has been found.
Acknowledgments This work is part of the research programme of the ’Stichting voor Fundamenteel Onderzoek de Materie (FOM)’, which is financially supported by the ’Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)’.
References 1. Rosca-Pruna, F., Vrakking, M.: Revival structures in picosecond laser-induced alignment of i2 molecules. Journal of Chemical Physics 116(15) (2002) 6579–6588 2. Zamith, S. Eur. Phys. J. D 12(255) (2000) 3. Tesch, C.M., Kurtz, L., de Vivie-Riedle, R. Chem. Phys. Lett. 343 (2001) 4. Friedrich, B., Herschbach, D. Phys. Chem. Chem. Phys. 2(419) (2000) 5. Hay, N. Phys. Rev. A 65(053805) (2000) 6. Stapelfeldt, H. Rev. Mod. Phys. 75(543) (2003) 7. Shir, O.M., Siedschlag, C., B¨ ack, T., Vrakking, M.J.: Evolutionary algorithms in the optimization of dynamic molecular alignment. In: 2006 IEEE World Congress on Computational Intelligence, IEEE Computational Intelligence Society (2006) 9817–9824 8. Shir, O.M., Siedschlag, C., B¨ ack, T., Vrakking, M.J.: The complete-basis-functions parameterization in es and its application to laser pulse shaping. In: GECCO ’06: Proceedings of the 8th annual conference on Genetic and evolutionary computation, New York, NY, USA, ACM Press (2006) 1769–1776 9. Beyer, H.G., Schwefel, H.P.: Evolution strategies a comprehensive introduction. Natural Computing: an international journal 1(1) (2002) 3–52 10. Hansen, N., Ostermeier, A.: Completely derandomized self-adaptation in evolution strategies. Evolutionary Computation 9(2) (2001) 159–195 11. Ostermeier, A., Gawelczyk, A., Hansen, N.: Step-size adaptation based on non-local use of selection information. In: Parallel Problem Solving from Nature (PPSN3). (1994)
