Freezing molecular vibrational energy flow with coherent control

Journal of Molecular Structure (Theochem) 589–590 (2002) 447–457 www.elsevier.com/locate/theochem

Freezing molecular vibrational energy flow with coherent control R.M. Bigwooda,1, M. Gruebelea,b,c,* a

Department of Chemistry, University of Illinois, Urbana, IL 61801, USA b Department of Physics, University of Illinois, Urbana, IL 61801,USA c Center for Biophysics and Computational Biology, University of Illinois, Urbana, IL 61801,USA Received 7 May 2001; revised 8 May 2002; accepted 8 May 2002

Abstract In a previous report, we examined the possibility of ‘freezing’ intramolecular vibrational energy redistribution (IVR) in organic molecules by using optimal coherent control. Here we describe the new methodology developed to achieve stabilization of initially excited nonstationary states. Our approach combines an approximate but full-dimensional molecular vibrational Hamiltonian, a frequency-domain wavelet representation of the electric field, a fast symplectic propagator to compute the IVR decay, and simulated annealing optimization of the electric field parameters. We find that the complexity of the vibrational wavepacket increases sufficiently slowly with time, so that with available pulse shapers, vibrational dephasing can be ‘frozen’ for 1 – 2 orders of magnitude in time beyond the usual IVR decay time. Slowing the IVR clock into the multi-picosecond regime may allow the natural selective reactivity (via Franck– Condon factors) of the initially excited nonstationary vibrational states to emerge. q 2002 Elsevier Science B.V. All rights reserved. Keywords: Franck–Condon factors; Thiophosgene; Wavelets; Strong field control; Symplectic propagator

1. Introduction Photoselective chemistry has moved tantalizingly closer to the goal of chemical reaction control. Two different approaches to this goal exist. Franck – Condon control relies on the shape of the molecular potential energy surfaces to guide the reaction [1 – 3]. Coherent control in addition makes use of interfering quantum phases [4 – 6]. Experimental examples range from interference between just two paths [7], to manipulation of molecular wavepackets using hundreds of phases [8,9]. Franck– Condon and coherent approaches can be combined [10], and optimal control * Corresponding author. Fax: þ 1-217-244-31-86. E-mail address: [email protected] (M. Gruebele). 1 Current address: Intel Portland.

theory provides the framework for finding optimal optical control fields [11]. Intramolecular vibrational energy redistribution (IVR) poses difficulties for both approaches to selectivity. Highly vibrationally excited wavepackets dephase rapidly [12 – 14], exploring parts of the potential surface undesirable from the point of view of reactivity, and losing the simple nodal structure of the initially prepared state so useful for coherent control. The IVR clock can be beaten by speeding up the reaction [3], or by slowing down IVR. Recent work has shown that IVR requires much longer to reach the statistical limit than predicted by exponential decay models [14 – 16]. Even at high energy, only a limited set of vibrational resonances mediates diffusion of the vibrational wavepacket in the 3N 2 6-dimensional vibrational state space. The

0166-1280/02/$ - see front matter q 2002 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 6 - 1 2 8 0 ( 0 2 ) 0 0 3 0 3 - 2

448

R.M. Bigwood, M. Gruebele / Journal of Molecular Structure (Theochem) 589–590 (2002) 447–457

Fig. 1. Number of control parameters for rephasing an interior state of SCCl2 plotted as a function of time. D(t ) was calculated using the basis set described in Ref. [19]. At short times, the number of control parameters increases exponentially with time, but then it continues to increase as a power law, dramatically lengthening the time over which rephasing is possible with a given number of control channels.

actual dimension of the vibrational wavepacket lies in the range d < 1 – 4 ! 3N 2 6; [17] and the survival probability PðtÞ ¼ lkCð0ÞlCðtÞll2 of an initially prepared bright state decays as a slow powerlaw [18,19]. Slow powerlaw dephasing should make it easier to control IVR. In one scenario, IVR transports energy to the reaction coordinate, which cannot be excited directly [20]; coherent control could enhance ‘funneling’ of the energy towards the reaction coordinate. In another scenario, the initially prepared state has the desired reactivity, but dephases rapidly; coherent control could ‘freeze’ the initially prepared state to preserve the desired reactivity [10]. The feasibility of freezing IVR by coherent control is related to the wavepacket dispersion DðtÞ < PðtÞ21 ; which measures the minimal number of states required to represent the vibrational wavepacket up to time t (Fig. 1) [19]. Because each state requires a phase and amplitude control parameter in the electric field [14], Number of control parametersðtÞ < 2DðtÞ:

ð1Þ

Fig. 1 shows an upper bound to D, for IVR of the molecule SCCl2 (E < 8,000 cm21). The initial exponential decay of D rapidly switches to a power law, lengthening the control time available with 103

Fig. 2. Energy level scheme for SCCl2 based on the data in Ref. [22]. The states lgl (impure ‘ground’ state at finite temperature),l0l (pure initial state selected with high-resolution laser), lbl (fs-excited bright state consisting of basis states lil or eigenstates lnl depending on the representation chosen) required for the control process are shown. Only the rCS coordinate of the six-dimensional potential energy surface of SCCl2 is shown because the initial state is mainly excited in this coordinate; subsequent dephasing occurs among all vibrational coordinates.

control channels by at least an order of magnitude compared to an exponential decay. We previously published results indicating that IVR can indeed be ‘frozen’ by coherent control [10]. Here, we focus on the new computational methodology we developed for our IVR control calculations, using a vibrational Hamiltonian for thiophosgene (SCCl2) as an example. The control field E(t ) is treated classically, and obtained from a spectral wavelet representation. This allows a realistic simulation of experimental amplitude/phase masks. We use a new version of the SUR symplectic propagator [21] to compute IVR dynamics, which are then optimized by simulated annealing to freeze the bright state for as long as possible. The target functional is chosen to freeze a state over an extended period of time, in contrast to the usual target functionals for reaching a target at a certain time or as t ! 1: We conclude that the vibrational Hamiltonian poses a much less severe limit to controlling reactivity than previously thought.

R.M. Bigwood, M. Gruebele / Journal of Molecular Structure (Theochem) 589–590 (2002) 447–457

449

Table 1 Vibrational and electric field parameters for the simulations in Section 3 Mode i

1 (CS str.)

2 (CCl str.)

3 (bend)

4 (umbr.)

5 (CCl str.)

6 (wag)

vi (cm21) xii (cm21) V ¼ 31.5 cm21 (28 cm21 in Fig. 7) Iint ¼ 1 J cm22 (Gaussian initial pulse)

1138 4 a ¼ 0:15 (0.17 in Fig. 7) m0b ¼ 2 Debye

500 0.5

289 1.5

471 0.2

818 1.5

301 0.6

2. Quantum-classical Hamiltonian and dynamics In this section we discuss the representation of the Hamiltonian and time propagator in detail. In the quantum-classical approximation, the Hamiltonian for a molecule interacting with a laser pulse can be written as HðtÞ ¼ Hmol þ Hint ðtÞ:

ð2Þ

Hmol is the time-independent molecular Hamiltonian. Hint is the time-dependent molecule-field interaction Hamiltonian with matrix elements of the type ðHint Þb0 ¼ kbl 2 m·EðtÞl0l:

ð3Þ

Here, l0l is from a manifold of initial states, and lbl is from a manifold of ‘bright’ states which carry the oscillator strength. With energy bandwidths typically accessible by a control pulse on the order of a few 100 cm21, only one or a small number of initial states l0l and bright states lbl is involved in IVR control, even though the density of eigenstates near lbl can be large. Fig. 2 summarizes the situation for the control of the n1 ¼ 8 state (CS stretching overtone) of SCCl2 ˜ electronic from a low-lying vibrational level in the A state. We will consider this example in detail here. For a 300 cm21 FWHM bandwidth of the control pulse, ˜ electronic only one vibrational state l0l on the A ˜ surface, and one bright state lbl ¼ l800000l on the X electronic surface need to be considered to a good approximation, although several hundred dark ˜ electronic surface mix vibrational states lIl on the X with lbl, causing the dephasing that is to be controlled. We pick SCCl2 because it combines a smalldimensional system tractable by full vibrational quantum dynamics with a high density of states that leads to a significant IVR decay to be controlled. The n1 ¼ 8 state has been found experimentally and computationally to lie above the IVR threshold [19,

22]. In our calculations, we use an effective vibrational Hamiltonian ð0Þ ð1Þ þ Hmol ¼ Hmol ¼ Hmol

6 X

  vi ni þ 12

i¼1

 þ xii ni þ

1 2

2

ð1Þ þHmol :

ð4Þ

ð0Þ Hmol comprises the harmonic vibrational energy and off-resonant anharmonic corrections. Its basis functions are anharmonic vibrational feature states. These ð1Þ states are mixed by the resonant couplings of Hmol : As a result, the spectrum of a bright feature state lbl breaks up into an IVR multiplet containing many ð1Þ eigenstates. We chose Hmol based on our averaged scaling model for the vibrational Hamiltonian, which successfully predicts IVR properties such as lifetimes and dilution factors [16,23,24]. Couplings between vibrational features in this model P are given by ð1Þ ðHmol Þn;n0 ¼ Va32n ; where n ¼ i¼1 – 6 lni 2 n0i l (2 should be added to n if n , 3). The model thus includes cubic and higher-order couplings. Optimal values for the vi, xii, V and a are based on Refs. [19, 22,23] and summarized in Table 1. The major qualitative shortcoming of Eq. (4) is its neglect of rotational structure, which affects both molecular energy levels and the m·E alignment. As the control pulse cycles population between l0l and the vibrational eigenstates overlapped with lbl, rotational ladder climbing and Coriolis couplings introduce additional dephasing, but these effects are small for heavy molecules or low J values. For now, our main question is whether the vibrational Hamiltonian precludes or allows control. Once the molecular Hamiltonian to be controlled has been defined, the interaction of the molecule with the control field given by Eq. (3) is in principle straightforward. In practice, keeping the laser pulse in the realm of what can be produced in the laboratory is

450

R.M. Bigwood, M. Gruebele / Journal of Molecular Structure (Theochem) 589–590 (2002) 447–457

Fig. 3. Summary of the control calculation. Input data are shown on the left. The iterative control cycle using a Daubechies wavelet representation for the field is shown on the right. The relationship among various quantities and algorithmic procedures is indicated by arrows.

critical. The maximum intensity is limited by ionization processes (unless they are desired). Currently available fast lasers and pulse shapers are limited in bandwidth, and in the number of independently variable phase and amplitude control channels [25,26]: starting with a near-Fourier transform limited pulse, amplitudes Ai and phases wi can be controlled in 100 –1000 frequency channels vi. In addition, it is well known from experiment and simulation that the time-dependence of control pulses is often complicated [8,9] Writing the electric field as EðtÞ ¼ 1ðtÞ cosðv0 t þ w0 Þ;

ð5Þ

in terms of a carrier wave at v0, the envelope function 1(t ) can be highly structured, even when the frequency domain representation is fairly smooth. A frequency domain representation that incorporates the experimental limitations is therefore desirable. To represent pulse shapers, we model the electric field in the frequency domain and then transform to the time domain for wavepacket propagation. A number of bases exist in which the electric field 1(v ) can be represented, including harmonic functions, Chebyshev polynomials, and wavelets. Daubechies-4 wavelets are ideal for this application (Fig. 3): the pyramidal algorithm to compute 1(v ) from wavelets is fast; wavelet locality allows a nearly physical representation of the phase and amplitude channels of experimental pulse shapers; wavelets’ hierarchical resolution allows optimization to proceed at both low and high spectral resolution [27]. Detailed

discussions of wavelets are available elsewhere [28]. Briefly, one step in the 4 coefficient wavelet transform can be represented in matrix form as 2

c0

6 6c 6 3 6 6 6 6 6 6 6 6 6 6 DAU4 ¼ 6 6 6 6 6 6 6 6 6 6 6 6 6 c2 4 c1

c1

c2

c3

2c2

c1

2c0

3

c0

c1

c2

c3

c3

2c2

c1

2c0

·

· ·

· c0

c1

c2

c3

2c2

c1

c3

c0

2c0

c3

7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 c3 7 7 7 2c0 7 7 7 7 c1 7 5 2c2

ð6Þ When applied to a vector representing the electric field as a function of frequency, the odd-numbered rows of DAU4 smooth the vector, while the evennumbered rows are sensitive to detail in the vector. To perform a wavelet transform, the matrix in Eq. (6) is applied to the vector as follows: the vector is multiplied by Eq. (6), and all of the smoothed rows are permuted to the first half of the vector. The process is repeated for the halved smooth data until only two smoothed components remain. The final hierarchical vector assembled from the 2 final smoothed and 2 þ 4 þ 8 þ … ‘detail’ components then contains the wavelet coefficients. These wavelet coefficients describe the original vector hierarchically from lowest

R.M. Bigwood, M. Gruebele / Journal of Molecular Structure (Theochem) 589–590 (2002) 447–457

Fig. 4. Initial form of the time-dependent Hamiltonian matrix, and final prediagonalized forms useful for iterative control calculations are shown. All results in Figs. 5 and 6 were computed using the fully prediagonalized basis where {lbl,lil} ! {lnl} and state l0l is connected to all states lnl via time-dependent matrix elements in a Golden Rule form.

to highest resolution. The inverse transformation is performed by applying the transpose of the matrix in Eq. (6 ) and reversing the sorting. Because of cusps in the Daubechies-4 wavelets, higher order wavelet transforms must be used if the function needs to be differentiable. For our IVR control application, differentiability is not required, and we prefer the higher locality of the Daubechies-4 wavelets, which mimic individual channels or groups of channels in spectral masks used experimentally. The electric field is modeled in the frequency domain by using two wavelet vectors, one for the phases wi and one for the amplitudes Ai. By adjusting the wavelet coefficients, the phase and amplitude are simultaneously optimized at all resolution levels. An inverse wavelet transform recovers the phases and amplitudes in each frequency channel vi. The real and imaginary parts of the electric field in each frequency channel, 1ðvi Þ ¼ 1ðrÞ ðvi Þ þ i1ðiÞ ðvi Þ are retrieved using the usual relationships tan wi ¼ EðrÞ ðvi Þ=EðiÞ ðvi Þ

ð7Þ

A2i ¼ EðrÞ ðvi Þ2 þ EðiÞ ðvi Þ2 :

ð8Þ

451

Finally, the electric field vector is Fourier transformed to the time domain EðtÞ ¼ EðrÞ ðtÞ þ iEðiÞ ðtÞ to time propagate the initial bright basis state. We now turn to the structure of the matrix representation into which Eq. (2) is cast for time propagation (Fig. 4). The molecular basis set contains three important classes of states: {l0l}, from which the bright state undergoing IVR is optically accessed, {lb}, the bright state which undergoes IVR, and {lil the manifold of dark feature states coupled to {lbl}. In the present example, {l0l and {lbll} each contain only one state, while lil contains m 2 2 states, leading to an m £ m matrix whose structure is shown in Fig. 4. The bulk of the computational effort comes from a large number of time propagations needed to test each trial pulse shape. In practice the {l0l,lbl,lil} representation of the Hamiltonian is inefficient because of the large number of off-diagonal coupling matrix elements in the IVR part of the Hamiltonian. When a large number of control channels is included in the optimization, requiring many P(t ) evaluations, a significant speed advantage is gained by prediagonalizing the molecular part of the Hamiltonian to form a ‘Golden Rule’ matrix which has only one state coupled to all others. Two useful representations are shown in Fig. 4. The {lil} manifold can be prediagonalized to yield uncoupled dark states {ln0 l} coupled to lbl, which in turn is coupled to l0l, leaving only m 2 1 off-diagonal matrix elements, one of which is time-dependent; or {lbl,lil} can be prediagonalized to yield the IVR eigenstates, all of which are coupled to l0l by matrix elements knlm·EðtÞl0l ¼ knlblkblm·EðtÞl0l:

ð9Þ

This representation will be used for all calculations presented here. It looks like a Golden Rule matrix, with state l0l coupled by m 2 1 time-dependent matrix elements to a prediagonalized manifold of m 2 1 states at energies En. The matrix-fluctuation-dissipation (MFD) theorem was used to compute the eigenvalues En and overlaps knlbl for the large matrices [29,30]. In the prediagonalized eigenstate representation, more than 10-fold increases in propagation speed are achieved. Finally, we consider the propagation algorithm itself. Symplectic integrators, such as the modified SUR algorithm described below, require a time step inversely proportional to the energy bandwidth, and will be slow without some manipulation of H(t ). In

452

R.M. Bigwood, M. Gruebele / Journal of Molecular Structure (Theochem) 589–590 (2002) 447–457

the present case, the energy bandwidth is approximately the greater of v0b or Dvib (the bandwidth of the states undergoing IVR). The former can be very large for visible or UV transitions, but is easily removed by making v0b the carrier frequency, then subtracting the carrier frequency from the excited states and making the rotating wave approximation. EðtÞ ¼ 1ðtÞðeivt þ e2ivt Þ in Eq. (5) is thus replaced by EðtÞ ¼ 1ðtÞð1 þ e22ivt Þ;

ð10Þ

in which we neglect the term oscillating at twice the transition frequency. The approximation is valid if the laser is nearly resonant with the transition, and the field strength is such that m·E ! v þ v0b : Because v 21 generally lies in the 0.5 – 2 fs range, rather large electric fields are still allowed. However, the A2 vector potential term in the Hamiltonian is neglected, so none of the calculations are in the strong field limit. To reduce the bandwidth problem caused by Dvib as much as possible, we further chose Eb ¼ 0 as the origin for all energies in the Hamiltonian matrix. The biggest problem for many propagators is that the shaped E(t ) in the rotating wave approximation is a complex time-varying function. Large-step propagators (e.g. Chebyshev) [31] cannot be used to advantage when H depends on time. Higher order integrators such as S3 [32] or the more general Bulirsch– Stoer algorithm [27] gain their accuracy by exploiting subtle error cancellations. These algorithms are destabilized by matrix elements that vary at a rate comparable to the time step. The SUR propagator is ideally suited for this type of problem [21]. It is the simplest of the symplectic integrators [33,34], numerically equivalent to the leapfrog algorithm, but requiring less storage and operating twice as fast. Symplectic propagators exploit the fact that the Schro¨dinger equation can be written as [33]

C_ r C_ i

! ¼

0

H^

^ 2H

0

!

Cr

! ð11Þ

Ci

in terms of the real and imaginary parts of the wavefunction. Using the basis set expansion ltl ¼ P ðrÞ ðiÞ j ðcj þ icj Þljl the SUR propagator for a real

Hamiltonian matrix can be written as [21] X ðrÞ ðiÞ ðrÞ Hij cj ðtÞ; cðrÞ i ðt þ DtÞ ¼ ci ðtÞ þ Dt j ðiÞ cðiÞ i ðt þ DtÞ ¼ ci ðtÞ 2 Dt

X

ð12Þ HijðrÞ cðrÞ j ðt þ DtÞ:

j

The inclusion of complex matrix elements requires that both the real and imaginary components of the vectors c appear in each step of Eq. (13). For coherent control calculations, we have implemented a more general version of SUR capable of dealing with timevarying complex matrix elements. The implementation can be written as o X n ðrÞ ðrÞ ðiÞ Hij þ kðrÞ cðrÞ i ðt þ DtÞ ¼ ci ðtÞ þ Dt ij ðtÞ cj ðtÞ j

þ Dt

X

ðrÞ kðiÞ ij ðtÞcj ðtÞ;

j ðiÞ cðiÞ i ðt þ DtÞ ¼ ci ðtÞ 2 Dt

Xn

o ðrÞ HijðrÞ þ kðrÞ ij ðtÞ cj ðt þ DtÞ

j

þ Dt

X

ðiÞ kðiÞ ij ðtÞcj ðtÞ;

j

ð13Þ in analogy to our original algorithm. Here kij ðtÞ ¼ mij 1ðtÞ is the product of the transition dipole matrix element and the field envelope 1(t ), and the frequency v0b must be subtracted from the states lbl and {lil} to dress them. The algorithm converges to the exact quantum dynamics as Dt ! 0; and provides excellent performance as long as the time step is chosen small enough to sample the substructure of 1(t ). Because of its one-step nature, the algorithm is also suited for variable Dt grids. This allows fine grids to be used while the electric field is turned on, and larger grids to be used while free propagation occurs. Other propagators may be competitive for IVR computations in the absence of a laser pulse, but SUR is uniquely suited for use in calculations that include a time-dependent molecule – radiation interaction together with IVR.

3. IVR control We now turn to the coherent control of IVR in the

R.M. Bigwood, M. Gruebele / Journal of Molecular Structure (Theochem) 589–590 (2002) 447–457

n1 ¼ 8 state of thiophosgene. The computational ingredients from Section 2 are summarized in Fig. 3. We will use it as a guide through the control calculation. The initial state for the field is a 300 cm21 FWHM Fourier-transform limited pulse. The integrated pulse energy is given in Table 1. For the remainder of this section, the maximum amplitude of the initial Gaussian pulse is normalized to 1 and all other amplitudes are referenced to it. The initial molecular state is the lbl ¼ 800000l anharmonic basis state of SCCl2. All other anharmonic basis states {lil} in a ^ 400 cm21 around lbl were generated. To construct the Hamiltonian matrix, states from the ensemble {lil} were selected by a tier criterion Lij ¼ ½1 þ ðDEij =Hij Þ2 21 . 0:003; where DEij is the energy gap between a state lil in tier n and a state ljl in tier n þ 1; and Hij is the corresponding coupling matrix element in the IVR Hamiltonian [24]. The bright state is in tier 0, and the matrix was computed up to 4 tiers deep, and then prediagonalized by MFD to yield the energies lnl and overlap matrix elements kblnl [30]. The usual target in optimal control is a wavepacket at a specific time t, or as t ! 1 [6]. In the present case, we are not trying to guide the wavepacket elsewhere, but rather to stabilize it in an initially prepared state which is not an eigenstate of Hmol. By adjusting 1(t ), the goal is to make lbl a near-stable state of H(t ). Ideally, one would like the time-evolving IVR wavepacket ltl to remain perfectly localized at the initially prepared bright state lbl, with the survival probability PðtÞ ¼ lkbltll2 unity up to a time tmax longer than the reaction time scale. Merely requiring that lkbltmax ll2 be unity does not take into account the behavior of P(t ) at intermediate times. The integrated survival probability ðtmax W¼ PðtÞdt ð14Þ 0

provides a much more robust measure of the effectiveness of the control for a given electric field. A few trivial optimization outcomes must be avoided. The first is over-weighting pulse shapes that maximize Eq. (14) by efficiently pumping a large population into the target wave-function at short times, without slowing the decay rate. While such pulses may be of interest in other applications, they are of no help in creating long-lived wave-packets for

453

selective chemistry experiments. On the contrary, a small population in lbl is perfectly acceptable as long as the remainder stays in the ‘unreactive’ state l0l. Normalizing W by ð1 2 c0 Þ; the population remaining in state l0l at long times, addresses this issue. Another uninteresting result is obtained when the laser bandwidth is reduced to the point of pumping a single, or a few closely spaced eigenstates. Eigenstates of Hmol are quasi-stationary at low fields so their P(t ) does not decay significantly, but this is not a very useful type of ‘freezing.’ Calculating P(t ) by projecting onto the eigenstate spectrum of the bright-state lbl, in conjunction with the above normalization, eliminates this possibility. The target functional becomes  n   2 X max   ðrÞ ðiÞ kblnl cn ðtÞ þ cn ðtÞ  ðtmax   T½1ðtÞ ¼ dt n¼1  2  2 :    ðiÞ  tmin 1 2 cðrÞ 0 ðtÞ 2c0 ðtÞ

ð15Þ

Excessive pulse intensities are avoided by letting the amplitude mask act only as an attenuator, just as in real experimental situations. Once the adjustable variables (wavelet coefficients used to represent the electric field), and the target function are defined, the optimization problem is completely defined. Most fast function optimizers accept a new set of parameters only if they produce result that is closer to the desired value. This makes them susceptible to falling into a local minimum. Simulated annealing (SA) and the genetic algorithm (GA) can be used to avoid this problem by allowing a move to a point that is farther away from the desired value [27]. Although GA provides satisfactory results and has been extensively used in coherent control calculations [8,9,35,36], here we test simulated annealing. In analogy to thermodynamic annealing, SA requires a definition of an ‘energy’ parameter, which can simply be the negative of the target functional, and a temperature. All parameter sets that result in a lower energy are automatically accepted. Those that result in an increase in energy are accepted if a Boltzmann weighting factor , e2E=kT exceeds a cutoff value. There are of course trade-offs between the rate of convergence and the ability to move out of local minima. Setting a high initial temperature will allow a

454

R.M. Bigwood, M. Gruebele / Journal of Molecular Structure (Theochem) 589–590 (2002) 447–457

Fig. 5. Top row: initial and optimal control field. Middle row: initial and optimal (slowest) IVR decay. Bottom row: real and imaginay parts of the control field, showing an inversion symmetry of the four-pulse sequence with respect to the pulse center.

large number of less optimal solutions to be selected, decreasing the probability of getting stuck in a local minimum. On the other hand, setting the temperature too high for the ‘roughness’ of the target functional merely wastes time. The algorithm is generally started with a fairly high temperature (here equal to the value of W obtained with a transform limited Gaussian pulse), and the temperature is gradually lowered as the optimization progresses. In our simulations, the temperature was allowed to decrease at a constant rate. The results of the iterative SA optimization runs are shown in Figs. 5 and 6. Sixty four field amplitudes and phases were adjusted, half the number available in the smallest commercial pulse shapers. Fig. 5 shows the initial and final time domain pulse envelopes at the top, and the initial and final IVR decays in the middle. About 70 iterations were required to achieve the final slowed-down IVR decay. When pumped with a Fourier transform-limited Gaussian pulse, the ‘uncontrolled’ IVR proceeds with a 100 fs initial decay, followed by a series of quantum beats that slowly die

off with a much slower rate.. The shaped pulse preserves at least 20% of the initial state for about 10 ps, and at the same time increases the maximal fraction of population transferred from the initial state by a factor of three. The lifetime of the l800000l feature state is thus extended by a factor of 100. In the time domain, the initially Gaussian electric field is replaced by a four-pulse sequence, whose real and imaginary parts are approximately antisymmetric with respect to the pulse center. Finally, note that some degree of control persists for a duration comparable to the shaped pulse width after the pulse sequence has been completed. Fig. 6 shows the optimized 64 amplitude and phase components in the frequency domain. The amplitude is structured, with approximately eight maxima in a 500 cm21 interval. The phase is considerably less structured, particularly if one takes into account that the strongest outliers w(vi) occur at small field amplitudes. In Fig. 6, the phase points are amplitude-coded by filled black circles so it is immediately

R.M. Bigwood, M. Gruebele / Journal of Molecular Structure (Theochem) 589–590 (2002) 447–457

455

Fig. 6. Amplitude and phase spectra of the optimal control field. The phase spectrum is coded with dots whose size is proportional to the field amplitude at the same frequency, to highlight the most important phase components.

evident which phase points are most important in determining the structure of the electric field. The most important phase variations are , 1 rad. Robustness is an important issue for any coherent control calculation. If the electric field is very sensitive to the underlying Hamiltonian, calculations can demonstrate feasibility of control, but cannot provide a useful initial guess for experimental work. IVR control turns out to be very robust, as can be seen by considering sensitivity either to variations in the control field, or to variations in the underlying Hamiltonian. Insensitivity to small changes in the field is evidenced by a rapid

approach to convergence as the coarse structure in Fig. 6 appears, and a slower rate of convergence as the fine structure builds in. Robustness is also evidenced by the fact that a reduction in channels still produce detectable ‘freezing’.

4. Discussion The results of Section 3 show that the vibrational Hamiltonian by itself allows coherent control of IVR. Polyatomic vibrational wavepackets can be frozen to

456

R.M. Bigwood, M. Gruebele / Journal of Molecular Structure (Theochem) 589–590 (2002) 447–457

Fig. 7. Tiers 1– 3 of the IVR Hamiltonian of SCCl2, coded from black (tier 1) to light gray (tier 3) to highlight the gateway state structure near the bright state at Dn ¼ 0: Each stick represents a basis state; the y-axis for a state ljlis the sum over the product of the Lorentzian cutoff criteria defined in Section 2, for all paths leading P from lbl to ljl:Lbj for a state ljl in tier 1, i Lbi Lij for a state ljl in tier 2, etc.

delay dephasing from sub-ps to several ps. The necessary electric fields have an amplitude and phase structure (Fig. 6) that is well within reach of available pulse shapers. The key feature of Hvib that allows IVR to be frozen is the hierarchical structure of the vibrational couplings, as modeled here by Eq. (4) [16,23,37]. Fig. 7 shows the corresponding computed ‘gateway’ structure of states through which the initial wavepacket evolves. The resulting powerlaw diffusion of the wavepacket in the vibrational state space is evident in the uncontrolled IVR decay in Fig. 5: the rapid initial decay is followed by a much slower phase. The slow phase has also been observed in much larger organic molecules, such as cyclohexylaniline, and appears to be generic [17,18]. The complexity of the vibrational wavepacket (D in Eq. (1) and Fig. 1) always increases slowly with time, lengthening the controllable period of time. This interpretation is also supported by the shaped electric field: the maximum amplitudes in Fig. 6 roughly correspond to the IVR gateway states in Fig. 7. Because tiers of gateway states describe IVR well [38,39], one expects the coarse structure of 1(t ) to dominate the dynamics, as most of the freezing occurs by manipulating the flow between the bright state and the few states with direct, strong couplings to it. Only vibrations have been modeled here. In real

molecules, there will be contributions to dephasing from the environment, from rotational and Coriolis effects in the gas phase, and from vibronic interactions. The condensed phase environment poses a formidable problem, but rotational and Coriolis effects can be sidestepped in the gas phase. Large molecules, for which IVR control becomes interesting, rotate slowly in a cold molecular beam [7]. Most organic molecules have rotational constants , 0.1 cm21. An IVR control pulse, which cycles population between l0l and lbl, will cause undesirable rotational ladder climbing along with vibrational control. Eventually the maximum energy gap DE <  2 will exceed the degeneracy criterion DE < BJ "=2trx : However, starting with J ¼ 0; and trx < 10 ps, J ¼ 5 is required before rotational dephasing sets in (i.e. five coherent cyclings between l0l and lbl are allowed). Other effects are even smaller. Coriolis effects are linear in J and generally  BJ  2 for J . 1: Vibrational mixing tends to zCoriolis BJk reduce the dispersion of the rotational constants of eigenstates within an IVR feature [40], and dB ¼ P  2 =n1=2 also makes a negligible contri½ ni ðBi 2 BÞ bution. Thus, if control can be achieved within a few pump – dump cycles, rotations do not present a problem. By eliminating vibrational dephasing, one may be able to control reaction dynamics. The time scale of unimolecular dissociation reactions and of collision complexes is on the order of picoseconds above threshold [41]. A crossover of IVR and reaction time scales is thus possible on the 1– 10 ps time scale achieved in our model. Freezing IVR alone does not lead to a nonstatistical reaction outcome. The control process relies on Franck– Condon factors in its second stage. For example, if two different predissociative local mode states are coherently stabilized, then unlike highly mixed eigenstates, they would overlap substantially different sets of continuum states. This is quite analogous to Franck– Condon control of direct dissociation reactions by wavelength tuning [3], or of bimolecular reactions by vibrational excitation [2], except that highly vibrationally excited states in bound wells can now be stabilized and used. Partial Franck –Condon control has one important advantage over ‘coherent control all the way’. Direct formation of product states via coherent control requires degenerate paths, so that product states can

R.M. Bigwood, M. Gruebele / Journal of Molecular Structure (Theochem) 589–590 (2002) 447–457

interfere as t ! 1: When Franck –Condon control takes over from coherent control of IVR on a time scale trx, that stringent degeneracy can be relaxed. Pairs of eigenstates with splittings up to DE < "=2trx can be considered as coherently interfering because their mutual dephasing is negligible on the time scale of the Franck –Condon controlled reaction. We close with some speculation relating the freezing of IVR to current strong vs. weak field experiments. In weak field control, the molecular potential energy surface and molecular states are in the limelight; in strong fields, the AC Stark shift and field gradient broadening control the dynamics to zeroth order. Our proposal, of rephasing or stabilizing an IVR feature to give Franck– Condon control time to work, is basically weak field control. However, control of IVR may also be at work in recent strongfield experiments creating multiple ionized products. It is quite remarkable how long and structured the optimal pulses in strong-field control experiments are [8,9]. Because specific products are formed even with shaped pulses longer than the IVR time scale (< 100 fs), nuclear motions are very likely also subject to direct control in these experiments, although the final formation of the ion is by definition a strong field event. This suggests that a two-pulse approach might be best: a lower energy pulse to help guide the nuclear motions in preparation for reaction, and a much stronger and shorter pulse to induce the final product formation.

Acknowledgments This research was funded by grants from the NSF (CHE 94-57970 and 99-86670). M.G. gratefully acknowledges a Camille and Henry Dreyfus Scholarship. The authors thank M. Kellman and R. Levis for stimulating discussions.

References [1] W. Fuss, S. Lochbrunner, A.M. Muller, T. Schikarski, W.E. Schmid, S.A. Trushin, Chem. Phys. 232 (1998) 161. [2] M.J. Bronikowski, W.R. Simpson, B. Girard, R.N. Zare, J. Chem. Phys. 95 (1991) 8647. [3] F.F. Crim, Annu. Rev. Phys. Chem. 44 (1993) 397.

457

[4] D.J. Tannor, S.A. Rice, J. Chem. Phys. 83 (1985) 5013. [5] P. Brumer, M. Shapiro, Acc. Chem. Res. 22 (1989) 407. [6] S.A. Rice, M. Zhao, Optical Control of Molecular Dynamics, Wiley, New York, 2000. [7] V.D. Kleiman, L. Zhu, X. Li, R.J. Gordon, J. Chem. Phys. 102 (1995) 5863. [8] A. Assion, T. Baumert, M. Bergt, T. Brixner, B. Kiefer, V. Seyfried, M. Strehle, G. Gerber, Science 282 (1998) 919. [9] R.J. Levis, G.M. Menkir, H. Rabitz, Science 292 (2001) 709. [10] M. Gruebele, R. Bigwood, Int. Rev. Phys. Chem. 17 (1998) 91. [11] H. Rabitz, W.S. Zhu, Acc. Chem. Res. 33 (2000) 572. [12] S.A. Rice, Adv. Chem. Phys. 47 (1981) 117. [13] T. Uzer, Phys. Rep. 199 (1991) 73. [14] M. Gruebele, Adv. Chem. Phys. 114 (2000) 193. [15] S. Schofield, P.G. Wolynes, J. Chem. Phys. 98 (1993) 1123. [16] R. Bigwood, M. Gruebele, Chem. Phys. Lett. 235 (1995) 604. [17] M. Gruebele, Proc. Natl Acad. Sci. USA 95 (1998) 5965. [18] S.A. Schofield, P.G. Wolynes, J. Phys. Chem. 99 (1995) 2753. [19] V. Wong, M. Gruebele, J. Phys. Chem. 103 (1999) 10083. [20] D.M. Leitner, P.G. Wolynes, Chem. Phys. Lett. 280 (1997) 411. [21] R. Bigwood, M. Gruebele, Chem. Phys. Lett. 233 (1995) 383. [22] R. Bigwood, B. Milam, M. Gruebele, Chem. Phys. Lett. 287 (1998) 333). [23] D. Madsen, R. Pearman, M. Gruebele, J. Chem. Phys. 106 (1997) 5874. [24] R. Bigwood, M. Gruebele, D.M. Leitner, P.G. Wolynes, Proc. Natl Acad. Sci. USA 95 (1998) 5960. [25] A.M. Weiner, D.E. Leaird, G.P. Wiederrecht, K.A. Nelson, Science 247 (1990) 1317. [26] M.A. Dugan, J.X. Tull, W.S. Warren, J. Opt. Soc. Am. B 14 (1997) 2348. [27] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Hannery, Numerical Recipes in Fortran, Cambridge University Press, New York, 1992. [28] C.K. Chui, An introduction to wavelets, Academic Press, New York, 1992. [29] M. Gruebele, J. Phys. Chem. 100 (1996) 12178. [30] M. Gruebele, J. Chem. Phys. 104 (1996) 2453. [31] H. Tal-Ezer, R. Kosloff, J. Chem. Phys. 81 (1984) 3967. [32] P.B. Visscher, Computers Phys. (1991) 596. [33] S.K. Gray, J.M. Verosky, J. Chem. Phys. 100 (1994) 5011. [34] S.K. Gray, D.E. Manolopoulos, J. Chem. Phys. 104 (1996) 7099. [35] C.J. Bardeen, V.V. Yakovlev, K.R. Wilson, S.D. Carpenter, P.M. Weber, W.S. Warren, Chem. Phys. Lett. 280 (1997) 151. [36] B. Amstrup, G.J. Toth, G. Szabo, H. Rabitz, A. Lorincz, J. Phys. Chem. 99 (1995) 5206. [37] R. Bigwood, M. Gruebele, ACH Models Chem. 134 (1997) 675. [38] M. Quack, Faraday Discuss. Chem. Soc. 75 (1983) 359. [39] A.A. Stuchebrukhov, R.A. Marcus, J. Chem. Phys. 98 (1993) 6044. [40] B.H. Pate, J. Chem. Phys. 110 (1999) 1990. [41] R.D. Levine, R.B. Bernstein, Molecular Reaction Dynamics and Chemical Reactivity, Oxford University Press, New York, 1987.