10

Φ8cmcsc10

CHAPTER 1

WAVE OPERATORS AND ACTIVE SUBSPACES; TOOLS FOR THE

SIMPLIFIED DYNAMICAL DESCRIPTION OF QUANTUM PROCESSES INVOLVING MANY-DIMENSIONAL STATE SPACES

GEORGES JOLICARD^{∗} AND JOHN P. KILLINGBECK^{†}

Abstract. Time-dependent and stationary wave operators are presented as tools to define active spaces and simplified dynamics for the integration of the time-dependent Schrödinger equation in large quantum spaces. Within this framework a new light is thrown on the duality between time-dependent and time-independent approaches and a generalized version of the adiabatic theorem is given. For the Floquet treatment of photodissociation processes, the choice of the relevant subspaces and the construction of the effective Hamiltonians are carried out using the Bloch wave operator techniques. Iterative solutions of the basic equations associated with these wave operators are given, based on Jacobi, Gauss-Seidel and variational schemes.

1. Molecular collisions: their computational representation and analysis..

1.1. Statement of the problem. A basic numerical problem for chemistry and molecular physics is the study of collisions. A general collision process leads from some entrance (initial) continuum to some exit (final) continuum. In the case of a full collision, energy exchanges appear between molecular bound states and molecular continua, leading to products characterized by new energy distributions, or by new chemical structures in the case of reactions and dissociations. Half-collision processes are those in which one continuum (entrance or exit) is a photon continuum and the other one is an ionization or dissociation continuum. Typical half-collision processes are photoionization and photodissociation and their reverse processes, radiative recombination and radiative association.

In most cases these half-collision processes can be described by a time-dependent Schrödinger equation [1].

ih̄\frac{∂}{∂t}Ψ(q, r, t) = H(q, r, t)Ψ(q, r, t) (1)

where r represents one or many dissociative coordinates, and q a group of bound coordinates. The Hamiltonian H is in part a differential operator with respect to q and r. It depends explicitly on the time when some of the molecular coordinates or the electromagnetic field are considered as classical variables [2, 3].

Many approaches involving the theory of wave packet dynamics have been proposed during the last ten years to treat equation (1) [4, 5, 6]. These treatments discretize the continua by working with a bounded range of radial coordinates in conjunction with finite basis sets of square integrable functions. They are well suited for numerical integration methods; the construction of both the partial propagator U(t + dt, t) and the global propagator U(t, t_{0}) involves the repeated formation of the matrix-vector product HΨ as the principal operation in the description of the time evolution process.

Ψ(q, r, t) = U(t, t_{0})Ψ(q, r, t_{0}) (2)

Traditional methods of calculating the partial propagator include the second order differencing scheme, the split operator method and the short iterative Lanczos method [7, 8]; the global propagator can be calculated using a scheme based on Chebyshev polynomial expansions [9]. These methods encounter two typical difficulties. The first arises from the discretization of the continua which effectively implies the addition of absorbing boundary conditions. The Hamiltonian operator then becomes nonhermitian, which often consequently reduces the efficiency of the propagators [10].

The second difficulty is the large CPU time which characterizes wave packet propagation computations. This often becomes prohibitively long even when Fast Fourier Transforms (FFT) are used to transform Discrete Variable Representation (DVR) into a Finite Basis Representation (FBR) . The difficulty can be partially overcome by deriving a time-independent treatment from the time-dependent one. For timedependent Hamiltonian this is unfortunately a difficult task, since the analytical form of U(t, t_{0}) is unknown. A solution has been proposed by Shirley within the framework of the Floquet theory for time periodic Hamiltonians [11]. Recently a promising approach, the (t, t′) method, has been proposed by Peskin and Moiseyev [12] for time non-periodic Hamiltonians.

∗ Laboratoire de Physique Moléculaire (U.R.A. 772 du C.N.R.S.). Faculté des Sciences et des Tecnniques - La Bouloie 25030 Besancon Cedex, France. New permanent address: Observatoire de Besançon. Université de Franche Comté, 41, Av. de l’Observatoire, BP 1615, 25010 Besançon Cedex, France

† Department of Applied Mathematics, The University, Cottingham Road HULL, HU6 7RX, United Kingdom 1

The study of collisions poses a second challenge to data storage and CPU time limitations with the calculation of resonance states (either purely molecular resonance states or resonances induced by the field-matter interaction) [13]. These states are eigenvectors of the stationary Hamiltonian having particular asymptotic features: they are purely outgoing and exponentially diverging functions and correspond to poles of the resolvent operator, lying in the complex lower half-plane [14]. Thus they can be characterized by the equations

H(q, r)Φ(q, r) = EΦ(q, r)
(3)  with lim_{r→∞} Φ(q, r) = ∑_{v} ϕ_{v}(q)A_{v} exp(ik_{v}r)

where the asymptotic fragment wavefunctions Φ_{v}(q) form a complete set. E and k_{v} are complex numbers with negative imaginary parts. A_{v} is a complex amplitude.

This concept is of great importance because the resonance states not only strongly influence full collision processes but also constitute an essential element for the description of half collisions (e.g. photodissociation and photoionisation, unimolecular reactions of van der Waals molecules ...) [15]. Unfortunately the large size of the basis set involved does not usually allow the use of standard matrix diagonalization programs so that this stationary eigenproblem, like the preceeding time-dependent problem, introduces severe difficulties in terms of central memory and CPU time requirements. The aim of our paper is to analyze these difficulties using the concept of stationary and time-dependent wave operators. It will be shown that these operators are efficient tools to define active spaces of reduced size, encapsulating the essential features of the dynamics [16]. It will also be shown that this approach casts a new light on the duality between time-dependent and time-independent approaches and that it offers a direct way to transform the time-dependent treatments into timeindependent ones by generalizing the Floquet operator concept. Particular attention will be paid to recursive methods for the solution of the basic nonlinear differential equations associated with the various wave operators which arise in the theory.

1.2. The boundary conditions. The full study of molecular collisions in a finite basis representation is usually carried out in three steps:

• The first step defines the finite L^{2} basis set representation and asymptotic absorbing boundaries. In this same volume Bacic introduces the DVR representation and presents associated strategies for reducing the dimensionality of the Hamiltonian matrix.

• The second step integrates the matrix differential equations obtained by projection of the relevant equations of motion onto the selected basis set.

• The third step makes an asymptotic analysis of the fragments resulting from

the interaction. These three steps are in fact closely correlated, and the merits or defects of each one have direct consequences for the others. The boundary conditions are, for example, of primary importance in the resonance state eigenproblem. This is because these states are nonstationary states that transfer quantum flux outside the region of interaction. The adjoint eigenvectors ϕ^{+}, constituting with the ϕ a biorthogonal basis set, are the complex conjugates ϕ^{∗} of the resonance eigenvectors, so that the incoming flux associated with ϕ^{+} is equal in amplitude to the outgoing flux associated with the corresponding ϕ states, and H is thus non-Hermitian with respect to these vectors [13],

ϕ^{†} = ϕ^{∗} ≠ ϕ (4)

This property is inconsistent with the use of a finite L^{2} basis set which produces an hermitian H matrix. One can overcome this difficulty by using a mathematical artifact, namely by introducing near the upper boundary an additional absorbing imaginary potential term (the optical potential) which includes a set of adjustable parameters λ [17],

H(q, r) → H(q, r) − iV_{opt}(λ, r) (5)

Many studies have shown that this model is appropriate for the resonance eigenproblem and that precise resonance state eigenvalues can be obtained through a stabilization procedure with respect to the parameters λ [18, 19]. A similar choice of boundary conditions can be made in the case of wave packet propagation but the shape and the position of the optical potential are correlated with the asymptotic analysis of the fragments. The most efficient methods of analysis give energy-resolved results over the full spectrum of the wave packet by analysing, at a fixed asymptotic position, the outgoing wave function or the quantum flux [20, 21]. However the asymptotic absorption of the wave packet has non-local effects so that the precision of the analysis essentially depends on the radial extent of both the optical potential and the outgoing wave packet [21]. This analysis cannot be considered for resonance states having an infinite extent but fortunately these states are characterized by an identical exponential time dependence of the flow in each open channel. Thus the ratios of the integrated flows are equal to the ratios of the instantaneous flows at a given arbitrary time. Further the branching ratios are simply expressed in terms of the asymptotic amplitudes of the resonances states in the various open channels [22, 23].

The preceding comments reveal that any new theoretical treatment such as the wave operator approach must possess some specific features in order to be applied successfully in collision theory. For example, it should be able to work with large non hermitian matrices and also be adaptable to permit variational stabilisation procedures, without repeating the full calculational procedure. In solving the eigenproblem it should supply the eigenvectors, which directly participate in the quantum flow calculations.

2. The wave operator theory.

2.1. The time-dependent wave operator concept. Let us consider again the time-dependent evolution equation (1) and let us assume that there exists some active subspace S_{0} with a dimension small compared to that of its orthogonal complement space S_{0}^{+} usually called the complementary space. The projector of the active space is then denoted by P_{0} and that of the complementary space by Q_{0}. By assuming that an orthogonal basis |i > ⊕ |α > spans the whole of the vector space required to describe the system we can set

P_{0} = ∑_{i} |i >< i| , Q_{0} = ∑_{α} |α >< α| (6)

The solution of the equation of motion (1) is given formally in equation (2) by introducing the propagator U(t, t_{0}), sometimes written more explicityl as U(t, t_{0}; H), since various different effective Hamiltonians are introduced in approximate treatments. The active subspace component of U(t, t_{0}; H) is P_{0}U(t, t_{0}; H)P_{0}. In general, using an obvious abbreviated notation, U will be the sum of four components

U = P_{0}UP_{0} + P_{0}UQ_{0} + Q_{0}UP_{0} + Q_{0}UQ_{0} (7)

where we can regared each operator as being represented in the full space by a partitioned square matrix with three zero blocks. If, for example, the full space is of the dimension N and the subspace S_{0} is of dimension M, then P_{0}UP_{0} will be represented by a matrix with a non-zero block of M ×M type and three zero blocks which complete the full N × N matrix. The M × M non-zero part of P_{0}UP_{0} will have an M × M inverse, which, when augmented by three zero blocks to form an N × N matrix, is denoted (in a slight abuse of the usual conventions) by (P_{0}UP_{0})^{−1}. From the matrix multiplication properties of partitioned matrices, represented by the projection operator properties P_{0}P_{0} = P_{0}, P_{0}Q_{0} = 0, etc, it follows that we then obtain the result

U(P_{0}UP_{0})^{−1} = P_{0} + Q_{0}UP_{0}(P_{0}UP_{0})^{−1} (8)

on remembering that in the full space S_{0} is represented by a unit M × M matrix plus three appropriate zero blocks. Reverting to the full notation, we denote the operator on the left of (8) by Ω(t, t_{0}) and call it the time-dependent wave operator [27], while the second operator on the right is called the reduced wave operator X(t, t_{0}) ; equation (8) then becomes

Ω(t, t_{0}) = P_{0} + X(t, t_{0}) (9)

with

Ω(t, t_{0}) = U(t, t_{0}; H) × (P_{0}U(t, t_{0}; H)P_{0})^{−1} (10)

X(t, t_{0}) is an off-diagonal part of Ω(t, t_{0}) ; in the N × N matrix representation it has non-zero elements only in the lower left (N − M) × M block. In numerical computation, of course, it is only the relatively small number of non-zero elements of X(t, t_{0}) and Ω(t, t_{0}) which need to be calculated. In the time-independent case, the formalism above shows that, if we can find a transformation of the basis vectors in the N dimensional space which will render X(ie Q_{0}UP_{0}) zero, then the operator U will not couple the space S_{0} and the complementary space, and we shall have produced an isolated M × M matrix block P_{0}UP_{0} which gives a more tractable problem of reduced dimension. A matrix form of this approach has been applied directly to matrix eigenvalue problems by some workers [24, 25], although the projection operator formalism is that predominantly used by workers in wave operator theory. In the time-dependent case, of course, the quantities (including the space S_{0}) are changing with time and a single fixed transformation will not serve to decouple the two spaces permanently by rendering X = Q_{0}U(P_{0}UP_{0})^{−1} zero at all times. It turns out that in principle one can construct a sequence of transformations, involving an effective Hamiltonian H_{eff}(t), which will keep on maintaining this decoupling. As one would expect, what is obtained is actually an equation of motion which must be satisfied to maintain the decoupling, since the formalism is essentially a recasting of the basic time-dependent Schrödinger equation. The starting point of the method is to combine equations (9) and (10) to give

U(t, t_{0}; H)P_{0} = (P_{0} + X(t, t_{0}))(P_{0}U(t, t_{0}; H)P_{0}) (11) which may be seen to be valid by recalling the partitioned matrix representation described above. Interpreting that representation in dynamical terms, P_{0}UP_{0} describes a first evolution step within the active space S_{0}, while the X term (involving Q_{0}UP_{0}) describes a transition from S_{0} to the complementary space without any reverse path (a reverse path would involve P_{0}UQ_{0}). The equations presented so far are purely formal insofar as the operators Ω and X are expressed via the evolution operator U. The term P_{0}UP_{0} in (10) is the S_{0} space projection of U, and the aim of the theory developed here is to regard it as the exact time development operator associated with an effective Hamiltonian H_{eff} via the standard Schrödinger equation of motion. Accordingly, we study a small time interval t_{n} → t_{n+1} and set

P_{0}U(t_{n+1}, t_{n}; H)P_{0} = U(t_{n+1}, t_{n}; H_{eff}) (12)

with

ih̄\frac{∂}{∂t}U(t, t_{n}; H_{eff}) = H_{eff}U(t, t_{n}; H_{eff}) (13)

The guiding idea behind the formalism is that at t = t_{n} we start with X(t) = 0, ie no coupling between S_{0} and its complementary space. As the time-development begins to produce this coupling, we arrange the transformations of the basis vectors to cancel this coupling. The formal mathematics behind the calculation is based on the identity

H = (H − ih̄\frac{∂X}{∂t}) + ih̄\frac{∂X}{∂t} (14)

Because of the particular block form of X, described earlier, it can be shown that the operator ih̄\frac{∂X}{∂t}, regarded as a kind of fictitious Hamiltonian, has a very simple development operator,

U(t, t_{0}; ih̄∂X/∂t) = 1 + X(t, t_{0}) (15)

The block form of the matrices also shows that X^{2} = 0, so that (1 − X)(1 + X) = 1 and 1 − X is the inverse of 1 + X. At this stage we make appeal to the theory of itnermediate representations [26], according to which the development operator U_{A+B} for an operator sum A + B can be written as

U_{A+B} = U_{B}U(U_{B}^{−1}AU_{B}) (16)

For the sum in equation (14) this gives the result

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

with

H_{eff}(t) = (1 − X(t, t_{0}))(H − ih̄\frac{∂X}{∂t})(1 + X(t, t_{0})) (18)

If we now specify X more precisely by demanding that it should be chosen to make Q_{0}H_{eff}P_{0} zero, we obtain an equation of motion for X ;

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

The quantity of principal interest now is the projection of H_{eff} within the active space S_{0}. The P_{0} projection operator on the left of equation (12 ) separates out the two left hand blocks. We do not need to specify more closely the other two parts of H_{eff} in the block matrix structure. Remembering the location of the various blocks in the matrix structure, we can write the active space effective Hamiltonian as

H_{eff}(t) = P_{0}H(t)(P_{0} + X(t, t_{0})) (20)

The three equations (17), (19) and (20) are the three basic equations of the timedependent wave operator theory, and are handled numerically by using small discrete time steps t_{n} → t_{n+1}, as noted in connection with equation (12). Equation (19) which defines the reduced wave operator can be written differently after some algebraic manipulations, [28], in particular:

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

Alternatively, by combining eqns. (18) and (21) we obtain

(H − ih̄\frac{∂}{∂t})Ω(t, t_{0}) = Ω(t, t_{0})H_{eff}(t)
(22)  = Ω(t, t_{0})H(t)Ω(t, t_{0})

The three different forms (19), (21) and (22) will be used in the next chapters and related to other theoretical treatments and concepts.

2.2. The stationary wave operator. Wave operator theory was largely developed within quantum chemistry to obtain new effective Hamiltonians which can provide accurate molecular energies and easily interpretable wavefunctions [29]. The basic idea is to define an effective Hamiltonian H_{eff} acting in a reduced model space S_{0} such that a subset of the exact energies of the Hamiltonian H coincide with those of H_{eff}. The choice made by Bloch [30, 31] corresponds to setting:

Ω = P (P_{0}P P_{0})^{−1} (23)

where P_{0} is, as previously, the projection operator of the model space S_{0} and P the projection operator of the corresponding target space. The eigenvectors |λ_{i}^{0}⟩ of H_{eff} are in this case the projections into the model space of the corresponding eigenvectors |λ_{i}⟩ of H and Ω appears in this context as the operator which transforms these eigenvectors projected into S_{0} back to the exact eigenvectors :

|λ_{i}^{0}⟩ = P_{0}|λ_{i}⟩
(24)  |λ_{i}⟩ = Ω|λ_{i}^{0}⟩

The theory of Bloch and its extensions [32], as developed in numerous quantum chemical studies, are purely stationary state methods; however, they can be derived from the general time-dependent wave operator theory presented in Section 2.1. Let us assume that the operator H is generated from the unperturbed operator H^{0} by adiabatically introducing the perturbation V (≡ H − H^{0}):

H(t) = lim_{γ→0^{+}}(H^{0} + V exp(γt))
H^{0} = H(t = −∞)
(25)  H = H(t = 0);

Fig. 1. Schematic representation of the effective Hamiltonian calculation by the wave operator transformation

then the derivative ∂X/∂t tends to zero at any time, and eq. (19) at t = 0 leads to :

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

This equation which is illustrated in Fig.1 is the basic equation of Bloch’s wave operator theory. Most of the existing literature presents this basic equation in the alternative form [29]:

HΩ = ΩHΩ (27)

However, (26) has a pedagogical advantage, since it shows that the reduced wave operator X is simply the non-diagonal part of a non-unitary and non-singular transformation (T = 1 + X) which possesses a trivial inverse (T^{−1} = 1 − X) [33]. This transformation isolates the S_{0} space by cancelling out the couplings which connect S_{0} to S_{0}^{+}, so that the eigenvalues originating within the space S_{0} are simply obtained by diagonalising H_{eff} = P_{0}HΩ.

3. Solution of the wave operator equations.

3.1. Definition of the active space. Equation (17) provides an exact description of the quantum evolution issuing from the space S_{0}. This equation is complicated when compared with the standard time-dependent Schrödinger equation, since Ω (orX) are the solutions of the non-linear equation (19). Thus its application must be restricted to processes which are confined within target spaces induced from relatively small active spaces. Fortunately, wave operator theory is an efficient tool for defining and constructing suitable active spaces, while permitting simultaneous generalization of the adiabatic theorem. A process will be said to be adiabatically evolving with respect to the active space S_{0} if the corresponding reduced wave operator at any time obeys the restriction that it has a suitably small time derivative,

∂X(t, t_{0})/∂t ∼ 0, (28)

i.e., if the time-dependent wave operator essentially converges to the instantaneous stationary one. When S_{0} at t_{0} consists of a single basis vector |λ^{0}⟩ which is itself one nondegenerate eigenvector of H(t_{0}), the adiabatic evolution process leads to :

Ψ(t) ∼ Ω exp(\frac{1}{ih̄} ∫_{t_{0}}^{t} E(t′)dt′) (29)

Fig. 2. Schematic representation in the dressed picture of the potential surfaces involved in the dissociation process of the HF molecule. The field intensity is I = 1T W/cm^{2} and the frequency (ω = 0.02914 a.u.) is tuned to the frequency of the transition (v = 0 → v = 2). The surfaces n = −1 and n = 1 are the dressed surfaces associated with the absorption and emission of one photon. The figure shows the bound and the continuum states retained by the wave operator selection algorithm and constituting the model space

The wave operator Ω is here the instantaneous eigenvector of H(t) at time t and E(t) is the corresponding eigenvalue. As Ω(t) satisfies an intermediate normalisation condition (⟨λ^{0}|Ω(t)⟩ = 1), E(t) posseses an imaginary component which ensures norm conservation when H(t) is hermitian. The wave operator formulation can also resolve more intricate situations in which the evolution, while appearing adiabatic with respect to the transitions between S_{0} and S_{0}^{+}, exhibit strong nonadiabatic or even sudden behaviour inside the active space. A direct definition of this active space can then be obtained by first defining the initial space spanned by the initial vector |λ^{0}⟩ (from a practical point of view the basis used should be chosen so that |λ^{0}⟩ has projections on only one or a few basis vectors). Then the wave operator associated with this initial space is derived in an iterative treatment (see section 3.2) and the vectors of the whole space are sorted in the order of decreasing absolute values of the components of this wave operator. Finally the active space, the dimension N_{0} of which is imposed, is built by taking the N_{0} first sorted vectors [34]. Accordingly the very concept of a reduced active space disappears when the number of the nonnegligible components becomes too large. This procedure of selection of an active space has been successfully applied by Wyatt et al. [35] to calculations of the overtone relaxation of benzene, where it appears simpler than the well known artificial intelligence methods [36], since search algorithms and evaluation functions are not required. The same procedure has also been used in the theory of photorelaxation and photodissociation processes [34, 37] . Figure 2 illustrates this for the dissociation of a Morse oscillator produced by a monochromatic field.

3.2. Integration of the stationary Bloch equation. The stationary wave operator Ω (or X ) plays a central role. It is used to define the active space concept (section 3.1). It directly participates in the construction of H_{eff} and in the relationships between the eigenvectors of H_{eff} and the exact eigenvectors of (24). It also participates in the reduced dynamical equations inside the target space (see section 4.1 ). Thus several linear and quadratic iterative methods have been proposed to solve the non-linear equation (26) [38]. (Some of them are intractable or approximate, e.g. they assume a diagonal structure of H in the complementary space S_{0}^{+}). The recursive distorted wave approximation (RDWA) [39] leads to an iterative process described by the equations

X = lim_{N→∞} X^{N}
(30)  X^{N} = X^{N−1} + ∑_{β∈S_{0}^{+}} ∑_{i∈S_{0}} |β >< i|(\frac{< β|H^{N}|i >}{< i|H^{N}|i > − < β|H^{N}|β >})

with

H^{N=0} = H
(31)  H^{N>0} = (1 − X^{(N−1)})H(1 + X^{(N−1)})

At each step, most of the required computational CPU time is devoted to forming the matrix-vector product HX. The trial value X^{N=0} = 0 is often used but in fact any arbitrary value can be selected; it is obvious that the choice of a trial operator X^{N=0} near to the converged solution increases the speed of convergence.

Divergence of the series (30) appears in the strong coupling regime, particularly when intruder states associated with avoided crossing are involved. To partially overcome these difficulties a treatment called the single cycle method (SCM) has recently been proposed in conjunction with the Gauss-Seidel method.[25]. The SCM method modifies the RDWA algorithm by using the fact that the family of special transformations 1 + X used is such that the simple product rule (1 + X)(1 + Y) = 1 + X + Y holds ; the result (1 − X) = (1 + X)^{−1} is a special case of this property. The SCM uses the decomposition,

1 + X^{N} = (1 + ΔX^{N})(1 + X^{N−1})
(32)  with ΔX^{N} = X^{N} − X^{N−1}

and then expands each transformation (1 + ΔX^{N}) as a series of elementary transformations:

1 + ΔX^{N} = ∏_{i∈S_{0}} ∏_{α∈S_{0}^{+}}[1 + |α >< i|ΔX_{α,i}^{N} ] (33)

where i is the column index of the wave operator in the model space and α the row index in the outer space S_{0}^{+}. The improvement comes from the fact that after each partial transformation, the H elements are modified before being used for the next transformation. This produces a larger radius of convergence but at the cost of a more complicated structure of the iteration. If a generalized iteration index (L ≡ (N, α, i)), incremented after each partial transformation, is introduced instead of N, the basic equation (30) remains identical but (31) is replaced by

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

It thus involves a step by step relationship (H^{L} = f(H^{L−1})) instead of a direct relationship H^{L} = f(H)). Fortunately the iteration remains tractable, since only a small part of the H matrix, namely the diagonal part and the band Q_{0}H^{L}P_{0}, is modified at each step and needs to be stored. Details of this iteration are given elsewhere [37].

3.3. Integration of the time dependent wave operator equation. The time-dependent wave operator is the main concept involved in our study insofar as it gives rise to the stationary wave operator. It appears in equation (19) as the solution of a non-linear differential equation. An iterative procedure, closely related to those used for the stationary wave operator, was proposed [40] to solve this equation. Unfortunately the slow convergence of this solution makes it not well adapted to the study of evolution over long time periods. A better approach is obtained by using second order differencing to approximate the time derivative in (19) and to preserve time reversal symmetry. With a variable time step, one gets :

X(t + ϵdt) = ϵ^{2}X(t − dt) + (1 − ϵ^{2})X(t)
(35)  + (ϵ + ϵ^{2})dt\frac{∂X(t)}{∂t} + 0(dt^{3})

with

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

The parameter ϵ connected with the time step variation is confined in a small interval centered on ϵ = 1. Its value is defined at each step with regard to the amplitude derivative ∂X/∂t. From a practical point of view it must be recalled that the major part of the CPU time used to integrate (35) is devoted to forming the matrix-vector product H × X appearing in the expression for ∂X/∂t.

3.4. From the stationary to the time-dependent problem. A basic feature of the wave operator theory is its capacity to establish relationships between the time-dependent and the time-independent formalisms. The fact that the Bloch wave operator has been derived from the time-dependent one, in the adiabatic limit, can for instance be used in stationary problems to overcome convergence difficulties of the iterative scheme in the strong coupling regime. As noted previously, the initial trial operator X^{N=0} introduced in these iterations is arbitrary; that derived from a time-dependent procedure by establishing the couplings nearly adiabatically can often be expected to be very close to the desired exact one.

The equation (35) is well suited for the time integration of a slowly varying evolution process because the features leading to large time derivatives at the adiabatic limit have been suppressed. Equation (29) shows for example that the quickly varying phase terms have been extracted and separated from the wave operator. Moreover, the large couplings which could appear near the avoided crossings are suppressed by this formulation if the states participing in these avoided crossings are included in the

Fig. 3. Lower part of the spectrum of two harmonic oscillators with linear and quadratic perturbations versus the perturbation parameter λ. The case studied by the wave operator formulation and presented in fig.4 corresponds to the value λ = 0.4 and to the 4^{th} group of quasidegenerate eigenvectors

active space S_{0}. A simple illustration is presented in Figs. (3) and (4). This corresponds to two harmonic oscillators in near resonances which are linearly perturbed and quadratically coupled. Fig. 3 presents the lower part of the energy spectrum of this system. The perturbation amplitude is established as a linearly increasing function of time on the interval [t = 0, t = N/ν] where ν is the smaller of the two vibrational frequencies. In the present case this amplitude varies from 0 up to 0.40. The wave operator corresponding to the 4^{th} group of quasidegenerate eigenvectors in the energy scale (see Fig. 3) has been calculated and the four corresponding approximate eigenvectors derived from eq.(24). Figure 4 presents the factors

log_{10} | < ϕ_{j}^{exact}|ϕ_{i}^{appr} > −δ_{ij}| (37)

for i = 1 to 4 versus the index j. These factors, which estimate the deviation from a purely stationary case, are unexpectedly small for N = 3 (a time interval with only three vibrational oscillations).

3.5. From the time-dependent problem to the stationary one. Equations 21 and 22 provide a means to transform the time-dependent dynamical problem into a stationary problem [41]. Equation (21) is the equivalent of the stationary equation (26) in which the Floquet-type operator H − ih̄∂/∂t is substituted for the molecular Hamiltonian H. In the case of Hamiltonians which are periodic in time (H(t) = H(t+ τ)), the Floquet theorem establishes the existence of stationary solutions, themselves

Fig. 4. The non-adiabaticity factors: log_{10} | < ϕ_{j}^{exact}|ϕ_{i}^{appr} > −δ_{ij}| for the four states (i = 1 to 4)(∘, +, ×, ∗ ) of the 4^{th} group of quasidegenerate states. The factors are shown versus the index j which denotes the harmonic basis vectors sorted in the order of increasing energy values.

periodic in time, i.e. [42] :

(H − ih̄\frac{∂}{∂t})Φ_{k}(t) = ϵ_{k}Φ_{k}(t)
(38)  with Φ_{k}(t + τ) = Φ_{k}(t)

The quasienergy ϵ_{k} is unique up to a multiple of 2πn/τ, and in accord with the Floquet theorem the quasienergy state can be written in the form

Ψ(t) = exp(−iϵ_{k}t/h̄)Φ_{k}(t) (39)

The two equations (38) and (39) can be identified with eqs. (22) and (17), respectively, for the case of a one-dimensional active space. In this case the effective Hamiltonian is the time independent scalar ϵ_{k} and the reduced wave operator itself becomes a quasivector obeying the periodic conditions :

X(t_{0} + nτ, t_{0}) = 0 (40)

However, the time-dependent wave operator theory adopts a more general point of view by generalizing the quasivector concept to degenerate active space situations. This generalization, which transforms the quasi-vector into a wave operator and the quasi-energy into a quasi-effective Hamiltonian operator improves the calculational scheme for the quasi-vectors when many quasi-vectors in near resonance are collected into the space S_{0}.

3.6. The concept of an evolutive nonorthogonal basis set. The wave operator formalism is confronted by difficulties when the initial state projects onto a large initial space. This is the case with wavepacket propagation calculations when DVR basis sets are used. The wave packet dilates during the propagation; it quickly spreads out in space and takes non-zero values on the whole grid. As soon as the dimension of the subspace S_{0} becomes non-negligible with respect to the dimension of S_{0}^{+}, the treatment presented in the subsection 3.3 is no longer useful. Nevertheless this dilemma can be overcome by defining, at the initial time of each propagation step, a new nonorthogonal basis set. If the primary orthonormal basis set is denoted by ϕ_{j}, j = 1 to N, then the new basis set (ϕ′) can be built up as follows:

ϕj^{′} = ϕ_{j}forj ≠ p
(41)  ϕp^{′} = Ψ(t_{0})

where Ψ(t_{0}) is the exact wave function at the new initial instant t_{0} (cf. eq. 1) and p is the index associated with the largest component < ϕ_{j}|Ψ(t_{0}) >. When the wave function is expanded in the new basis set:

Ψ(t) = ϕ′A′(t) for t ≥ t_{0} (42)

the resulting matrix equation is :

ih̄\frac{∂}{∂t}A′(t) = (B^{−1}HB)A′(t)
(43)  with B = ϕ^{†} × ϕ′

The difficulty generated by the large initial space has disappeared in this new scheme since A′(t_{0}) has only one nonzero component :

Aj^{′}(t_{0}) = δ_{jp} (44)

Moreover the matrix B which define the new Hamiltonian matrix has a simple form. It equals the identity matrix except for the p^{th} column, for which B_{jp} =< ϕ_{j}|Ψ(t) >. The matrix B^{−1} has the same structure but with:

(B^{−1})_{pp} = 1/B_{pp}
(45)  (B^{−1})_{jp} = −B_{jp}/B_{pp} for j ≠ p

Thus going from H to B^{−1}HB only requires a matrix-vector product.

4. Illustrative numerical examples : The photodissociation of molecules by an intense laser field. Molecular photodissociation is a process which, while apparently complicated, can often be described theoretically using a small target space formed by only a few quasi-vectors. The wave operator formalism, adapted in accord with Floquet theory, is then well suited for its description.

4.1. Projected dynamical equation in the quasienergy formulation. Consider a molecular system driven by a periodic external field (i.e. a continuous field) with frequency ω and period τ(= 2π/ω) . Let us introduce the Floquet state nomenclature |α, n > where |α > designates the elements of a complete orthonormal set of eigenfunctions of the unperturbed molecular Hamiltonian and |n > are the Fourier

Fig. 5. Schematic representation of the time relaxation before and after the wave operator transformation

vectors such that < t|n >= exp(inωt). One can then express the transition probability (α → β) averaged over the initial phases of the field seen by the molecular system, as

P_{α→β} = ∑_{n} | < n, β| exp(−iHΔt/h̄)|0, α > |^{2} (46)

where H(≡ H −ih̄∂/∂t) is the Floquet matrix. This exact expression for the evolution can be simplified by first defining a reduced active space S_{0}, as explained in section (3.1); Figure 2 illustrates this selection of S_{0} for a model system. Then the wave operators Ω_{P} and Ω_{Q} associated with the two complementary spaces S_{0} and S_{0}^{+} are introduced and the non-unitary transformation

T = Ω_{P} + Ω_{Q} (47)

is applied to the evolution operator (exp(−iHΔt/h̄)P_{0}), as shown in Fig. 5. The simplified relaxation scheme is obtained by neglecting the relaxation path through the space S_{0}^{+}. A detailed demonstration in reference [34] leads to the final equation :

P_{α→β}(Δt) = ∑_{n} | ∑_{η∈S_{0}} < n, β|Φ_{η} >< Φ_{η}|0, α > exp(−iϵ_{η}Δt/h̄)|^{2} (48)

in which the overlap matrices and the quasi eigenvalues have simple expressions in terms of the wave operator Ω_{P}. Eq. (48) reveals that the wave operator formulation simplifies the dynamical description by confining it within the space spanned by the quasi-vectors derived from space S_{0}. An interesting situation, corresponding to the model presented on Fig. 2, is shown on Fig. 6. The conditions leading to appreciable dissociation probabilities from the ground state (v = 0) correspond to the use of doubly degenerate target spaces in which the quasi-vector state Φ(n = 0, v = 0) mixes with a second quasi-vector Φ(n = −1, v) in near resonance [37].

Fig. 6. Photodissociation of the HF molecule in the frequency range [0.015a.u., 0.048a.u.] calculated at t = 10^{9}a.u. for a continuous wave laser of intensity I = 1T W/cm^{2}. Upper panel: The non negligible overlaps between the initial state |n = 0, v = 0 > and eigenvectors |Φ_{n,v} > of HF. Middle panel: Imaginary part of the eigenvalues of the states |Φ_{n,v} >. Lower panel: The dissociation probabilities.

4.2. Extension to photodissociation processes using short laser pulses. The wave operator-Floquet formalism can be generalized to non-periodic Hamiltonians to describe dissociations induced by short laser pulses of form

E(t) = F(t) cos(ωt) (49)

In a discretized model the range of variation [0, F^{max}] is divided into N elementary intervals of the same width ΔF and a corresponding partitioning of the time interval is produced (t_{i}) with F(t_{i}) = iΔF, i ≤ N and F(t_{i}) = (2N − i)ΔF, N < i ≤ 2N. The continuous function is then approximated by the discontinuous function which takes

Fig. 7. Photodissociation of the H_{2}^{+} molecule: Complex energies of the central Floquet block (n = 0) without field-matter interaction. A group of (2 × 150) functions has been used to span the two electronic states (1sσ_{g})^{1}Σ_{g}^{+} and (2pσ_{u})^{1}Σ_{u}^{+} which participate. (+ denotes the bound states of the g potential surface. ∘ and ∗ denote discretized continuum states of the potential surfaces g and u respectively).

the constant value F(t_{i}) when t is between the two adjacent points t_{i} and t_{i+1}. An active space S_{0} is defined as for the case of a continuous laser and for each discretized value of the field magnitude a target space is associated with the active space :

S_{0} → S^{(j)} (50)

By assuming that at the discontinuity between two adjacent values of the electric field, namely j and j + 1, the wavefunction located in the subspace S^{(j)} is completely projected into the subspace S^{(j+1)}, the wave function obtained at any time t is [43]:

Ψ(t) = ∑_{γ} |Φ_{γ_{n}} > exp(−iE_{γ}^{n}Δt_{n}/h̄) ∑_{β} < Φ_{γ_{n}}|Φ_{β_{n−1}} > exp(−iE_{β}^{n−1}Δt_{n−1}/h̄)
(51)  … exp(−iE_{μ}^{1}Δt_{1}/h̄) < Φ_{μ_{1}}|0, α >

where |0, α > denotes the initial molecular state, and the ∑ summations are limited to the target space S^{(j)} previously constructed. A recent study has analysed the photodissociation of H_{2}^{+} within the framework of this theory [43] and confirmed the results of previous studies [44]. Two main points have been revealed. First, the dissociation process is correctly described by collecting in the evolutive target spaces the resonance quasi-states characterized by Gamov boundary conditions, thus neglecting the diffusion states which possess small lifetimes (see fig. 7). Second, it is possible to repeat fast calculations of dissociation probabilities for varying pulse shapes if some data has been stored during the initial construction of the target spaces: in particular, the eigenvalues of the vectors constituting the basis set of S^{(j)}: E_{γ}^{j}, the overlap matrices of basis sets of consecutive target spaces < Φ_{γ_{j+1}}|Φ_{β_{j}} > and the asymptotic amplitudes of the basis vectors Φ_{j}(r_{∞}). This allows a detailed analysis of the influence of the pulse shape on the dissociation probabilities and on the kinetic energy spectrum of the fragments [28].

5. Conclusion. Quantum dynamical problems involving coupled molecular and photon continua require the use of discrete representations which require sets of square integrable functions of large dimension. Various strategies have been used to reduce the dimension of these basis sets (cf the article by Bacic in this volume). In most cases, this dimension reduction is only a first step and must be accompanied by other simplification procedures at the level of the dynamics. The Bloch wave operator theory presented in this chapter proposes a procedure based on the selection of an active space within which the dynamics is confined, together with a projection of the dynamical equations into this space. In recent years several formulations of dynamical theory ivolving active subspaces have been given ; the selection of such a space is often made in the context of artificial intelligence methods. However, the theory of Bloch wave operators and its numerical implementation offers several specific advantages for handling problems of reactive dynamics:

• The selection of the active space and the construction of the projected dynamical equations are both based on the same central concept, that of the wave operator.

• The theory is equally applicable to treatments of the Hamiltonian operator which lead to either real or complex matrices; this flexibility is essential in the treatment of molecular resonances using discrete finite bases.

• The eigenvectors which are associated with the active space, and which contain all the relevant dynamical information, can be obtained from the wave operator calculation.

• The fundamental non-linear equations which constitute the heart of the method can be solved by several iterative treatments in which the basic computational operation, requiring the largest portion of CPU time at each iteration, is the formation of the matrix-column product HΨ, where H is the fixed Hamiltonian matrix and Ψ is a column which changes at each iteration. This feature of the calculations is obviously suitable for parallel computations.

• The iterative nature of the calculations permits the introduction of adjustable parameters into the Hamiltonian and the study of the parameter dependence of results without the need to repeat ab initio the whole of the calculation. This feature has been exploited in the stabilization approach to molecular resonances, where the variable parameter is the amplitude of the optical potential term, and in the study of laser-induced photodissociation, where the variable parameter is the electric field amplitude.

• From a didactic point of view, the wave operator approach provides a deeper understanding of reactive phenomena ; effects which at first sight involve several continua often reduces to the study of a few discrete resonance states.

The second essential point about our formulation is the use of the general concept of a time-dependent wave operator which has the stationary wave operator of Bloch, widely utilized in quantum chemistry, as its adiabatic limit. Application of the fully time-dependent wave operator are as yet little developed, but we have high hopes of its utility. Two formal features of the operator are already noteworthy. The concept allows a new interpretation of the adiabatic theorem. In the presentation of Messiah [26] the theorem concerns only the strict adiabatic limit for spaces which are purely non degenerate. The wave operator equations permit the study of more general evolution processes characterized by two time scales, one relating to a rapid internal evolution of the active space, the other relating to quasi-adiabatic exchanges between the active space and its complementary space. Equations (21) and (22), which define the timedependent wave operator, transform the time-dependent dynamical equations into stationary equations in which the original Hamiltonian H(t) is transformed into a Hamiltonian H(t) − ih̄∂/∂t of Floquet type. This approach gives a changed formal status to the time variable and suggests some new routes for the treatment of the Schrödinger equation.

REFERENCES

[1] R. Kosloff,Time dependent quantum mechanical methods for molecular dynamics, J. Phys. Chem., 92 (1988), pp. 2087-2100.

[2] G.D. Billing, Classical path approach to inelastic and reactive scattering, in Super computer Algorithms for Reactivity, Dynamics and Kinetics of small Molecules, Reidel Publ. Co. North-Holland, 1989

[3] G.D. Billing, The dynamics of molecule-surface interaction, Comp. Phys. Rep., 6 (1990), pp. 383-450

[4] R. Heather, An asymptotic wavefunction splitting procedure for propagating spatially extended wavefunction : application to intense field photodissociation of H_{2}^{+}, Comp. Phys. Comm, 63 (1991), pp. 446-459.

[5] D. Neuhauser, D.M. Baer, R.S. Judson and D.J. Jouri, The application of timedependent wavepacket methods to reactive scattering, Comp. Phys. Comm, 63 (1991), pp. 460-481.

[6] D.J. Kouri and R.C. Mowrey, Close Coupling wave packet formalism for gas phase non reactive atom-diatom collisions, J. Chem. Phys, 86 (1987), pp. 2087-2094

[7] C. Leforestier, R.H. Bisseling, C. Cerjan, M.D Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.D. Meyer, N. Lipkin, O. Roncero and R. Kosloff, A comparison of different propagation schemes for the time dependent Schrödinger equation, J. Comp. Phys, 94 (1991), pp. 59-80

[8] T.J. Park and J.C. Light, Unitary quantum time evolution by iterative Lanczos reduction, J. Chem. Phys, 85 (1986), pp. 5870-5876

[9] Y. Sun, D.J. Kouri, D.W. Schwenke, Time-dependent wave packet algorithm for inelastic molecule-molecule scattering, Computer Phys. Commun, 63 (1991), pp. 51-62

[10] R. Kosloff and D. Kosloff, Absorbing boundaries of wave propagation problems, J. Comp. Phys., 63 (1986), pp. 363-376

[11] J.H. Shirley, Solution of the Schrödinger equation with a Hamiltonian periodic in time, Phys. Rev., B138 (1965), pp. 979-987

[12] U. Peskin and N. Moiseyev, The solution of the time dependent Schrödinger equation by the (t, t′) method. Theory, computational algorithm and applications, J. Chem. Phys., 99 (1993), pp. 4590-4596

[13] B.R. Junker, Recent computational developments in the use of complex scaling in resonance phenomena, Advan. At. Mol. Phys., 18 (1982), pp. 207-263

[14] R. Lefebvre,Siegert quantization, complex rotation and molecular resonances, J. of Phys. Chem, 88 (1984), pp. 4839-4844

[15] W. Reinhardt, Complex coordinates in the theory of atomic and molecular structure and dynamics, Ann. Rev. Phys. Chem., 33 (1982), pp. 223-255

[16] G. Jolicard and J.Périé, The theory of stationary and time-dependent wave operators and their applications in photochemistry, Advan. Multiphoton Process and Spect., Vol. 8 (1994), pp. 1-126

[17] G. Jolicard and E.J. Austin, Optical potential method of calculating resonance energies and widths, Chem. Phys., 103 (1986), pp. 295-302

[18] M. Garcia Sucre and R. Lefebvre, Localization of the wave functions of resonance states through the use of an optical potential, J. Chem. Phys., 85 (1986), pp. 4753-4754

[19] N. Rom, N. Lipkin and N. Moiseyev, Optical potential by the complex coordinate method, Chem. Phys., 151 (1991), p. 199-204

[20] G.G. Balint Kurti, R.N. Dixon and C.C. Marslon, Time-dependent quantum dynamics of molecular photoframentation processes, J. Chem. Soc. Faraday Trans., 86 (1990), pp. 1741-1749

[21] J.Périé and G. Jolicard, A spectral analysis method using the optical potential as a wave absorber, J. Phys. B., 26 (1993), pp. 4491-4502

[22] U. Peskin, N. Moiseyev and R. Lefebvre, Partial widths by asymptotic analysis of the complex scaled resonance wave function, J. Chem. Phys., 92 (1990), pp. 2902-2909

[23] J.Périé, G. Jolicard and J.P. Killingbeck, An iterative Bloch approach to the resonance state problem, J. Chem. Phys, 98 (1993), pp. 6344-6351

[24] J.A.R. Coope, Partitioning without implicit eigenvalues, Mol. Phys, 18 (1970), pp. 571-574

[25] J. Killingbeck and G. Jolicard, Perturbation methods for the matrix eigenproblem, J. Phys. A. Math. Gen., 25 (1992), pp. 6455-6459

[26] A. Messiah, Mecanique Quantique Part 2, Dunod, Paris 1964

[27] G. Jolicard, A recursive solution to the time dependent Schrödinger equation using a generalized quasidegenerate Bloch formalism, J. Chem. Phys., 90 (1989), pp. 2320-2327

[28] G. Jolicard, The solution of the time dependent Scrhödinger equation by the timedependent wave operator method, Phys. Rev. A to be published

[29] Ph. Durand, Direct determination of effective Hamiltonians by wave operator methods. I. General formalism, Phys. Rev., A28 (1983), pp. 3184-3192

[30] C. Bloch, sur la théorie des perturbations des états liés, Nucl. Phys., 6 (1958), pp. 329-349

[31] J. Des Cloizeaux, Extension d’une formule de Lagrange à des problèmes de valeurs propres, Nucl. Phys., 20 (1960), pp.321-346

[32] J.P. Malrieu, Ph. Durand and J.P. Daudey, Intermediate Hamiltonians as a new class of effective Hamiltonians, J. Phys. A Math. Gen., 18 (1985), pp. 809-826

[33] J.P. Killingbeck, Microcomputer algorithms. Action from algebra, Adam Hilger, Bristol 1991

[34] G. Jolicard and A. Grosjean, on the validity of partitioning techniques for molecular infrared-multiphoton excitation, J. Chem. Phys., 95 (1991), pp. 1920-1927

[35] R.E. Wyatt, C. Iung and C. Leforestier, Quantum dynamics of overtone relaxation in benzene I. 5 and 9 modes for relaxation from (CH(v=3), J. Chem. Phys., 97 (1992), pp. 3458-3476

[36] J. Chang and R.E. Wyatt, Preselecting paths for multiphoton dynamics using artificial intelligence, J. Chem. Phys, 85 (1986), PP. 1826-1839

[37] G. Jolicard, J.P. Killingbeck, Ph. Durand and J.L. Heully, A wave operator description of molecular photodissociation processes using the Floquet formalism, J. Chem. Phys., 100 (1994), pp. 325-333.

[38] J;C. Barthelat and Ph. Durand, Theory of effective Hamiltonians and applications, in computational Chemistry structures interactions and reactivity Part A, Elsevier, Amsterdam, 1992

[39] G. Jolicard, Effective Hamiltonian theory : an intermediate representation method for the wave operator calculation, Chem. Phys., 115 (1987), pp. 57-68

[40] G. Jolicard and E.J. Austin, On the recursive integration of stationary and timedependent Bloch wave operator equations, Chem. Phys. Letters, 180 (1991), pp. 503-508

[41] T. Tung Nguyen Dang, Adiabatic time evolution of atoms and molecules in intense radiation fields, J. Chem. Phys., 90 (1989), pp. 2657-2665

[42] S.I. Chu, Generalized Floquet theoretical approaches to intense-field multiphoton and nonlinear optical processes, Advan. Chem. Phys, 73 (1989), pp. 739-798

[43] G. Jolicard and G.D. Billing, The extension of wave operator-Floquet formalism to molecular photodissociation processes with short laser pulses, J. Chem. Phys. (submitted for publication)

[44] G. Giusti-Suzor, X. He, O. Atabek and F.H. Mies , Above threeshold-dissociation of H_{2}^{+} in intense laser fields, Phys. Rev. Lett., 64 (1990), pp. 515-518 ; G. Jolicard and O. Atabek, Above threshold dissociation dynamics of H_{2}^{+} with short intense laser pulses, Phys. Rev., A46 (1992), pp. 5845-5855