Fully quantum coherent control

Chemical Physics 267 (2001) 33±46

www.elsevier.nl/locate/chemphys

Fully quantum coherent control M. Gruebele * Departments of Chemistry, Physics, and Biophysics, Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Received 18 September 2000

Abstract The coherent control problem is formulated in a second-quantized basis of system±®eld states (CC2Q). In principle, many system±®eld states are needed to compute the control target functional T. Here we show that a new and very simple formula for the fully quantized control functional can be obtained. As a result, T can be evaluated exactly without explicit knowledge of the system±®eld states. Two examples demonstrate the spirit of the approach, which can be scaled to larger problems by applying combined Lanczos and Hellman±Feynman methodologies. Although CC2Q, as a full quantum theory, does not scale as favorably as the semiclassical approach, its symmetrical treatment of system and ®eld promises new insights because it not only views the ®eld as a shaper of molecular coherence, but also the molecule as a shaper of ®eld coherence. Ó 2001 Elsevier Science B.V. All rights reserved.

1. Introduction The semiclassical formulation of coherent control treats the system and control ®eld on an unequal footing. The full quantum coherence of the system is taken into account, while the control ®eld enters classically as a time- or frequency-dependent parameter set which is to be optimized [1]. Besides controlling molecular reactions, quantum control is of interest in many other applications, such as dephasing/decoherence theory [2], quantum computing [3], and for probing the molecular Hamiltonian via input±output transformations of a probe ®eld [4]. In many of these applications information is coupled both into and out of the system via a ®eld, and it would be desirable to have a symmetric representation of the system±®eld

*

Fax: +1-217-244-3186. E-mail address: [email protected] (M. Gruebele).

which takes into account modi®cations of the ®eld structure by the control process. This becomes particularly important in applications with very low photon ¯uxes where semiclassical theory fails. At the same time, it would be useful if one could optimize the system±®eld interaction without the dicult constraint of a Schr odinger equation with a time-dependent Hamiltonian. Here we present a simple new formalism which provides a general control target functional for a system interacting with a quantum ®eld. The approach starts with the standard second-quantized system±®eld Hamiltonian. Because of the computational need for discretization, the quantum ®eld is treated in a ®nite volume. The resulting Hamiltonian is time-independent, and also no longer depends on the control parameters, which are now encoded in the initial and target density matrices. As a result, a surprisingly simple form of the control target functional is obtained: it is completely symmetrical in the initial and target density

0301-0104/01/$ - see front matter Ó 2001 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 1 - 0 1 0 4 ( 0 0 ) 0 0 4 0 4 - 3

34

M. Gruebele / Chemical Physics 267 (2001) 33±46

matrices, it is no longer represented by a trace, and it can be evaluated exactly without explicit knowledge of the system±®eld eigenstates or timeevolving density matrices. In order to put this new approach within computational reach for reasonably large system±®eld density matrices, we adapt procedures based on the Hellman±Feynman theorem, which have already proven themselves in large scale computations of molecular vibrational dephasing (IVR) [5]. To aid practical computational implementations, the Lanczos method [6] and basis set optimization by resonantly tiered ®ltering [2,7] are discussed. Two examples illustrate the method. A twolevel system is treated analytically to illustrate the spirit of the formalism. A larger numerical calculation shows how the method can be applied to more realistic systems to optimize control ®elds or extract molecular information in the presence of spontaneous emission.

matrix which started at q0 , and the target density matrix qT . Although the integration is taken between 1, the ®eld e…t; c† in practice limits the duration of the control. The general approach of feedback control which can be applied to Eqs. (1)± (3) has been well developed in engineering and molecular systems [4]. The semiclassical formulation can also be expressed in the spectral domain [9±11]. In that form it explicitly recognizes that coherent control in the limit t ! 1 requires multiple paths connecting three or more molecular energy levels, two or more of which generally belong to a degenerate product manifold. Rather than using e…t; c†, the control ®eld in the spectral representation is de®ned by amplitudes and phases of the radiation mode …k; k†, where k is the wave vector and k the polarization index: e ˆ e…A…ki ; ki †; u…ki ; ki ††;

…4a†

2. Second-quantized control

e ˆ e…A…xi †; u…xi †† ˆ e…Ai ; ui †:

…4b†

2.1. The basic formalism

Case (4b) arises when the light is collimated and in one polarization in vacuo, such that xi ˆ cjki j suciently speci®es the mode. The set of control parameters is ‰cŠ ˆ ‰Ai ; ui Š. The frequencies xi (or more generally ki and ki ) are usually de®ned by resonance, phase or other matching conditions, and not treated as control parameters, although one may do so (see Section 3.1). The target problem in the spectral case is usually treated by a perturbation theory (e.g. third order for 1 ‡ 3 photon interference), and the enforcement of resonance allows a ÔdressedÕ state approach. A fully quantized version of the coherent control problem can be obtained by replacing Eq. (1) by the well-known second-quantized Hamiltonian. At a ®rst glance, it would seem that the resulting target function requires prohibitively large basis set calculations due to the many quantum ®eld modes. Instead as we shall see, the target functional can be manipulated into a surprisingly simple form which is amenable to computation with fairly sizeable systems. As will be seen, these advantages arise because the second-quantized Hamiltonian satis®es energy conservation and no longer depends on the control parameters. The resulting for-

The common formulation of the coherent control problem uses a time-dependent semiclassical Hamiltonian [1,8] H …c† ˆ Hsys

le…t; c†;

…1†

where Hsys is the quantum Hamiltonian of the system to be controlled, l is the system coordinate-dependent dipole operator, and e…t; c† is the control ®eld with control parameters c. The parameters c are adjusted to maximize the target functional T ‰H …c†Š,   Z t  Z 1 T ‰H …c†Š ˆ dt Ot Tr qT exp i dt0 H …t0 ; c†= h 1 1   Z t 0 0  q0 exp i dt H …t ; c†= h ; 1

…2†

subject to a set of constraints on the parameters c, ffi …c† ˆ 0g;

…3†

which can be enforced by Lagrange multipliers. T is the overlap between the time-evolving density

M. Gruebele / Chemical Physics 267 (2001) 33±46

mulation is time-independent and the full system± ®eld states are ÔdressedÕ, but unlike the standard multiple path-interference formulation of the previous paragraph, it treats the molecular and ®eld degrees of freedom on a completely equal footing. We begin with the standard second-quantized Hamiltonian analogous to Eq. (1) [12]: …1†

…2†

H ˆ Hsys ‡ Hfield ‡ Hint ‡ Hint X ˆ Hsys ‡ hxk;k …ayk;k ak;k ‡ 12†  k;k

1=2 X qn  2p h pn ek;k …e m vxk;k n;k;k n ‡e

ik
rn

ik
rn y ak;k

…2†

ak;k † ‡ Hint :

…5a†

In this form, all types of excitation and emission processes can be treated, and phase matching is explicitly taken into account via the factors eik
rn : in principle, one could compute the three-dimensional structure of the output ®eld based on the molecular trans±rot±vibronic wave functions and a full ®eld basis. For later computational convenience, we adopt a formulation in a cavity of ®nite volume v. For the purpose of molecular control, it is often convenient to make the dipole approximation X H  Hsys ‡ hxk;k …ayk;k ak;k ‡ 12†  X ‡i k;k

k;k

1 …vxk;k †

1=2

^Šek;k …ayk;k ‡ ak;k †: ‰Hsys ; l

…5b†

^ is the system dipole opIn this representation, l erator, k, k sum over wave vectors and polarizations, and xk;k and ak;k are the ®eld frequencies and ladder operators. Neglecting the quadratic vector …2† potential term Hint means that even-photon scattering processes are not exactly represented. Setting the phase factors eik
rn  1 means that only near-®eld (``thin ®lm'') phenomena can be computed, unless phase matching is imposed as an additional constraint. When required, the additional matrix elements needed for Eq. (5a) are straightforward to evaluate; none of what follows in the derivation of the second-quantized control target functional depends on the choice of Eqs. (5a) vs. Eq. (5b). Note that Eqs. (5a) and (5b) are

35

not in an interaction representation, so there is no explicit time dependence. We now turn to the derivation of the target functional T using Eqs. (5a) and (5b). Eigenstates jNi of H will be used for the purpose of deriving a formal expression for T. In practice, these eigenstates are not available for systems of interesting size: if fjsig are the eigenstates of Hsys and fjnk;k ig are number eigenstates of Hfield , a reasonable representation of the ®eld degrees of freedom requires a large product basis fjsijnk;k ig. Full diagonalization of H in such a basis to obtain eigenvectors is not possible. The issues associated with circumventing the use of eigenstates in an exact formalism, and with the choices of ®eld bases (e.g. number eigenstates jnk;k i vs. coherent states j n; ui) will be addressed in Sections 2.2 and 2.3. In terms of the Hamiltonian (5a) and (5b) Eq. (2) for the target function can be rewritten as Z 1 T ‰q…c†Š ˆ dt TrfqT …c†eiHt=h q0 …c†e iHt=h g: …6a† 1

There are two important di€erences between Eqs. (2) and (6a). The Hamiltonian is no longer timedependent, so application of the propagator is straightforward and requires no time ordering operator Ot . The Hamiltonian no longer explicitly depends on amplitude or phase control parameters c, which are now associated with the initial and target density matrices. For simplicity, from now on we will write T ‰cŠ, with the understanding that the c-dependence comes through q. (H still explicitly depends on the frequencies xk;k , if those are considered to be control parameters.) Although one may customarily think of the parameters c as associated with the initial density matrix q0 , the same parameters appear in the target density matrix qT . The control problem in this representation is thus completely symmetrical, as can be seen by rewriting Eq. (6a) in the equivalent form Z 1 T ‰cŠ ˆ dt Trf…e iHt=2h qT …c†eiHt=2h † 0

Z ‡

0

 …eiHt=2h q0 …c†e 1

dt Trf…e

iHt=2 h

 …eiHt=2h qT …c†e

iHt=2 h

†g

q0 …c†eiHt=2h † iHt=2 h

†g:

…6b†

36

M. Gruebele / Chemical Physics 267 (2001) 33±46

It is also clear from Eq. (6a) or (6b) that the target functional is independent of time due to the integration, and that the control problem should be expressible directly in terms of the parameters c, without any time evolution. Eq. (6a) can be expressed in terms of Green operators by using the Laplace transform representation of the propagator e

iHt= h

I 1 h…i hk H † dk ekt i 2pi I 1 dk ekt i ˆ hG…ihk† 2pi ˆ

1

…7†

which results in Z 1 Z 1 Z 1 1 0 T ‰cŠ ˆ 2 dt dE dE0 e i…E E †t=h 4p 1 1 1  TrfqT …c†G… E0 †q0 …c†G…E†g: …8† Integration over t yields a d-function which eliminates the integration over E0 such that Z 1 1 T ‰cŠ ˆ dE TrfqT …c†G… E†q0 …c†G…E†g: 2p 1 …9† The last integral can be evaluated by inserting full product basis identities into Eq. (9) to yield T ‰cŠ ˆ

X 1 hNjqT …c†jMihMjq0 …c†jN i 2p N;M Z 1 1  ; dE …E E †… E EM † N 1

…10†

and by taking care to use the proper asymptotic expressions, which eliminate one of the two summations. The ®nal result obtained formulates the target functional in terms of the initial and ®nal density matrices and the total energies T ‰cŠ ˆ

1 X hNjqT …c†jNihNjq0 …c†jNi : 2 N jEN j

…11†

Eq. (11) is the principal result of this paper. It has several features that will prove useful:

· T ‰cŠ depends on the absolute value of the total eigenenergies EN , not on the system energies. This makes it particularly easy to enforce resonance conditions via overall energy conservation. In the limit of a continuous quantum ®eld energy is rigorously conserved and the denominator may be taken out of the sum. It may also be shifted by an o€set without disturbing the maximum of T ‰cŠ. · T ‰cŠ is no longer a trace; it depends on the product of diagonal elements of the density matrices in the system±®eld basis. This is critical for practical applications (Section 2.2). · Coherent interference is encoded in Eq. (11) because q0 and qT in general do not commute with H, so that the phase parameters in c enter into T ‰cŠ. · T ‰cŠ preserves the complete initial state±target state symmetry already noted in Eq. (6b). · T ‰cŠ depends on the control parameters without the constraints of forward or backward temporal propagation, so searches for optimum ®eld amplitudes and phases are greatly simpli®ed compared to optimal control of the time-dependent Schr odinger equation. Other types of constraints (such as a penalty for excessive ®eld amplitudes) can still be implemented by Lagrange multipliers. · The system±®eld information is now solely encoded in the density matrices in terms of molecular parameters and ®eld control parameters c. Nonlinear or even powerful linear variational methods can therefore be used to e€ect optimization. · The above derivation was performed using discrete-state notation, but the generalization to a continuum where the upper end of the sum is replaced by an integration is straightforward. In practice, the continuum could be replaced by a discretized or ®nite-degeneracy representation, just as in the case of standard time-dependent or independent approaches. 2.2. Computing the density expectation values Eq. (11) provides a very compact formulation for the target functional T ‰cŠ. However, there is a

M. Gruebele / Chemical Physics 267 (2001) 33±46

road block to its application: it apparently requires knowledge of the eigenfunctions of H to permit evaluation of the density matrix expectation values. Even in a discretized model, a very large number of photon states jnk;k i or jai has to be included in the photon basis. While it is possible to calculate eigenvalues of matrices with N > 106 basis states relatively easily, repeated computation of all the eigenvectors is out of the question. A solution requiring eigenvalues only must therefore be obtained, at least for the optimization step. A related problem occurs in the vibrational dephasing of molecules at high energies (IVR): large numbers of vibrational basis states interact to produce the envelope of spectral features, making the dephasing dynamics dicult to calculate. The IVR problem can be solved by the MFD theorem [5,13,14], an application of the Hellman± Feynman theorem hf j

oEf0 oH 0 jf i ˆ ; ok ok

…12†

where jf i are the eigenfunctions of H 0 . MFD evaluates the spectral envelope without knowledge of the eigenstates. In combination with Lanczos diagonalization, this method has evaluated the dephasing of systems with 6±63 degrees of freedom and >105 basis states [7,15,16]. It has also recently

37

been used to assign quantum numbers to spectra without knowledge of the eigenfunctions [17]. Consider how Eq. (12) can be used to obtain information useful for systems undergoing complex quantum dynamics as a result of interacting with a radiation ®eld. (a) MFD theorem: in its original form, it applies to a single excited state carrying oscillator strength, which mixes with a set of dark background states coupled amongst one another. The resulting coupling matrix has the form shown in Fig. 1. Consider the modi®ed Hamiltonian H 0 ˆ H ‡ k…PHQ ‡ QHP †;

…13†

where P ˆ j0ih0j is the bright state projector and Q is its complement. Inserting into Eq. (12) yields 1 1 oEf If ˆ jhf j0ij2 ˆ …14† 2 Ef E0…0† ok kˆ0 which expresses the intensities in the spectrum in terms of the eigenvalues only; a simple ®nite difference derivative with two eigenvalue evaluations is sucient to calculate the entire spectrum, and formulations with a single evaluation are also possible. Eq. (14) is particularly useful because its structure resembles a quantum version of the ¯uctuation±dissipation theorem. Any Hamiltonian ^ will yield the expectation of the form H ‡ kO

Fig. 1. (A) Partitioning of H to compute the spectrum of a state j0i, which interacts with a coupled manifold of states QHQ, via the MFD theorem. P is the projector for j0i and Q its complement. (B) Partitioning of H if the full Hamiltonian is blocked, for example into electronic states. The block Hint in that case contains optical transition operators (and nonadiabatic interactions, if needed. Other bi and tripartitionings are recognized as the Tannor±Rice or Brumer±Shapiro schemes discussed in the text.

38

M. Gruebele / Chemical Physics 267 (2001) 33±46

^ For example, the even simpler modivalues of O. ®ed Hamiltonian H ‡ kP yields If ˆ oEf =dk by modifying a single matrix element; the modi®ed Hamiltonian H ‡ k^ ljgihgj^ l connecting a state jgi to a manifold of strongly coupled states fjnig yields the full absorption spectrum from jgi to fjnig without any knowledge of the eigenstates jf i other than their eigenvalues. (b) Density matrices: by choosing H ‡ kq0 ‡ k0 qT , one can likewise reduce Eq. (11) to T ‰cŠ ˆ

1 X …oEN ‰cŠ=ok†…oEN ‰cŠ=ok0 † ; 2 N jEN ‰cŠj

…15†

eliminating all reference to the eigenstates jNi of H. Eq. (15) is the principal computationally oriented result of this paper: the target functional T ‰cŠ can now be evaluated by a small number of eigenvalue computations (3±5 depending on the method used). This is particularly useful because the density operators in question are usually very sparse operators in a basis of the type fjnijnk;k ij n; uig, but not very sparse in the eigenbasis. Multiple target evaluations up to N  105 required by optimization algorithms are thus feasible on modest computers, and >106 on large computers. 2.3. Large systems: Lanczos algorithm, basis sets and basis ®ltering Although a detailed discussion of all the computational approaches that can be brought to bear on Eq. (15) is beyond the scope of a single paper, some comments on practical aspects of working with Eq. (15) are required. For computations on large systems, ecient eigenvalue algorithms and ecient representations of the basis space are required. We discuss these in turn. The eigenvalue algorithm must be of high accuracy, so the derivatives required by Eq. (15) can be computed reliably. The Lanczos algorithm satis®es this requirement [18]. Other algorithms we have tested [19] are suitable for the calculation of spectra, but do not provide sucient accuracy for the evaluation of central di€erence derivatives with k  0:001 as required by Eq. (15). Lanczos iteration su€ers from one potential drawback in con-

nection with Eq. (15). Because the total system± ®eld basis states generally conserve energy by construction, the eigenvalue problem is very sti€ in the limit of weak ®elds, where spectral degeneracies are not signi®cantly lifted. A large number of iterations is required for convergence in such cases. This problem may be alleviated somewhat by using a spectral ®lter such as G0 …E† ˆ …E Hsys 1 Hfield † to precondition the Lanczos procedure [20]. For weak ®elds, one can also use degenerate perturbation theory to compute the necessary energies. As written in Eqs. (5a) and (5b), the Hamiltonian is completely general. In many cases, it is convenient to partition the Hamiltonian, for example when the Born±Oppenheimer approximation de®nes electronic states. Eq. (11) lends itself very easily to this kind of partitioning. For two electronic states, Fig. 1b shows the general form of the Hamiltonian matrix in a zero-order basis (e.g. fjnijnk;k ig from Section 2.1): two blocks H1 and H2 near the diagonal contain the couplings within system states 1 and 2, while the two o€-diagonal blocks contain the system±®eld couplings connecting 1 and 2. In the zero-order basis, the density matrices q0 and qT are similarly disposed: q0 typically has nonzero elements only in the H1 block, while qT may have nonzero elements in both blocks depending on the target to be reached. In the case of a Tannor±Rice scheme [21], q0 would be nonzero only in the thermally accessible block of H1 , qT only in the upper energy block of H1 ; the H2 block of q0 and qT would be zero in the zero-order basis, although all hNjqjNi in the eigenbasis would be nonzero because the optimization process accesses zero-order states from the H2 block as ÔintermediatesÕ. In a Brumer±Shapiro scheme [10], H could be conveniently tripartitioned with q0 nonzero in the H1 block, qT nonzero in the H3 block, and the H2 block containing the ÔintermediatesÕ which allow interference. A numerical example is given in Section 3.2. For a spectral representation of H, the choice of basis becomes very important. When many ®eld modes are involved in a control scheme, the dimensionality of the Hilbert space required by Eq. (15) can take on astronomical proportions. There are two interesting limits that can be distinguished.

M. Gruebele / Chemical Physics 267 (2001) 33±46

At high photon ¯ux, the usual case in molecular control, the number of potential basis states could be enormous. As seen below, certain propensity rules become nearly rigorous selection rules in this limit, allowing drastic pruning of the space of system±®eld basis functions. Low photon ¯ux is of interest in experiments involving quantum information transfer into and out of systems, unobtrusive measurement, and single or few molecular measurements such as in near-®eld coherent studies of molecules on surfaces. In such applications, propensity rules are rather weak and fairly complete bases must be included, especially if the quantum ®eld e€ects observable only at low photon ¯uxes are to be properly taken into account. Fortunately, the low photon numbers mean that small bases can span each degree of freedom, allowing an adequate discretized representation. Most of the discussion of basis sets and basis ®ltering will therefore concentrate on the high-®eld case. A problem also arises when uncollimated light sources such as near-®eld apertures are used. Unlike high-®eld collimated laser experiments, near-®eld interactions do not involve a coherent input ®eld with a single wave vector k, and large phase shifts can occur in the ®eld upon interaction. Several basis sets have already proved convenient in literature discussions of quantum radiation [12]. In all of these basis sets, individual ®eld states have a photon amplitude n1=2 and a phase u associated with the mode …k; k†. In practice a mixture of di€erent basis sets, adapted to the properties of di€erent ®eld degrees of freedom, will yield the best results. In the following discussion we will label modes by either …k; k†, or by x ˆ cjkj, or suppress the label altogether when there is no danger of confusion. The simplest ®eld basis is the number basis. These are eigenstates of Hfield with coupling matrix elements hn…k; k†  1jaykk ‡ akk jn…k; k†i p ˆ …n…k; k† ‡ 1=2†  1=2:

…16†

Their ®rst advantage is orthogonality. Their second advantage is that they allow full control over phase and amplitude of aP®eld mode via linear combinations of the type cn jni. The cn are

39

complex amplitude/phase control parameters, one for each basis function in the sum. The amplitude jcj which adjusts the contribution of each basis state jni is not to be confused with the photon amplitude n1=2 of the mode …k; k†, which is given by P n ˆ njcn j2 . The number basis is very convenient for low photon ¯uxes, but very inconvenient at high photon ¯uxes because individual number states do not even approximately represent a high intensity laser ®eld. Rather, they correspond to a coherent radiation state highly squeezed in amplitude at the expense of phase. To represent unsqueezed input laser ®elds, a particular linear combination of number states is very useful, the coherent state j n…x†; u…x†i ˆ

1=2 n 1 X … n…x† eiu…x† † p jn…x†i n! nˆ0

…17†

which o€ers minimum uncertainty in both phase and amplitude, and corresponds to a normal monochromatic laser beam in the classical limit. The properties of this basis are well documented [12, 22], but those particularly pertinent to the evaluation of Eq. (15) are worth reiterating. The main concern is that the coherent states are not orthogonal, and that the operator H…1† couples each coherent state to an in®nite manifold of other coherent states, 1=2

h n0 u0 jay ‡ ajnui ˆ … n1=2 eiu ‡ n0 e e

n0 =2 n=2‡ei…ux

iu0

†

u0x †

p nn0

:

…18†

Fortunately, the situation simpli®es greatly for large n. First, H …1† couples only coherent states with n0 ˆ n  1, corresponding to energy conservation in the system plus ®eld upon emission or absorption. This still leaves the possibility of large phase jumps u ! u0 ˆ u ‡ Du in the coherent state. Fig. 2 shows the magnitude of the coupling matrix element (18) as a function of n and Du. The coupling increases as the root of n, and at the same time sharpens rapidly in phase as e

nDu2 =2

:

…19†

For large photon numbers the phase change is thus suppressed, and if energy conservation is enforced in the basis, only one coherent state is

40

M. Gruebele / Chemical Physics 267 (2001) 33±46

Fig. 2. Magnitude of the matrix element of H…1† coupling two coherent states as a function of n and Du. The matrix element obeys a propensity rule Du ˆ 0 for large n. The size distributions at n ˆ 100 and 1000 are shown as thicker curves superimposed on the contour surface.

necessary for each transition. This automatically solves the orthogonality problem also, because orthogonality of the overall system ®eld product basis is maintained by the system wave functions. The phases u…x† for di€erent modes, as well as the mode intensities n…x† can of course be adjusted as control parameters even in this simple case. For very low photon ¯uxes, several coherent states must be included for each system state to allow for phase jumps; in that limit, it may be more useful to express the entire problem in terms of number states. Otherwise, a generalized Lanczos diagonalization HM ˆ MSK must be solved, where S is the nondiagonal overlap matrix caused by the coherent basis states. It is also possible to represent pulses with arbitrary phase and amplitude characteristics using coherent states. In general jP …Ai ; ui †i ˆ

N Y j n…xi †; u…xi †i;

…20†

iˆ1

where each frequency xi (or generally …ki ; ki †) corresponds to a di€erent ®eld degree of freedom,

resulting in a product form. In practical calculations the product can include only a ®nite number of terms, or even a very sparse ®ltered basis (see below) becomes unmanageable. A ®nite product (20) corresponds to a train of pulses, rather than a single pulse. In experiments, pulse shapers in any case allow control over only a ®nite number of frequency components of the ®eld. Note that the product form of Eq. (20) implies that the orthogonality problem is not exacerbated in pulses compared to monochromatic states because each frequency component belongs to a separate degree of freedom. Pulses are not limited to expansion in terms of coherent states. Number states, or any other type of squeezed state, can also be used in the basis depending on the application and type of input ®eld. Lastly, but very importantly especially in the low ¯ux limit, comes the representation for spontaneous emission modes, which can be treated by the quantum ®eld Hamiltonian in Eq. (5). An exact treatment requires a spontaneous emission manifold jspi ˆ

Y k;k

…sp†

jnkk i 

N Y jn…sp† xi i

…21†

iˆ1

compared to which the laser manifold has small or zero measure, in keeping with the classical intuition that entropy should increase when laser radiation is converted by absorption and spontaneous emission from well-de®ned …k; k† to wave vectors distributed over 4p steradians and all polarization components. Although the large number of degrees of freedom involved in the spontaneous emission component (21) of the total wave function generally presents diculties, it has one advantage: the manifold fjspig and ÔlaserÕ degrees of freedom like jni or j n; ui are essentially orthogonal (except for a set of measure zero), and can be treated as such. The manifold in Eq. (21) cannot be properly reduced to one degree of freedom if long-time dynamics are required. The usual treatment is to develop a master equation for the density matrix, which decays exponentially according to Wigner± Weisskopf [22]. Such an approach could indeed be used if Eqs. (6a) and (6b) were applied to the

M. Gruebele / Chemical Physics 267 (2001) 33±46

control problem. For an approximate solution using Eq. (11), jspi can be approximated by a small number or even a single product basis function, as shown on the right of Eq. (21). If the spontaneous emission lifetime for a system transition from jni to jn 1i is known, then we can approximately replace the spontaneous emission part of the coupling matrix element by  1=2 N X h  …1† …sp† …sp† hnxi 1jHint jnxi i ˆ …22† 2pssp qeff iˆ1 where the xi uniformly cover the linewidth of the jni to jn 1i transition and hqeff ˆ …xi‡1 xi † 1 . The matrix element (22) is taken to be independent of n…sp† because in the real manifold jspi the probability of exciting the same …k; k† twice is vanishingly small. Eq. (22) is approximately valid only for control ®elds that achieve their goal in the temporal representation in t <  hqeff =2; in other words, the Hermitean conjugate matrix element of Eq. (22) will eventually put spontaneously emitted photons back into the system, unlike the manifold on the left of Eq. (21), which corresponds to an in®nite product wave function and a continuous radiation bath. A more elegant explicit continuum solution to the problem posed by spontaneous emission may well exist, but for now Eq. (22) will suce to deal with spontaneous emission. A product basis for computing with Eq. (11) must include system, input/output modes and spontaneous modes, and thus has the general form 8 9 ( ) 0 N N < Y = Y X jsi cij j nj …xi †; uj …xi †i jn…sp† i …23a† xk : iˆ1 ; j kˆ1 if (for example) coherent states are used to represent the control ®eld. In this form the c ˆ ‰cij Š can serve as linear control parameter, and the nj and uj are ®xed. A very ecient linear variational principle can be applied for optimization of the cij , but at the expense of basis set size. Alternatively, one can reduce Eq. (23a) partially or completely to a form 8 9 N N0 < Y = Y jsi j n…xi †; u…xi †i jn…sp† …23b† xk i : : iˆ1 ; kˆ1

41

In this case, the control parameters c are the n…xi † and u…xi † themselves, and generally a nonlinear variational principle must be applied to optimize the ®eld. One can also consider the frequencies, or even the molecular constants in Hsys as control parameters, although they are usually ®xed by resonance conditions or knowledge of the Hamiltonian. Nevertheless, the possibility of modifying Hsys or Hfield shows that information can be coupled into and out of both the ®eld and the system. For example, modifying the input ®eld or changing Hsys (e.g. via the Stark e€ect) can provide detailed information about Hsys in the output ®eld. To illustrate the need for basis set ®ltering, consider the six-level example shown in Fig. 3. The ®eld frequencies are labeled by the system states they nominally connect. A minimal basis could include: · N ˆ 6 resonant frequencies in the control ®eld. (This assumes x23 , x24 , x34 , x15 , x16 lie outside the bandwidth of the control ®eld; it also assumes that a symmetry enforces x25 ˆ x26 , x35 ˆ x36 , x45 ˆ x46 , at the ®eld strengths of interest.) · N 0 ˆ 18 modes for the six di€erent resonant spontaneous emission frequencies if spontaneous emission is to be treated approximately. (This assumes that x23 , x24 , x34 , x15 , x16 emission is negligible because of long wavelength or selection rules, and that x25 ˆ x26 , x35 ˆ x36 , x45 ˆ x46 , but gives the remaining six emission frequencies, three xi in Eq. (22) to mimic coupling of the system to the radiation continuum.) We have 24 radiation degrees of freedom total. If the spontaneous emission degrees of freedom are treated with only two basis functions each, the system with six, and each control degree of freedom with four (to allow for nonzero Du for example), there are already 61  46  218 > 6  109 product basis functions, impractical for iterative optimization on any present computer. To prune this product basis drastically without a great loss of accuracy, two principles can be invoked: energy conservation, and the selection and propensity rules discussed above.

42

M. Gruebele / Chemical Physics 267 (2001) 33±46

Fig. 3. Basis selection by tiering for a six-level system. Starting with the initially populated state in tier s ˆ 0, tier s ‡ 1 contains all states coupled to states newly generated in tier s which satisfy the selection and propensity rules, and which lie in a speci®ed energy window. Sample states generated at each step are shown at the bottom. By the third tier, a minimal basis set allowing all possible processes at least once has been constructed. The procedure may be continued in the same manner to enlarge the basis. In the text, it is assumed that x25 ˆ x26 , etc. This assumption may have to be relaxed by adding three additional ®eld modes, if the system±®eld interaction can break the symmetry which causes the degeneracy of j5i and j6i in the absence of the ®eld. Note that in principle upwards Ôspontaneous emissionÕ could also occur in this purely conservative model. Its suppression is enhanced by the choice of interaction operator in Eq. (22) and by increasing the number of spontaneous ®eld modes in the calculation.

Energy conservation holds strictly for the eigenstates jNi of the full Hamiltonian. They form a continuum of system±®eld states in the limit v ! 1 in Eqs. (5a) and (5b), and an initial state at energy EN can only interact with degenerate continuum states in the limit t ! 1. In that case, one can replace Eq. (11) by T ‰cŠ ˆ

1 X hNjqT …c†jNihNjq0 …c†jNi: 2jEj N

…24†

For computational purposes, this principle approximately applies also to the basis states which are eigenstates of Hsys ‡ Hfield . For ÔweakÕ ®elds, the eigenvalues are not shifted very much compared to the zero-order energies. One can therefore truncate the basis into an energy window DE which can be of the order jl
Ej. Considerable pruning can be achieved as long as DE   hxij . ÔWeakÕ is a relative term here: dynamical Stark shifts even for rather strong ®elds are generally much smaller than the energy of near-IR or higher frequency photons. For example, a two-level system interacting with a single photon yields a product basis fjs ˆ 0ijn ˆ 0i, j0ij1i, j1ij0i, j1ij1ig, but energy conservation allows this to be pruned to fj0ij1i, j1ij0ig for all but the strongest ®elds because the molecule is

either in its ground state and the photon in the ®eld, or the photon is gone and the molecule in its excited state. Pruning by energy conservation is most easily enforced together with applicable selection and propensity rules by considering the system±®eld interaction as a tiered process. Each tier is associated with a new set of states connected to the previous tierÕs states by H…1† , and produces a new set of control parameters. The example in Fig. 3 illustrates this procedure. Let the initial state be …sp† j1iPj nij , uij iPj0kl i. The system is in its ground state interacting with a coherent ®eld of 6 modes, and no photons are in the spontaneous background ®eld yet (i.e. xij  kT ). In the ®rst tier one can directly generate three more basis functions which satisfy energy conservation, selection rules, and the Du ˆ 0 propensity rule. (3k ‡ 3 basis functions can be generated if k phase shifts kDu0 =2 . . ., Du0 , Du0 ,. . .kDu0 =2 are allowed per mode). From there, we can generate second and third tier basis functions. At this point all possible processes have been included once, so the minimal basis has only about 100 states …k ˆ 0†. Even for k ˆ 4, the basis has only about 1000 basis functions, compared to >6 109 for the direct product basis. If Du 6ˆ 0, some of the coherent

M. Gruebele / Chemical Physics 267 (2001) 33±46

state phase shifts that can occur during multiple ®eld cycles are already taken into account by the ®rst three tiers. To increase the basis size systematically further, the cycle of tiers can be continued with the most recently generated set of new states to include additional phase shifts, or k, N or N 0 can be increased to test for convergence. 3. Two examples 3.1. Interaction of a two-state system with resonant number states We start with a very simple analytical example which illustrates the procedure for a well-known system±®eld [12]. Consider a two state system with eigenstates fj0i, j1ig and energies …0; hx0 †. Let this system interact with a number state ®eld fjnx ig of frequency x polarized along the dipole orientation of the system. If we take into account only a single cycle of interaction with a near-resonant monochromatic ®eld, the full basis set fj0i; j1ig fjnx ig can be ®ltered down to only two basis functions, fj0ijnx i, j1ijnx 1ig. In essence, this is analogous to the single frequency rotating wave approximation in the semiclassical case (although the number ®eld is highly squeezed). The Hamiltonian matrix is then given by

43

Let the system initially be in the ground state, and let the goal be maximal transfer of population by the number state ®eld:     1 0 0 0 and qT ˆ : …26† q0 ˆ 0 0 0 1 The Hamiltonian matrix modi®ed according to Section 2.2(b) then becomes  p  hDx=2 ‡ k le0 nx p H0 ˆ …27† le0 nx hDx=2 ‡ k0 with eigenvalues k k0 1 E ˆ ‡  …k2 ‡ k02 ‡ h2 Dx2 2khDx 2 2 2 1=2 ‡ 2k0 hDx 2kk0 ‡ 4nl2 e20 † : From this follows q 1 h2 Dx2 ‡ 4nl2 e20 jE jjk;k0 ˆ0 ˆ 2 oE =okjk;k0 ˆ0 ˆ hjq0 ji ˆ

1 1 hDx  q 2 2 h2 Dx2 ‡ 4nl2 e2 0

oE =ok0 jk;k0 ˆ0 ˆ hjqT ji ˆ

1 1 Dx  q : 2 2 h2 Dx2 ‡ 4nl2 e2 0 …29†

…1†

H ˆ Hsys ‡ Hfield ‡ Hint     0 0 hx…n ‡ 1=2†  0 ˆ ‡ 0 hx0 0 hx…n 1=2† p   0 le0 nx ‡ p le0 nx 0 p   hDx=2  le0 nx ˆ ^ ; …25† p le0 nx hDx=2  where Dx ˆ x0 x, and the average energy has been subtracted from the diagonal in the ®nal expression. The two ®eld control parameters are c ˆ ‰nx ; xŠ. The amplitude control parameter nx enters the matrix through the choice of basis; x is often ®xed by a resonance condition and not considered a control parameter, but here we will also use it as a control parameter. There is no phase parameter in this simple example.

…28†

Inserting into Eq. (15) the target function becomes T ‰nx ; xŠ ˆ

2nl2 e20 h2 Dx2 ‡ 4nl2 e20


3=2 ;

…30†

which is obtained without any explicit consideration of the eigenstates. T ‰cŠ is plotted in Fig. 4. For ®nite nle20 , the target is maximized when Dx ˆ 0. This corresponds to the usual resonance condition. At nle20 ˆ 0, there is a node in T because no population transfer can be achieved without a dipole moment. Although T ˆ 0 when nle20 ˆ 0, the maximum target is achieved for in®nitesimal values of nle20 . There are two reasons for this: T drops o€ rapidly as nle20 increases because resonance is

44

M. Gruebele / Chemical Physics 267 (2001) 33±46

e€ected by a single-phase monochromatic number ®eld (a p pulse is necessary, which requires a larger ®eld basis space). As a general precaution, maximization of T does not guarantee that qT is reached, only that it is maximally overlapped by the density matrix given the available control parameters. The actual density matrix can be evaluated by a single propagation of Eq. (6a), for example. 3.2. A more realistic system Fig. 4. Target surface T ‰n; DxŠ for the two level system in Section 3.1. There is a node at zero n (or l), but in the absence of spontaneous emission yet presence of a dynamical Stark shift, the maximum is reached for in®nitesimal n.

Fig. 5. Eigenvalue surface for the Hamiltonian in Eq. (25). The dressed states become nondegenerate as a function of Dx or the Rabi frequency.

destroyed by the dynamical Stark e€ect, as illustrated by the plot of E in Fig. 5; and because no spontaneous emission is allowed in this example, even a weak ®eld can induce maximum population transfer. It should be noted that when T is maximized at Dx ˆ 0 and nle20 ! 0, qT is not actually reached by application of the simple number eigenstate ®eld. Instead, p the eigenfunctions of H become ji ˆ …1= 2†…j0i  j1i† and the corresponding density matrix has at best a 0.5 overlap with qT : full population transfer can simply not be

For larger systems, the basis set ®ltering and eigenvector-free computation of T ‰cŠ as a function of the control parameters have to be performed numerically. The F95 code CC2Q has been written to implement the ideas described in Section 2. Calculations for the six-level system of Fig. 3 will be discussed here. Table 1 summarizes the system energies, transition moments, and the ®eld control parameters for each coherent input mode. The initial state was chosen at kT ˆ 0 for the system. The target qT was chosen as a diagonal matrix containing all basis states with jWsys i ˆ j6i with equal weight, irrespective of the ®eld modes. This target will optimize the population in state j6i without giving preference to any speci®c output ®eld. A triply iterated basis with N ˆ 6 and N 0 ˆ 6 containing 30 basis functions was generated for the optimization as discussed in Section 2.3, Eq. (23b). It was assumed that the relative phases of the 1 ! 2 and 1 ! 4 transitions could be adjusted independently, but the frequencies of the two modes were held constant. A nonlinear iterative simplex search of the parameter space was applied to maximize the target as a function of u12 and u14 . More powerful search algorithms could be used to better avoid false minima; even linear optimization could be used with a larger basis (see Section 2.3, Eq. (23a)). Local minima are not a problem in this case. For this six-level system, the CC2Q method allows the entire coherent control target surface to be simply fully evaluated. Fig. 6 shows a twodimensional cut through the full optimization surface T near the optimum in ui1 and u14 . As expected, the target function is sensitive to the phase parameters of ®elds involved in di€erent

M. Gruebele / Chemical Physics 267 (2001) 33±46

45

Table 1 Energies, transition dipoles and ®eld parameters for the calculation in Section 3.2a ~ sys; x

i

1

1 2 3 4 5 6

0 cm 498 cm 1 501 cm 1 502 cm 1 1020 cm 1 1020 cm 1

I0 ˆ 1 J cm2 V ˆ 0:001 mm3 a b

i

j

lij …D†

I=I0

u (rad)

1 1 1 2 2 3 3

2 3 4 5 6 5 6

0.5 0.3 0.4 0.2 0.2 0.3 0.3

0.4 1.0 0.8 0.1 0.1 0.2 0.2

1.57b 0 1.24b 0 0 0 0

4 4

5 6

0.4 0.4

0.3 0.3

0 0

Also given are the intensity and mode volume required for numerical evaluation using the Hamiltonian in Eq. (5). Parameters scanned in Fig. 6.

Fig. 6. Target for excitation of a six-level system with two degenerate ®nal states as a function of two of the ®eld phases shown in Table 1.

branches of this Brumer±Shapiro scheme. The computation required about 6 s of CPU time on a 500 MHz Pentium Linux workstation using a Fujitsu-Lahey F95 compiler. The optimization in Fig. 6 is in the usual spirit of molecular coherent control. With the output ®eld now available as a target, one can also ask questions about how sensitive the output ®eld is to the input ®eld or to molecular parameters. Fig. 7 explores two di€erent cases. In Fig. 7A, the initial state was kept the same as before, and the ®eld state j n12

1; u ˆ 1:57ij n13 ; u ˆ 0ij n14 ; u ˆ 1:24i

 j n25

1; u ˆ 0ij n35 ‡ 1; u ˆ 0i Y jn…sp† i  j n45 ; u ˆ 0i

…31†

was chosen as a target, with equal weight given to all other degrees of freedom in the target density

matrix. Eq. (31) corresponds to a three photon pump±pump±dump process. Fig. 7A shows the target as a function of the transition dipole l26 . The di€erent components of the output ®eld are sensitive to the tuning of the transition dipole moment, as exempli®ed by the component shown in Fig. 7A. Unlike simple pump±probe or absorption spectroscopies based on intensities, detection of output ®eld components, which is becoming technically feasible at least for large ®eld amplitudes, promises a much more detailed window onto the quantum system being probed, and the CC2Q approach can be used to analyze such experiments. Fig. 7B shows how the control ®elds can in¯uence spontaneous emission, a calculation not possible by semiclassical means. The same initial state as before was chosen, but this time the target was a density matrix with equal weight for all states having jWsys i ˆ j1i and a photon in the spontaneous emission degrees of freedom corresponding to the 2 ! 1, 3 ! 1 or 4 ! 1 transitions. This spontaneous emission yield would be maximized if the transitions from jWsys i ˆ j2i, j3i, or j4i to jWsys i ˆ j5i or j6i were shut down, allowing the slower spontaneous emission process to return the system from j2i, j3i, or j4i to j1i. As expected, Fig. 7B shows that the maximum target is reached as the photon number n25 approaches zero (with similar behavior for the other modes pumping states 2, 3, 4 ! 5, 6). It is worth reiterating the caveat of Section 2.3 here: with a discretized representation of the emission continuum, the results will be quantitatively correct only if the control ®eld does not allow too many

46

M. Gruebele / Chemical Physics 267 (2001) 33±46

Fig. 7. (A) Target as a function of one of the system transition dipoles. The optimal pulse from Fig. 6 was used, and a three photon pump±pump±dump output component of the ®eld monitored. (B) Target as a function of the intensity of the 2 ! 5 coherent input ®eld. The target density matrix attempts to maximize spontaneous emission from system states 2, 3, 4 back to 1. This occurs when the 2, 3, 4 ! 5, 6 ®elds have been shut down.

photons to enter the spontaneous degrees of freedom. 4. Outlook and applications The full quantum ®eld treatment of coherent control with a time-independent Hamiltonian is conceptually and computationally a simpler problem than time-dependent optimal control. This is mainly due to the form of Eq. (11), which depends only on diagonal elements of density matrices computable by the Hellman±Feynman theorem, and allows enforcement of energy conservation. The size of tractable system±®eld combinations is limited by the need for many ®eld degrees of freedom, but energy conservation as well as propensity and selection rules allow the treatment of interestingly large systems by tiered basis set pruning. The quantum ®eld approach rigorously uni®es time-dependent and independent approaches through the connection between Eqs. (6) and (11). Acknowledgements This work was supported by the National Science Foundation CHE 9986670 and by a UIUC University Scholarship.

References [1] S.A. Rice, M. Zhao, Optical Control of Molecular Dynamics, Wiley, New York, 2000. [2] M. Gruebele, R. Bigwood, Int. Rev. Phys. Chem. 17 (1998) 91. [3] F. Troiani, U. Hohenester, E. Molinari, Phys. Rev. B 62 (2000) R2263. [4] H. Rabitz, W.S. Zhu, Acc. Chem. Res. 33 (2000) 572. [5] M. Gruebele, J. Chem. Phys. 104 (1996) 2453. [6] J.K. Cullum, R.A. Willoughby, Lanczos Algorithms for Large Symmetric Eigenvalue Computations, Birkhauser, Boston, 1985. [7] M. Gruebele, Adv. Chem. Phys. 114 (2000) 193. [8] D.J. Tannor, S.A. Rice, J. Chem. Phys. 83 (1985) 5013. [9] M. Shapiro, P. Brumer, J. Chem. Phys. 84 (1986) 4103. [10] P. Brumer, M. Shapiro, Acc. Chem. Res. 22 (1989) 407. [11] M.N. Kobrak, S.A. Rice, J. Chem. Phys. 109 (1988) 1. [12] D.F. Walls, G.J. Milburn, Quantum Optics, Springer, Berlin, 1994. [13] M. Gruebele, J. Phys. Chem. 100 (1996) 12178. [14] M. Gruebele, J. Phys. Chem. 100 (1996) 12183. [15] M. Gruebele, Proc. Nat. Acad. Sci. USA 95 (1998) 5965. [16] R. Bigwood, M. Gruebele, ACH Models Chem. 134 (1997) 675. [17] R. Chen, H. Guo, Chem. Phys. Lett. 308 (1999) 123. [18] G. Charron, J.T. Carrington, Mol. Phys. 79 (1993) 12. [19] V.A. Mandelshtam, H.S. Taylor, J. Chem. Phys. 107 (1997) 6756. [20] S.W. Huang, T. Carrington, J. Chem. Phys. 112 (2000) 8765. [21] D.J. Tannor, R. Koslo€, S.A. Rice, J. Chem. Phys. 85 (1986) 5805. [22] J.I. Musher, Am. J. Phys. 34 (1966) 276.