EFFECTIVE HAMILTONIAN THEORY AND

MOLECULAR DYNAMICS

Georges Jolicard Observatoire de Besançon. Université de Franche-Comté,

41, Av. de l’Observatoire, B.P. 1615, 25010 Besançon Cedex, France February 14, 1995

1 INTRODUCTION: THE NEED FOR THE

EFFECTIVE HAMILTONIAN APPROACH

IN MOLECULAR DYNAMICS

The theory of effective Hamiltonians was developed within nuclear physics in the sixties by Bloch[1] and Des Cloizeaux[2] as a method for determining the effective interactions between nucleons. Subsequently it was developed and largely used in quantum chemistry[3, 4] to overcome some limitations of both ab-initio and semi-empirical methods[5, 6] A more recent application of the theory is in molecular dynamics, where the basic processes studied are full collisions, in which energy transfers and structural modifications appear between molecular bound states and molecular continua [7], and half-collisions, in which one participating continuum is a photon continuum and the other one is an ionization or dissociation continuum [9]. Most of the recent theoretical treatments of collisions discretize the molecular continua by adding absorbing spatial boundaries and by working with a bounded range of radial coordinates in conjunction with finite basis sets of square integrable functions. Wave packet propagation methods have often been used to integrate these systems using discrete representations[8]. As the required C.P.U. time rapidly becomes prohibitively long when wave packets are propagated on multidimensional grids, new attemps have been made to derive time-independent treatments from the time-dependent ones. A solution within the framework of the Floquet theory has been proposed by Shirley [9] for time-periodic Hamiltonians [10] and has been generalized recently by Peskin and Moiseyev [21] to time non-periodic Hamiltonians.

The discretization approach to molecular continua and to the

time variable appears appropriate for implementation on computers working with discrete digital algebra and it has close links with recent progress and developments in the area of Discrete Variable Representations (D.V.R.). Discretized variable methods require the solution of matrix algebraic problems, in particular the study of complex eigenvalue problems. Good illustrations of such technical problems are provided by the study of purely molecular resonance states, which influence the full collision processes, and of the resonances induced by the field-matter interaction [11, 12, 13] which constitute the basic element for the description of half collisions [14].

From a technical point of view these eigenvalue problems are

characterized by two principal features. The sizes of the vector spaces involved are very large (typically 10^{4} ≤ N ≤ 10^{6}) but the number of eigenvectors effectively contributing to the dynamical processes is often small. This effective confinement of the dynamics to a small subspace is observed in the theory of unimolecular relaxation of polyatomic molecules, when local vibration modes are strongly decoupled from the surrounding bath [15], and also appears more systematically in calculations describing photodissociation experiments. In this last case the dynamical processes are dominated by a few quasi-energy eigenstates of the Floquet matrix [10].

These features prompt the use of projection techniques or ef-

fective Hamiltonian methods to describe the evolution of the vectors which dominate the dynamical processes. The vectors belong to a subspace of relatively small dimension, which is usually called the active space (A.S) The literature offers various formulations for selecting these active spaces and for projecting the quantum dynamics into them. Probably the best known approaches are artificial intelligence techniques [16, 17, 18] for subspace selection and the recursive residue generation method [19], in which transition amplitudes are computed from a recursive calculation which uses the Lanczos algorithm to give the residues and poles of several Green’s functions.

In this volume we will describe both stationary and time-

dependent adaptations of the Bloch wave operator theory [33], paying particular attention to the connections with the (t,t’) theory [21].

2 STATIONARY AND TIME-DEPENDENT AS-

PECTS OF THE BLOCH HAMILTONIAN

THEORY

To aid clarity the theoretical development of this section and the applications presented in the next sections are limited to the field of molecular photodissociation phenomena.

2.1 discretization of space and time

Within the framework of the dipole approximation and using the length gauge, the Hamiltonian of a molecule subjected to a laser pulse is :

H = H_{o}(q, r) + V (q, r) − μ(q, r)E(t)  (1)

H_{o} depends on q, which denotes the bound vibration-rotation coordinate of the isolated fragments, and on the interfragment distance r. It includes the internal energy of the fragments and the kinetic energy of the relative translation. The matter-field coupling involves the dipole moment μ and the electric field amplitude E. Usually it is given as the product of a time-dependent shape function ϵ and a real cosine function with a peak frequency ω.

E(t) = ϵ(t) × cos(ωt)  (2)

Other gauges can be used to express the field-matter interaction [22] but it is not our aim here to analyze their respective merits or defects. The finite representation of the Hamiltonian used in the wave packet propagation algorithm and in the present matrix wave operator theory requires a discretization of the molecular continua. However, the finite basis set of (L^{2}) square integrable functions involved in this discretization introduces artificial periodicities. This drawback has drastic implications for photodissociation calculations since it is not then possible to reproduce the asymptotic behaviour of the resonance states which are generated by the field-matter interaction and which are the major contributor to the dissociation process. These defects can be removed by using similarity transformations on the dissociative radial coordinates in the complex plane [23]. The most simple similarity transformation leads to the conventional complex coordinate method [24]:

H(q, r) ⟶ S^{−1}H(q, r)S = H(q, r exp(iθ))  (3)

More elaborate choices lead to the exterior complex scaling method [24]. Another possible technique is the optical potential model [27]. In this case an imaginary optical potential is introduced asymptotically along the radial axis :

H(q, r) ⟶ H(q, r) − iV_{opt}(r)  (4)

A variational procedure with respect to the parameters included in V_{opt} results in complex plane trajectories which have stationary points close to the exact resonance eigenvalues

In the following it will assumed that the Hamiltonian has been

modified by a similarity transformation or by the addition of an optical potential so as to reproduce correctly the dissociation dynamics using a D.V.R. Nevertheless, for brevity it will still be denoted by the symbol H(q, r).

The discretization of the time variable is another important

problem which has been analyzed recently by Peskin and Moiseyev [21]. The difficulties of constructing the time evolution operator in the case of a timedependent Hamiltonian (eqs 1 and 2) can be partially removed by considering the time as an extra coordinate in a generalized Hilbert space. The time dependent Schrödinger equation.

H(x, t)Ψ(x, t) = ih̄\frac{∂}{∂t}Ψ(x, t)  (5)

where x designates collectively all the molecular coordinate (q, r), can then be replaced by [21]:

Ψ(x, t) = ψ(x, t^{,}, t) |_{t^{,}=t}  (6)

where :

ψ(x, t^{,}, t) = exp^{(−i/h̄)H(x,t^{,})(t−t_{o})} ψ(x, t^{,}, t_{o})  (7)

and H(x, t^{,}) is a Floquet-type operator

H(x, t^{,}) = H(x, t^{,}) − ih̄\frac{∂}{∂t^{,}}  (8)

2.2 The time-dependent wave operator

We consider the time evolution described by (5), under the provisional assumption that, correlated to an active space S_{o} of dimension M, there exists some target subspace S of dimension M in which the dynamics mainly takes place. M is supposed small compared to the dimension (N − M) of the orthogonal complement subspace S_{o}^{+}, and the projection operators of the active space and its orthogonal complement are denoted by P_{o} and Q_{o},respectively. A simplified description of the relaxation issuing from space S_{o} can be derived by separating the evolution localized inside the subspace S_{o} and the transitions into the complementary space S_{o}^{+}. We associate with the Hamiltonian H and with the active space S_{o} [28] a time-dependent wave operator Ω with the defining equation

Ω(t, 0) = U(t, 0; H) × (P_{o}U(t, 0; H)P_{o})^{−1}.  (9)

By using the matrix multiplication properties of partitioned matrices, represented by the projection operator properties P_{o}P_{o} = P_{o}, P_{o}Q_{o} = O, we obtain the result :

Ω(t, 0) = P_{o} + X(t, 0)  (10)

The operator X(≡ Q_{o}ΩP_{o} ≡ Q_{o}UP_{o}(P_{o}UP_{o})^{−1}), called the reduced wave operator, is an off-diagonal part of Ω; in the N × N matrix representation it has non-zero elements only in the lower left (N − M) × Mblock. The introduction of this scheme into the equations describing the evolution issuing from space S_{o} leads to :

U(t, 0; H)P_{o} = (P_{o} + X(t, 0))(P_{o}U(t, 0; H)P_{o})  (11)

which exhibits an evolution divided into two parts. The first part is represented by P_{o}UP_{o} and is inside the space S_{o}. The second part is described by the operator X and represents the transition from the space S_{o} to the complementary space S_{o}^{+} without any reverse paths ( a reverse path would involve P_{o}UQ_{o}).

Equation (11) performs an important simplifying factoriza-

tion, since the global evolution described by the operator U(t, 0; H)P_{o} is, by contrast, an infinite series of back and forth transitions.At this stage of the discussion, equation (11) is still purely formal insofar as the operators ω and X are expressed via the evolution operator, even when the aim of this theory is to regard the projected term P_{o}UP_{o} as a time development operator associated with an effective Hamiltonian H_{eff}. This connection with effective Hamiltonian theory has been set out in ref [28]. A pedagogical presentation can also be found in ref [29]. These developments leads to the equation

U(t, t_{o}; H)P_{o} = (1 + X(t, t_{o}))U(t, t_{o}; H_{eff})  (12)

where

H_{eff}(t) = P_{o}H(t)(P_{o} + X(t, t_{o}))  (13)

and where the associated equation of motion for the reduced wave operator X is

ih̄\frac{∂X(t, t_{o})}{∂t} = Q_{o}(1 − X(t, t_{o}))H(t)(1 + X(t, t_{o}))P_{o}  (14)

The price to be paid for the reduction of the dimension of the space is the requirement to solve a non-linear differential equation.

2.3 The Bloch wave operator as the adiabatic limit of the

time-dependent wave operator

This theory and such subsequent developments as intermediate Hamiltonian theory [30] were developed in Quantum Chemistry to obtain new Effective Hamiltonian which can provide accurate molecular energies and easily interpretable wavefunctions while avoiding difficulties (e.g. intruder state effects) which appear in perturbation formalisms. The basic idea is to define an effective Hamiltonian H_{eff}, acting in a reduced model space S_{o}, such that a subset of the exact energy levels of the Hamiltonian H coincides with the energy levels of H_{eff}. The non-unitary transformation which changes the Hamiltonian into an effective Hamiltonian is the wave operator. The choice made by Bloch corresponds to setting the wave operator equal to

Ω = P (P_{o}P P_{o})^{−1}  (15)

The theory of Bloch in its initial version is thus a purely stationary state method; however, it can also be derived as a limiting case of the time-dependent wave operator theory. Let us assume that the operator H is generated from the unperturbed operator H_{o} by adiabatically introducing the perturbation V (≡ H − H_{o}) :

H(t) = lim_{γ→o^{+}}(H_{o} + V exp(γt))
H_{o} = H(t = −∞)
H = H(t = 0)  (16)

it can then be proved that the time derivate ∂X/∂t tends to zero at any time [33], so that eq. (14) takes the form

Q_{o}(1 − X)H(1 + X)P_{o} = 0  (17)

This equation is the basic equation of Bloch’s wave operator theory, although it is presented in the literature in a different but equivalent form

HΩ = ΩHΩ  (18)

Equation (17) has the advantage over eq. (18) of showing clearly that the wave operator is a non-unitary and non-singular transformation (T = 1 + X) which possesses a trivial inverse (T^{−1} = 1 − X)[31] (see fig.1). This transformation cancels the couplings which connect S_{o} to S_{o}^{+}. The procedure followed to derive the eigenvalues originating within the space S_{o} and the associated eigenvectors is unexpectedly simple. The eigenvalues are obtained by diagonalising H_{eff} = P_{o}HΩ inside the space S_{o} :

(λ^{o})^{−1}(P_{o}HΩ)λ^{o} = P_{o}EP_{o}  (19)

where E denotes the diagonal eigenvalue matrix. The eigenvectors | λ_{i}^{o}⟩ of H_{eff}(eq. 19) are the projections into the model space of the corresponding exact Figure 1: Schematic representation of the effective Hamiltonian calculation by the wave operator transformation

eigenvectors | λ_{i}⟩. The wave operator Ω tranform these projected eigenvectors back into the exact eigenvectors in the full space

| λ_{i}^{o}⟩ = P_{o} | λ_{i}⟩
| λ_{i}⟩ = Ω | λ_{i}^{o}⟩  (20)

These developments show that the wave operator Ω(or its non-diagonal part X) plays a central role. It participates in the construction of H_{eff} and connects the eigenvectors of H_{eff} with the exact eigenvectors (eq. 20). Thus solving the eigenvalue problem requires us to calculate this operator and consequently to solve the non-linear equation (17). For this task several linear and quadratic methods have been proposed in the literature [3]. However, some of them, while seeming simple and attractive, introduce drastic implicit and unjustifiable approximations, such as assuming a diagonal structure of H in the large complementary space S_{o}^{+}. An iterative solution technique called the Single Cycle Method (S.C.M.) has recently been proposed in conjunction with the Gauss-Seidel method [32]. This method expands the transformation 1 + X as a series of partial transformations

1 + X = ∏_{N} ∏_{i∈S_{o}} ∏_{α∈S_{o}^{+}}(1+ | α⟩⟨i | ΔX_{α,i}^{N})  (21)

The three indices N, α and i are incremented successively. For example the row index α of the reduced wave operator in the outer space S_{o}^{+} is incremented first, then the column index i and finally the index N which denotes the number of cycles of the iteration. After each partial transformation the H elements are modified before being used for the next transformation. Namely, by introducing a generalized iteration index (L ≡ (N, α, i)) incremented after each partial transformation :

H^{L} = (1 − ΔX^{L})H^{L−1}(1 + ΔX^{L})  (22)

with

ΔX^{L}(≡ ΔX_{α,i}^{N}) = \frac{⟨α | H^{L−1} | i⟩}{⟨i | H^{L−1} | i⟩ − ⟨α | H^{L−1} | α⟩}  (23)

The iteration, initiated with H^{L=0}, is stopped when ΔX^{L} is negligible. A more simple iteration, which is less efficient from the point of view of the convergence properties, is obtained by modifying the Hamiltonian after each complete cycle N → N + 1 and not after each partial transformation. This method known as the Recursive Distorted Wave Approximation (RDWA) [33] involves a direct relationship (H^{L} = f(H)) instead of a step by step relationship (H^{L} = f(H^{L−1})).

2.4 The Quasiadiabatic evolution of the instantaneous eigen-

vectors

The principal aim of the wave operator theory is to reduce the dimension of the vector spaces directly participating in the dynamics by introducing the concept of an effective Hamiltonian and deriving corresponding projected equations of motion. Altough the time-dependent version is devoted to time-dependent processes and the stationary version to general eigenproblems, the two versions of this theory are not independent. As we have indicated previously, the Bloch wave operator is the adiabatic limit of the time-dependent one, and promising relationships between the time-dependent and time-independent formulations can be derived from this feature.

In a more pragmatic fashion the two theories can be used to-

gether in stationary problems to overcome convergence difficulties within iterative solution schemes in the strong coupling regime. This possibility arises from the fact that the iterative scheme (21-23) accepts any arbitrary trial value and not only ΔX^{L=0} = 0. The choice of a trial operator near to the converged solution increases the speed of convergence and can even make troublesome divergences disappear. Such an optimal trial operator can be derived from a time-dependent procedure by etablishing the couplings almost adiabatically, i.e.

H(t) = H_{o} + V F(t), −T ≤ t ≤ 0
F(−T) = 0, F(0) = 1  (24)

If the switching time T tends to infinity and if simultaneously the time derivative ∂F/∂t tends to zero at any time, the propagated wave operator becomes identical to the stationary one. When a finite switching time T is used, the corresponding wave operator, which exhibits some deviations from the purely stationary case, can still be used as a trial operator in the stationary iterative procedure. The time propagation can be accomplished using eq.(14) together with second order differencing to approximate the time derivative and to preserve time reversal symmetry.

X(t + dt) = X(t − dt) + 2dt\frac{∂X(t)}{∂t}  (25)

with :

\frac{∂X(t)}{∂t} = \frac{1}{ih̄}Q_{o}(1 − X(t))H(t)(1 + X(t))P_{o}  (26)

The central question concerning this double procedure is: what is the advantage of using the non-linear equation (14) and deriving the quasi-adiabatic solution in place of using the linear Schrödinger equation directly? The main advantage is the small magnitude of the time derivative in eq.(26). Indeed, the wave function in the adiabatic limit, T → ∞, is a product of two terms.

Ψ(t) =| λ_{i}(t)⟩ exp(∫_{−∞}^{t} −iE_{i}(t^{,})/h̄dt^{,})  (27)

where | λ_{i}(t)⟩ is an instantaneous eigenvector, the components of which are slowly varying time functions, and the second factor is a purely phase term involving the integral of the instantaneous eigenvalues. It is this second factor which produces the largest contribution to the time derivative ∂Ψ/∂t present in the Schrödinger equation. If one now considers the one-dimensional space S_{o} spanned by the vector (| i⟩ ≡ | λ_{i}(−∞)⟩) and forms the corresponding nondegenerate wave operator, one obtains:

Ω(t) = ∑_{α} | α⟩⟨i | (⟨α | Ψ(t)⟩/⟨i | Ψ(t)⟩)
= ∑_{α} | α⟩⟨i | (⟨α | λ_{i}(t)⟩/⟨i | λ_{i}(t)⟩)  (28)

In this last expression the quickly varying phase term appearing in eq.(27) has disappeared and the time integration is consequently easier. This feature persists when degenerate active spaces are used. Moreover these degenerate active spaces can be introduced in a way which suppresses the large couplings which appear near crossings and avoided crossings in the quasiadiabatic propagation; the states participating in these avoided crossings must simply be included in the active space S_{o} in order to achieve this.

A simple illustration is presented in figs. 2 and 3. This corre-

sponds to a system of two harmonic oscillators in near resonance, with linear and quadratic couplings. The Hamiltonian in creation and annihilation operator form is

H = (a_{1}^{+}a_{1}+ \frac{1}{2})hν_{1}+(a_{2}^{+}a_{2}+ \frac{1}{2})hν_{2}+λ[(a_{1}^{+}+a_{1})+(a_{2}^{+}+a_{2})+(a_{1}^{+}+a_{1})(a_{2}^{+}+a_{2})]

(29)

with hν_{1} = 1., hν_{2} = 1.05 and 0 ≤ λ ≤ 0.6. The lower part of the spectrum is presented in Fig.2. A large part of this spectrum has been reproduced successfully with the quasiadiabatic procedure. For this, the coupling para-Figure 2: Lower part of the spectrum of two coupled harmonic oscillator versus the perturbation parameter λ (see eq.29)

Figure 3: (a): The nonadiabatic factor (eq.30) for the four states (i = 1 to 4)(∘, +, ×, ∗) of the quasidegenerate subspace with 3 quanta. The factors are shown versus the index j of the harmonic basis vectors sorted in the order of increasing energy values. The switching time of the perturbation is equal to 3 vibrational periods. (b):In this case the switching time of the perturbation is equal to 30 vibrational periods

meter λ is represented as a linearly increasing function of time on the interval [t = 0, t = N/ν_{1}]. Then a quasidegenerate active space is selected (in the present case a collection of the unperturbed states corresponding to a given number of vibrational quanta) and the propagation equations (25-26) are integrated. The trial wave operator thus obtained is subsequently introduced as initial,i.e. zeroth order, element in the iterative SCM procedure (eqs(21-23)). Finally the eigenvalues are derived by diagonalizing the effective Hamiltonian (eq 19), and the wave operator is applied to these projected eigenvectors to transform them back into the exact eigenvectors (eq.20). An interesting point is the analysis of the errors produced by using a finite switching time interval rather than taking the purely adiabatic case. The resulting deviations are represented in fig.3 by the factors

log | ⟨λ_{j}^{exact} | λ_{i}^{app}⟩ − δ_{i,j} |  (30)

for the quasidegenerate subspace spanned by the 4 vectors with 3 quanta : (3, 0; 2, 1; 1, 2; 0, 3) and for a coupling parameter λ = 0.4 (eq.29). The switching time interval corresponds to three vibrational oscillations (N = 3) in fig.3_{a} and to thirty oscillations in fig.3_{b}. As expected one obtains a decrease of the deviations by increasing the switching time; nevertheless it is surprising to note the relatively small deviation obtained with the short switching time of three vibrational periods.

2.5 The connection with (t, t^{,}) theory

The wave operator theory is an efficient tool for investigating the formal relationships between the stationary and the time-dependent treatments of molecular processes. Its careful analysis allows us to connect the present theory to the matrix Floquet theory for time-periodic Hamiltonians and also to the (t, t^{,}) theory, which is in a sense the generalization of the Floquet theory to non-periodic Hamiltonians

By introducing the Floquet-type operator H which is associated with the Hamiltonian H and which is defined as

H(t) = H(t) − ih̄∂/∂t  (31)

one can transform the basic equation (14) of the time-dependent wave operator theory. After some algebric manipulations, one obtains [34] :

Q_{o}(1 − X(t, t_{o}))H(t)(1 + X(t, t_{o}))P_{o} = 0  (32)

Alternatively, combining eqs.(32) and (13) yields the result

H(t)Ω(t, t_{o}) = Ω(t, t_{o})H_{eff}(t)
= Ω(t, t_{o})H(t)Ω(t, t_{o})  (33)

The most prominent feature is the close resemblance between the pair of equations (32) and (33) and the equations (17) and (18) of the stationary theory; the Floquet type Hamiltonian H must simply be taken in place of H to generalize the stationary effective Hamiltonian theory to the time-dependent one. When the active space S_{o} is a one-dimensional space, eq.(33) becomes an eigenvalue equation : the wave operator is then an instantaneous eigenvector of the Floquet-type Hamiltonian H(t) and H_{eff}(t) is the corresponding instantaneous eigenvalue. To take advantage of the similarity between the stationary and the time-dependent formulations it is necessary to give a new status to the time coordinate, as is done in (t, t^{,}) theory, by introducing a new time variable t^{,}, acting as a time coordinate in a generalized Hilbert space [35]. In this generalized Hilbert space, ”box” normalization is used for x and for t^{,} (0 ≤ t^{,} ≤ T) and the inner product of two functions is defined as :

⟨⟨ϕ_{i} | ϕ_{j}⟩⟩ = \frac{1}{T} ∫_{0}^{T} dt^{,} ∫_{−∞}^{∞} dxϕ_{i}^{∗}(x, t^{,})ϕ_{j}(x, t^{,})  (34)

T is the time period in the case of purely time-periodic Hamiltonians; it can also be identified with the finite duration of the field-matter interaction for non-periodic cases. In this new framework any instantanenous eigenvalue can be determined from the 1 × 1 effective Hamiltonian (eq.33) by expanding the initial wave function (here assumed to be identical to a basis vector ϕ_{1} of the generalized Hilbert space) in term of a complete set of quasienergy states i.e stationary solutions of the Floquet Hamiltonian H :

| Ψ(t_{o})⟩ = ∑_{k} c_{k} | λ_{k}(x, t^{,})⟩ |_{t^{,}=t_{o}}
H | λ_{k}⟩ = E_{k} | λ_{k}⟩  (35)

By introducing the wave operator corresponding to this one dimensional active space, Ω ≡| ψ⟩/⟨ϕ_{1} | ψ⟩ and by using the (t, t^{,}) theory to propagate the wave function ψ, one obtain the result

Ω(t, t_{o}) = ∑_{j} | ϕ_{j}⟩⟨ϕ_{1} | (\frac{∑_{k} c_{k} exp(−iE_{k}(t − t_{o})/h̄)⟨ϕ_{j} | λ_{k}(x, t^{,})⟩}{∑_{k} | c_{k} |^{2} exp(−iE_{k}(t − t_{o})/h̄)})_{t^{,}=t}
E_{eff}(t) = \frac{∑_{k} E_{k} | c_{k} |^{2} exp(−iE_{k}(t − t_{o})/h̄)}{∑_{k} | c_{k} |^{2} exp(−iE_{k}(t − t_{o})/h̄)}  (36)

This equation reveals that the effective eigenvalue is a complex eigenvalue. If the expansion (35) is reduced to one term k_{1} the effective eigenvalue becomes a quasienergy and the column which represents the one-dimensional wave operator becomes identical to the corresponding quasivector | λ_{k_{1}}⟩.

Equation (36) has the pedagogical value of revealing the origin

of the time-dependent effective eigenvalue but is not of much pratical use. In the search for solutions of the time-dependent Schrödinger equation it is certainly more efficient to represent the Floquet operator by using the set of Fourier basis functions: (exp(2πint^{,}/T) ; n = −∞, …, +∞). In most of the photodissociation experiments the initial state is a t^{,}-independent vector function ψ(x, t_{o}). Consequently the time evolution operator has the simple form [9]

U(t, t_{o}) = ∑_{n=−∞}^{∞} exp^{2πint/T} [exp^{−iH(x)(t−t_{o})/h̄}]_{n,0}  (37)

where 0 and n refer to the Fourier basis functions.

The evolution equation (37) can be further simplified when de-

scribing most photodissociation experiments, since the photodissociation processes, while apparently complicated, can often be described using a small target space formed by only a few quasi-vectors. Several applications which illustrate this feature will be represented in the next three sections. We must now choose how to define the active space which is expected to reproduce correctly the dynamic features of the relaxation process and also how to project the dynamics into this active space.

Choosing a model space for a static eigenvalue problem is usu-

ally not too difficult, but the situation is quite different in the present case. The model space is not a fixed space into which the reduced quantum information is projected, but rather an evolving active space within which the dynamics is confined; this makes its definition more intricate. However, our previous demonstration that stationary and time-dependent processes can be closely related through the time-dependent wave operator concept and the (t, t^{,}) theory enables us to adopt a systematic selection procedure for an active space

Let | i⟩ ≡| n = 0, α⟩ designate the initial state of the matter-

field space and Let S_{o} be the trial one-dimensional active space spanned by this vector. The reduced stationary wave operator corresponding to this space and to the Floquet type Hamiltonian is constructed. To do this the S.C.M. (eqs.21-23) is iterated up to a high order N_{iter} and vectors of the whole space are sorted in the order of decreasing values of the matrix component moduli of the wave operator

| ⟨γ | X | i⟩ |  (38)

with the initial state | i⟩ being arbitrarily placed in first position. With a sufficiently large iteration order N_{iter}, all the states of the whole space are ordered and the active space, the dimension N_{o} of which is imposed, can then be formed by taking the N_{o} first vectors. At first sight it might seem inappropriate to use a stationary formulation to select an active dynamical space but the procedure appears more reasonable when one recalls that the stationary wave operator associated to the Floquet Hamiltonian is identical to the time-dependent wave operator (eqs 32-33).

In this partitioning of the vector space, the relaxation produced

by the term exp(−iH(t − t_{o})/h̄) in eq.(37) can be separated into two parallel and independent paths by applying a nonunitary transformation T based on the two reduced wave operators X_{P} and X_{Q} [4]. This transformation supresses the couplings generated by H between the two subspaces S_{o} and S_{o}^{+}. We have

T = 1 + X_{P} + X_{Q}  (39)

with

T^{−1} = P_{o}P + Q_{o}Q  (40)

The usual notations P_{o}, Q_{o} and P, Q for the projection operators of S_{o} and S_{o}^{+} and of the associated target spaces S and S^{+} are used.

The projection procedure is achieved by neglected the relax-

ation path through the space S_{o}^{+}. The weakness of the contribution from this path results from the active space selection rule which has been adopted and which amalgamates in S_{o}^{+} isolated states which possess a weak overlap with the initial state, groups of states of high density whose collective influence disappears through interference effects or, more simply, states with small lifetimes and states which are energetically forbidden. This leads to the result

exp(−iHΔt/h̄)P_{o} = ΛP_{o} exp[−i(P_{o}EP_{o})Δt/h̄](P_{o}ΛP_{o})^{†}  (41)

where the quasivectors projected into S_{o}, i.e. P_{o}ΛP_{o}, also form the (N_{o} × N_{o}) matrix which diagnalizes the effective Hamiltonian P_{o}HΩ_{P}. The few columns ΛP_{o} are reconstituted from the projected quasivectors using eq.(20), namely

ΛP_{o} = Ω_{P} (P_{o}ΛP_{o})  (42)

3 SELECTION OF AN ACTIVE SPACE

The efficiency of the wave operator formulation in choosing active spaces and its capicity to correctly reproduce the relaxation dynamics by working in small vector subspaces can be tested by looking at the time evolution of some simple closed systems which involve the interaction of a field continuum with bound molecular states and of some open systems which involve the interaction of an initial field continuum and a final molecular continuum.

3.1 Rigid rotator in a strong laser field

We take as an exemple of a closed system the energy relaxation of a rigid rotator with a permanent dipole moment μ under the influence of a coherent monochromatic electric field with frequency ω and polarized along the z direction. With (J, m)denoting the rotator quantum numbers, the corresponding Floquet matrix, diagonal with respect to the projection quantum number m has the elements

⟨J, m; n | H | J^{,}, m; n^{,}⟩ = [Bh̄J(J + 1) + (n + \frac{1}{2})h̄ω]δ_{JJ^{,}}δ_{nn^{,}}
− \frac{μE}{2} (δ_{n,n^{,}+1} + δ_{n,n^{,}−1})⟨J, m | cos θ | J^{,}, m⟩(43)

where B is the rotation constant, E the amplitude of the electric field and θ the angle between the rotator dipolar moment axis and the fixed orientation of the field. The case of the HF rotator (with Bh̄ = 21cm^{−1}, μ = 1.91 Debye) in a monochromatic field of intensity 1T W, with h̄ω = 1500cm^{−1}, has been analyzed using a basis of 220 states consisting of the direct product of the rotational basis (| J, m = 0⟩, 0 ≤ J ≤ 19) with the 11 Fourier states (−5 ≤ n ≤ 5). The intensity of the field in this case is large enough to produce non perturbative effects and so to provide a severe test for the wave operator subspace selection procedure. The quantity which is of experimental interest is the rotational transition probability averaged over the initial time t_{o}. Equations (37) and 41) yield the result

P_{J→J^{,}}(Δt) = ∑_{n} | ∑_{η∈S_{o}}⟨J^{,}, 0; n | λ_{η}⟩⟨λ_{η} | J, 0; 0⟩ exp(−iE_{η}Δt/h̄) |^{2} (44)

The wave operator formulation applied here plays a double role in the projection of the Schrödinger equation into reduced subspaces of fixed size. First it selects the active space which is best adapted to describe the exact relaxation process. Second, it constructs the target space and the corresponding effective Hamiltonian which are correlated to this active space (a more coarse grained approximation would consist of using the primary effective Hamiltonian simply projected into the active space).

Fig.4_{a} presents the results for the calculated survival proba-

bility J = 0 → J = 0 using two active spaces of increasing dimension N_{o} = 11, N_{o} = 17. The results as compared to an exact integration are good, in view of both the small dimension of the model spaces (11 and 17 ≪ 220) and the high intensity of the laser field. Fig.4_{b} illustrates the advantage of correctly selecting the active space dynamically using the wave operator procedures, as opposed to using a fixed basis covering a band of energies around the initial energy.

Figure 4: (a): The transition probability P (J = 0 → J = 0) as a function of the time, measured in laser fields periods. Continuous line: exact integration. Dotted and dashed lines: results using active spaces, of size N_{o} = 11 and N_{o} = 17 respectively, as selected by the wave operator algorithm. The wave operator results often coincide with the exact results, making it difficult to distinguish them. (b): Same as fig.4_{a} but using active spaces of dimension N_{o} = 11 and N_{o} = 39 with basis sets forming energy windows centered on the initial state. The results in these static active spaces diverge from the exact result after a few optical cycles.

3.2 Photodissociation of a diatomic molecule

An example of active space selection for an open system is presented in fig.5. The process considered is the direct photodissociation of the HF molecule, the molecular potential used being the Stone and Noid Morse potential [36]. This potential surface is drawn in figs.5 in a dressed picture involving 10 Floquet blocks (−6 ≤ n ≤ 3). The surface n = 0 corresponds to the initial state, while the surfaces numbered n and −n are the dressed surfaces associated with the absorption and emission of n photons. Each Floquet block is represented using a basis set of 150 field-free eigenstates, composed of 6 bound states and 144 complex eigenvectors which discretize the continuum. Fig.5, also shows the states selected by the wave operator algorithm when the dimension N_{o} = 12 of the active space is imposed, the initial state being the ground state (v = 0, n = 0). The figures reveal that only the adjacent blocks (n = 1, −1, −2) are involved in the reduced subspace and that the selection procedure gives priority to various bound states in the adjacent blocks and to the diffusion states of the final Floquet block which are in near resonance with the initial state. An analysis of these figures gives us a rough understanding of the photodissociation mechanism. With a laser frequency equal to 0.60 a.u. the dissociation proceeds via a direct one-photon process from the ground state directly to the continuum. When the laser is tuned to the bound-bound transition (v = 0 → v = 2), i. e. ω = 0.02914, this is, by contrast, a sequential two-photon process which involves the intermediate state | n = −1, v = 2⟩ in resonance with the inital state. A Figure 5: (a):Eigenvectors selected by the wave operator algorithm to describe the photodissociation process for the HF molecule. The field intensity is I = 1T W/cm^{2} and the frequency ω = 0.060a.u.. The dimension of the active space N_{o} = 12 is imposed. (b) The same as fig.5_{a} but for a laser frequency equal to ω = 0.02914 detailed analysis of the selected eigenvectors reveals that the increase of the dissociation probability at ω = 0.02914 is produced by a mixing of the two resonances | λ_{n=0,v=0}⟩ and | λ_{n=−1,v=2}⟩ giving rise to a significant decrease in the lifetime of the former.

4 THE RESONANCE STATES OF VAN DER

WAALS MOLECULES

The resonance state problem plays a central role in describing the vibrational energy redistribution and photofragmentation processes of van der Waals (v.d.W.) molecules. The calculation of resonance states is difficult to handle and many indirect treatments have been proposed in the literature. Resonance state information can be derived from inelastic transition probabilities in the study of the collision of the fragments [37], while resonance energies and widths are extracted from correlation functions in the Prony method [38]

However, the resonance state problem is essentially an eigen-

value problem, and so the wave operator formulation provides a third approach which is particularly suitable. The main advantage of this direct technique is that it gives the eigenvector together with the eigenvalue and thus provides access to the branching ratio by means of an asymptotic analysis of the outgoing fluxes in the various open channels. The relationships between the quantum fluxes and the partial widths have been analyzed in refs. [39, 40]. Defining by Γ_{α}^{(n)} the partial width of resonance state α on channel n and taking into account the common exponential time dependence of the flow produced by the resonance on each channel, one obtains the result

Γ_{α}^{(n)} = Γ_{α} \frac{F_{n}}{∑_{n} F_{n}}
= Γ_{α} \frac{∂F_{n}/∂t}{∑_{n} ∂F_{n}/∂t}  (45)

where the instantaneous flux ∂F_{n}/∂t possesses a simple expression in terms of the components of the resonance state when a D.V.R. is used [40].

We present here some results relating to the vibrational pre-

dissociation of a triatomic molecule XI_{2}(B^{3}Π), X = Ne, Ar. The study used a colinear model, involving two coupled Morse oscillators, by taking as coordinates the internal coordinate of the I_{2} diatom and the distance between X and the centre of mass of I_{2}. Our treatment assumes that the resonance states of the XI_{2} complexes arise from the uncoupled eigenstate products ϕ_{v}(R_{I_{2}})×ϕ_{l}(R_{XI_{2}}). This assumption allows us to identify the resonance states by using the same pair of indices (v, l). which are used for the free states.

Fig.6 shows the behaviours of the local partial widths Γ^{n}, n =

v − 1 and n = v − 2 of the states ϕ_{v=2,l=4}^{res} as a function of the dissociative coordinate R_{XI_{2}}. These curves reveal that NeI_{2} and ArI_{2} dissociate essentially by transferring one vibrational quantum from I_{2}. The variation of the partial widths with respect to the dissociate coordinates is an indication of the lack of completness of the basis set used. The less satisfactory stability in the case Figure 6: (a): The local partial width v − 2 (solid line) and v − 1 (dashed line) of the (v = 2, l = 4) resonance state of NeI_{2} vs the dissociative coordinate R_{XBC}. (b): Same as case (a) but for ArI_{2}.

of ArI_{2} results from a larger reduced mass, which corresponds to shorter de Broglie wavelenths and so to the requirement for much finer spatial grids.

5 PHOTODISSOCIATION PRODUCED BY SHORT

LASER PULSES

Most of the studies of the photodissociation processes induced by short and intense lasers pulses use wave packet formulations. A current interest with experimental relevance is the search to understand how the pulse shape can influence the photodissociation process and, more precisely, how it can close or open some dissociation channels and modify the branching ratios. This search involves comparing the results of many numerical integrations, all corresponding Figure 7: Photodissociation of H_{2}^{+}: proton relative population versus kinetic energy. The solid curves correspond to the channel 1sσ_{g}. The three spectra of increasing width correspond to electric field gaussian envelopes of decreasing lifetime: τ =200 a.u., 100 a.u., 50 a.u. The dotted lines correspond to the channel 2pσ_{u}.

to the same initial molecular state and the same field amplitude but using various pulse shapes. The standard wave packet methods then require large C.P.U. time, even when the use of DVR basis sets and Fast Fourier Transforms (F.F.T.)accelerate the calculations. The wave operator theory can appreciably reduce this time and also show the process in a new light, if one introduce a discretized model of the electric field and the concept of an evolutional target space. If the electric field time dependence is written E(t) = F(t) cos(ω_{o}t) then the discretization model divides the whole range of variation of the envelope function between F = 0 and F^{max} into N elementary intervals of the same width ΔF = F^{max}/N. The horizontal line F_{i} = i × ΔF cuts out the continous function F(t) in an even number of points t_{i}. The continous function F(t) is then approximated by the discontinuous stepped function F (t) which takes the constant value F(t_{i}) when t is between the two adjacent point t_{i} and t_{i+1}. In this discretized scheme the evolutional target space model is set up by first defining an active space, S_{o} selected by using the wave operator selection procedure. For each discretized value of the field magnitude, a target space is associated with the active space S_{o} → S^{(j)}. The vectors which span S^{(j)} are the Floquet eigenvectors | λ_{η}⟩ derived from the relevant free vectors | η⟩ by adiabatically switching on the field-matter interaction.

Within this discretized model the field-matter interaction is a

series of discontinuous constant couplings. The propagation produced by these constant couplings is trivially calculated; it corresponds to diagonal effective Hamiltonians inside the successive target spaces. Finally, it is assumed that at the discontinuity between two adjacent field values, namely j and j + 1, the wave function located in the subspace S^{(j)} is completely projected into the subspace S^{(j+1)}, i.e. the projection outside this subspace is regarded as negligible. The wave function obtained at time t can then be written

Ψ(t) = ∑_{γ}^{,} | λ_{γ_{n}}⟩ exp(\frac{−iE_{γ}^{n}Δt_{n}}{h̄}) ∑_{β}^{,}⟨λ_{γ_{n}} | λ_{β_{n−1}}⟩
× exp(\frac{−iE_{β}^{n−1}Δt_{n−1}}{h̄}) … exp(\frac{−iE_{μ}^{1}Δt_{1}}{h̄})

× ⟨λ_{μ_{1}} | 0, α⟩  (46)

where| 0, α⟩ denotes the initial molecular state and the ∑^{,} summations are limited to the target spaces.

The interest of this expression is that it is dependent on the

pulse shape only via the successive time intervals Δt_{n}. Thus, the investigation of various pulse shapes does not require us to repeat the full dynamical integration if some data has been stored during the construction of the target spaces; namely for each discretized value of the field:

1. the eigenvalues of the Floquet states spanning the target spaces;

2. the overlap matrices of Floquet eigenvectors for consecutive target spaces.

Figs.7 and 8 show the results of applying this model to calculate the photodissociation of the (v = 0, j = 1) rovibrational level of the H_{2}^{+} ground state subjected to short electromagnetic pulses. Fig.7 corresponds to gaussian field pulses of increasing width. The calculations were made by using the A.P gauge, invoking the dipole approximation and neglecting the A.A term. The maximum of the laser intensity is I_{max} = 0.4 10^{13}W/cm^{2} and the laser wavelength λ = 329, 7nm is chosen to introduce an energetically favorable vertical three-photon transition from the bound (1sσ_{g}, v = 0, j = 1) to the continuum (2pσ_{u}) state. The figure displays the proton relative populations as a function of kinetic energy. This spectrum is obtained by performing a half-time Fourier transform of the asymptotic amplitude of the outgoing wave function in the two exit channels [41] .It Figure 8: (a): Photodissociation of H_{2}^{+} induced by the absorption of one photon of wavelenght λ = 113.7nm. Lower panel: the electric field envelope versus time in a.u. Upper panel: the dissociation probability (solid line) and the difference between 1 and the squared modulus of the wavefunction (dashed-dotted line). (b): Proton relative population versus kinetic energy (in a.u.)

should be noted that the widths of these spectra are in perfect agreement with the theoretical predictions of ref.[42], thus indicating a mechanism for the twophoton peaks which involves a three-photon absorption, followed by an emission of one photon which occurs during the dissociation process itself. In fig.8 the maximum of the laser intensity is I_{max} = 0.62 10^{13}W/cm^{2} and the wavelength λ = 113.7nm induces a direct one-photon transition to the continuum. This figure reveals that by using a double gaussian pulse, it is possible to give an oscillating structure to the kinetic spectrum and thus to make some values of the velocities of separation of the fragments disappear from the spectrum.

6 CONCLUSION

Most recent treatments of dynamical problems involving coupled molecular and photon continua require discrete representations using bases of large dimension. The wave operator formulation proposes a new strategy to reduce the dimension of the dynamical vector space needed without losing the essential details of these processes. Various applications, represented here by photodissociation processes and in Wyatt et al. by the overtone relaxation of benzene [43] indicate that this formulation is one of the most efficient to select the relevant active spaces. It does not yet have such a well-established reputation as a tool to project the dynamics into the active space and to construct the effective Hamiltonian, but we have high hopes of its future utility in this area. Promising lines for further investigation include the study of the close links with the (t, t^{,}) theory of Peskin and Moiseyev (sketched in section 2.1) and the exploitation of the links between the time-dependent and the time-independent formalisms; for example, it might be possible to study the quasiadiabatic evolution of the time-dependent wave operator as a first step in the construction of a reasonable trial stationary wave operator which can then be further refined in an iterative procedure.

Acknowledgement: The author would like to thank John P. Killingbeck for a careful reading and correcting of the manuscript and for valuable comments.

References

[1] Bloch, C. 1958. Nucl. Phys. 6: 329

[2] Des Cloizeaux, J. 1960. Nucl. Phys. 20: 321

[3] Barthelat, J. C., Durand, Ph. 1992. In Computational Chemistry. Structures-interactions and Reactivity. Part A: p. 26. ed. S. Fraga. Amsterdam: Elsevier

[4] Durand, Ph. 1985. J. of Mol. Struc. 120: 443; Durand, Ph., Malrieu, J. P. 1987 Ab Initio Methods in Quantum Chemistry I ed. K.P. Lawley. New York: J. Wiley and Sons Ltd.

[5] Salem, L. 1966 The molecular orbital theory of conjugated systems. New York: Benjamin

[6] Freed, K. F. 1977 Semi empirical methods of electronic structure calculation, Part A: Techniques ed. G.A. Segal. New York: Plenum Press

[7] Levine, R. D. 1969 Quantum mechanics of molecular rate processes. Oxford: Clarendon Press

[8] Kosloff, R. 1988. J. Phys. Chem. 92: 2087

[9] Shirley, J. H. 1965. Phys. Rev. B. 138: 979

[10] Chu, S. I. 1989. Adv. Chem. Phys. 73: 739 [11] Junker, B. R. 1982. Advan. At. Mol. Phys. 18: 207

[12] Lefebvre, R. 1984. J. Chem. Phys. 88: 4839

[13] Kukulin, V. I., Krasnopol’sky, V. M., Horacek, J. 1989. Theory of resonances. Dordrechr: Kluwer Academic Publishers

[14] Reinhardt, W. P. 1982. Ann. Rev. Phys. Chem. 33: 223

[15] Dübal, H. R., Quack, M. 1984. J. Chem. Phys. 81: 3779; ibid 72: 342; Peyerimhoff, S. D., Lewerentz M., Quack, M. 1984. Chem. Phys. Lett. 109: 563; Bagott,J.E., Chuang, M. C., Zare, R. N., Dübal, H. R., Quack, M. 1985. J. Chem. Phys. 82: 1186; Segall, J., Zare, R. N., Dübal, H. R., Lewerentz, M., Quack, M. 1985. J. Chem. Phys. 82: 1186

[16] Chang, J., Wyatt, R. E. 1986. J. Chem. Phys. 85: 1826

[17] Lederman, S. M., Marcus, R. A. 1988. J. Chem. Phys. 88: 6312 [18] Tietz, J. V., Chu, S. I. 1983. Chem. Phys. Lett. 101: 446

[19] Nauts, A, Wyatt, R. E. 1983. Phys. Rev. Lett. 51: 2238; 1984. Phys. Rev. A 30: 872

[20] Jolicard, G., Périé, J. 1993. Advances in multi-photon processes and spectroscopy Vol.8, p. 1-126 Singapore: World Scientific

[21] Peskin, U., Moiseyev, N. 1993 J. Chem. Phys. 99: 4590 [22] Peskin, N., Alon, O. E., Moiseyev, N. 1994 J. Chem. Phys. 100: 7310 [23] Moiseyev, N., Hirschefelder, J. O. 1988 J. Chem. Phys. 88: 1063 [24] Baslev, E., Combes, J. M. 1971. Commun. Math. Phys. 22: 280;

Aguilar, J., Combes, J. M. 1971. Commun. Math. Phys. 22: 265; Simon, B. 1972. Commun. Math. Phys. 27: 1

[25] Simon, B. 1979. Phys. Lett. A 71: 211

[26] Rom, N., Engdahl, E., Moiseyev, N. 1990. J. Chem. Phys. 93: 3413 [27] Jolicard, G., Austin, E. J. 1985. Chem. Phys. Letters. 21: 106; Joli-

card, G., Austin, E. J. 1986. Chem. Phys. 103: 295

[28] Jolicard, G. 1989. J. Chem. Phys. 90: 2320

[29] Jolicard, G., Killingbeck, J. P. 1994 InDomain-based prallelism and problem solving decomposition methods in computational science and engineering S.I.A.M. Book (unpublished)

[30] Malrieu, J. P., Durand, Ph., Daudey, J. P. 1985. J. Phys. A Math. Gen. 18: 809

[31] Killingbeck, J. P. 1991 Microcomputer algorithms-Action From Algebra. Bristol: Adam Hilger

[32] Killingbeck, J. P., Jolicard, G. 1992. J. Phys. A Math. Gen. 25: 6455 [33] Jolicard, G. 1987. Chem. Phys. 115: 57 [34] Jolicard, G., Moiseyev, N. J. Chem. Phys. (unpublished)

[35] Howland, J. S. 1974. Math. Ann. 207: 315

[36] Stine, J. R., Noid, D. W. 1979. Opt. Comm. 31: 161

[37] Beswick, J. A., Jortner, J. 1981 Adv. Chem. Phys. 47: 363

[38] Gray, S. K. J. Chem. Phys. 1992 96: 6543

[39] Moiseyev, N., Peskin, U. 1990. Phys. Rev. A42: 255

[40] Jolicard, G., Billing, G. D. 1992. J. Chem. Phys. 97: 997

[41] Balint-Kurti, G. G., Dixon, R. N., Marston, C. C., Mulholland, A. J. 1991. Comp. Phys. Comm. 63: 126

[42] Atabek, O., Jolicard, G. 1994. Phys. Rev. A 49: 1186

[43] Wyatt, R. E., Iung, C., Leforestier, C. 1992. J. Chem. Phys. 97: 3458; ibid 97: 3477