Improved L2-stabilization theory to compute resonances under multichannel conditions

2 February 1996

CHEMICAL

PHYSICS LETTERS

ELSEVIER

Chemical Physics Letters 249 (1996) 237-243

Improved L2-stabilization theory to compute resonances under multichannel conditions R.F. Salzgeber, U. Manthe, Th. Weiss, Ch. Schlier Fakultfit ffir Physik, Universiti~t Freiburg, Hermann-Herder Strafle 3, D-70104 Freiburg, Germany

Received 25 July 1995; in final form 6 November 1995

Abstract

The stabilization theory for the calculation of resonance energies and widths in the version of Mandelshtam, Taylor, Ryaboy and Moiseyev has been applied to collinear ABA molecules in which two Morse potentials are kinetically coupled. This system has about 40 resonances in 10 channels. The results are compared with converged computations by the complex scaling method. Subtracting the background phase of the uncoupled problem (but including the full potential) improves the stabilization calculation appreciably, allowing easier identification of the resonances and better fits of the scattering phase to the arctan-sigmoids, from which the resonance position in the complex plane can be computed.

1, I n t r o d u c t i o n Much effort has been put into the computation of molecular vibrational resonances. This interest comes from new experimental techniques for the measurement of lifetimes of single predissociating states [1,2], from new possibilities for recording precise photodissociation spectra [3] and recent work on the dynamics underlying unimolecular decay [4-9]. Several methods compete to be the best for such compurations; we mention only wavepacket propagation [10], several variants of complex scaling [11-15] and the stabilization method [16]. The latter was recently revived [17] and shown to be a general method, allowing the computation of any observable which depends on the resonances [14,17-21]. Whereas the literature has many papers applying different methods to systems having just one or two channels and only a few resonances (notable exceptions are the Y-matrix method of Refs. [22,23] and resonances in Coulomb problems [24,25]), we feel

that it is also necessary to implement these methods in a way which makes it possible to investigate rich spectra of metastable states with all their complexity, e.g. multichannel decay, overlapping of resonances and threshold effects. Insight into the behavior of such large sets can then help to understand the statistical behavior of resonances [4,7-9,26,27] or the correlation of resonance parameters to the system's classical dynamics [28,29]. Different from bound states, resonances are always mixed with (representations of) continuum states. To achieve the computation of large sets of them one should not only have an effective method of computing the data set containing them, but also be able to automatically extract them from a much larger continuum background. We feel that this task has largely been neglected, since it demands replacing eye and hand by difficult heuristic programming. We have more or less solved this problem for the complex scaling method, and are working on it for the stabilization method. In this Letter we describe

0009-2614/96/$12.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0009-2614(95)01390-3

238

R.F. Salzgeber et a l . / Chemical Physics Letters 249 (1996) 237-243

an improvement of the latter, which we think will aid in the automatic discrimination of resonances from the continuum background. We start from the recent version [20] of the stabilization method, which uses the eigenphase sum ~7(E) to compute the resonances. ~lere, the authors apply the method to the full Hamiltonian yielding a phase shift which diverges for L ~ oo. To avoid this, one must subtract a background phase, which is generally taken as that of the free particle Hamiltonian. It could only be neglected because in any practical application L remains finite, In this Letter, the idea is to take seriously the subtraction of the background phase, and to choose it as the true difference between the coupled and uncoupled Hamiltonians. This reduces the accumulated background phase at the highest resonance of the system more than tenfold, and has two consequences: broad resonances stand out against the background more clearly, and the fit of the arctanjumps in r/(E), which identify the resonances, becomes much easier, more accurate and safer. Thus some resonances that cannot be identified in the straightforward unsubtracted eigenphase calculation can be seen without problems after background subtraction has been implemented. In Section 2 we give a brief summary of the two versions of stabilization theory discussed here. In Section 3 we describe our application and analyze the errors. In Section 4 a short discussion follows,

2. The stabilization method The key idea of stabilization theory is that resonances are localized near the potential well. Putting the whole system in a box of size L, much larger than the range of this well, results in a new system which should cover all the localized states but does not have an infinitely dense continuum. Once the finite basis is converged, and the range L of this basis is much larger than the well region, nothing dramatic should happen to the resonant states as L is further increased. All other states (i.e. the representation of the continuum) must change with L. To study resonance parameters with the modified stabilization theory [20] one needs the set of energy eigenvalues versus box size Ej(L). Therefore, one

must diagonalize a real-valued Hamiltonian in a given basis for many box sizes L i. Since the energy eigenvalues for eigenstates in an empty box behave as Ej(L) ~ 1/L 2, they must cross an existing resonance energy E r at a certain box size L. It can be shown that near a resonance the slope of the eigenvalues Ej(L) is smaller than far away from it, which means that pieces of the Ej(L) curves 'stabilize' around E r. In Refs. [18,20] Mandelshtam et al. have demonstrated how to extract the resonance parameters from these curves. One starts from the change in the density of states caused by the interaction

Ap(E)=Tr[$(E-H)-

B ( E - H0)],

(1)

where H 0 is the Hamiltonian without interaction. For an isolated resonance, the density of states gives directly the resonance parameters E r and F. Applying Eq. (1) to the box Hamiltonian and averaging over boxes extending from Lmi n tO Lmin + A L = Lmax yields [18] the approximation (1 j~. ] ~ [ - t Ap(E) = lim a I.-, ~ ~ ~ L)= e

-/90( E, L m i n A L ) ) .

(2)

where P0 is the averaged (over boxes extending from Lmi n tO Lmax) density of states of the non-interacting Hamiltonian. Lj(E) are the inverted Ej(L) curves. The index j numbers the curves resulting from different eigenvalues, and sums over all j for which Ej(L) = E with L ~ [ Lmin, Lmax ]. To avoid the differentiation in Eq. (2), one can [20] also compute the eigenphase sum ~/(E) of the system, which is related to the density by [30] r/(E) =

~r["Ap(E ') dE'.

(3)

"o Now, provided we are able to compute r/(E) from our data and assuming that the phase shift near a resonance with energy E r and width F is described by an integrated Breit-Wigner formula with a nonresonant background ( E -- E r / 7 ( E ) = '0background(E) + arctan / ~ , ], (4)

R.F. Salzgeber et al./ Chemical Physics Letters 249 (1996) 237-243

one can fit this Ansatz to the approximated eigenphase and obtain the resonance parameters from the fit. In Ref. [18] Mandelshtam et al. showed that the eigenphase sum can be approximated by: "q(E)=

239

In the case of the kinetically coupled double Morse oscillator, which describes a collinear A-B-A molecule, this gives 1 H 0 = .~_~( p~ + p2) + VMor~e( r l ) + VM.... ( r 2 ) ,

lim (~b(E, Lmin+AL)

AL---~~

-- (DO( E , tmi n --I--A t ) ) ,

(5)

HI = - - - P l P 2 ' ],z

(7)

with

where

(]>( E , Cmi n "Jr a n )

WMorse = D ( e -t~(r-r°) - 1) 2.

1 =w-~]~Lj(E) +N(E, Lmin). (6) J Here qS(E, Lmi n "k- A L ) is the phase shift caused by the full Hamiltonian and ~b0(E, Lmi n + A L) the phase shift caused by the Hamiltonian without interaction, N(E) is the number of states for which Ej(Lmin) < E and j numbers the eigenvalues as in Eq. (2). Note that qb(E, Lmin, AL) and q~0(E, Lmin, AL) diverge for AL ~ ~, but r/(E) remains finite, and is

In the calculations we used [31] / z = 30/31 amu, D = 0 . 7 6 eV, 13= 1 ao I and r 0 = 1.4 a 0, which yields 10 bound states in each single bond, and took ~ = 1/7, 0.35 and 0.5 as examples. All matrix elements have been calculated in a basis consisting of symmetrized products of the functions / 2

independent of the underlying basis set. Therefore, "r/(E) is a quantity which is much easier to handle. In Ref. [20], arguing that po(E) 'cannot yield any reasonable information', ~b0(E) was incorporated into the background phase "qbackground(E), defined as the non-resonant contribution to ~(E). This leads to a purely numerical 'background', which increases if one tries to improve the calculation by a larger basis or larger A L. In this Letter we avoid this effect, and reduce the non-resonant background to its physical minimum.

3. Subtraction of the background phase shift and accuracy of the method To compute the unique quantity rl(E) in Eq. (3), we must compute ~b(E, Lmi n + AL) and the background phase shift ~b0(E, Lmin + AL). The latter is determined by the choice of the unperturbed Hamiltonian. Minimizing r/b~kg~o~,a profoundly simplifies the identification of resonances employing Eq. (4). Therefore, in this Letter we compute ~b0(E, Lmin, AL) as the phase of the uncoupled rather than the free particle Hamiltonian.

~.(kr) =

V

(n + 1)(n + 2) 2rk3/eL2n(2kr)'

(8)

(9)

where the L2.(x) are the generalized Laguerre functions, defined by 1 d" - e X x - " d x" (e-Xx"+~). L~ ( x ) = -n!

(10)

The momentum-, coupling- and coordinate-matrices are known analytically, and the potential matrix is conveniently computed with the standard DVR method [32,33]. We used 1830 basis functions (corresponding to 60 basis functions in each bond) and varied their range not by putting the Hamiltonian in a box as suggested above, but by varying k in Eq. (9) in the interval 3 ~< k ~< 10. This causes some inconvenience as explained in the next paragraph, but we wanted to use the same basis as in the complex scaling calculations. Since all eigenfunctions necessarily vanish where the largest basis function ~6Ü(kr) in any degree of freedom approaches zero, one may associate their range to the boxes required by the stabilization theory. Thus, the above values of k correspond to box sizes of 1.14 ~
240

R.F. Salzgeber et aL / Chemical Physics Letters 249 (1996) 237-243

continuum, i.e. it is not possible to converge the non-resonant eigenvalues for a given basis q t ( k r ) to the eigenvalues ° f a b°x ° f size x~aX" In ° a t case' if one increases the basis size N keeping k constant, the eigenvalues get denser rather than converging to a constant number as in a box Hamiltonian. It follows that the density of states and the total phase shift at a given energy depend on the choice of the basis, meaning that they get larger and larger if one increases the number of basis functions. Thus the identification of a sudden jump of ~ in the total phase shift 4,(E) becomes more and more difficult. However, if we subtract the uncoupled phase 4,o to get the 'pure' resonance phase 77(E), this undesired effect vanishes almost completely, and t h e r e s o n a n c e s pop out clearly. Since our goal is to study the method, and we did not want it to be polluted by additional numerical problems, (e.g. high order interpolation of E(L)), we used 700 different box sizes (i.e. 1400 diagonalizations) to get the stabilization data. A further small modification is necessary because the phase of the uncoupled system, defined by H 0, contains the phase shifts due to its bound states. For convenience one may call these states 'resonances' of 6 = 0. Their lifetime is infinite and their width is zero, therefore, they appear as discontinuous jumps in the 4,0(E) curve. These resonances must not appear in T/(E), since they have no corresponding features in the total phase. Their elimination is easy, since their phase shift is exactly "rr and their positions are known analytically. We eliminated these gaps by numerically pasting the pieces of 4,0(E, Lmin, AL) together. In Fig. 1 we show a part of.the unmodified total phase 4 ' ( E ) and the 'pure' resonance phase r/(E). In both cases the same energy range covering ten resonances and one threshold is shown. The arctan shape of the resonances is clearly visible for all resonances in the 'pure' phase, but can be much harder and in some cases barely be seen in the total phase (viz. the resonance at -- 1 eV). In Fig. 1 only a part of the data, which actually covers all energies up to the three-panicle dissociation threshold, has been shown. The error of the method compared with complex scaling data (computed with a basis set with N = 1830 and k = 6.9 ao 1) is shown in Fig. 2. Here the relative error for each resonance energy and

7° I . . . . . . . . . ~. 60 baekg. . . . d i n c l u d e d ~5° I i 40 ~ 30 r~ . . . . . . . . . . . . . . . . . . . . . 085 09

12 . . . . . . . . . . . . .

.

~ : ~ ; /if ~ ' ~' ' ' !

J

~ I

i

0.95

energy [eV]

~

1o5

.....

~. • ,~.

x0 ~~ ~ ~ } 8~ ! b.ekg.... d ~ b t r . ~ t ~ ? ~ ~ i ~u° ~ '2~ - " ~ ~"~ 0 F" . . . . . . . . . . . . . . . . . . . . . . 0.B5 0.9 0.95 I .... gy [.v]

"

1 05

Fig. 1. Comparison of the phase shifts of the standard stabilization

theory

(upper plot) with the modified one, i.e. ~b0 subtracted

(lower plot), for 6 = 1 / 7 . Note the different scales for the phase

shift. Resonance positions are indicated by arrows, the dashed line marks a threshold. Only a part of the energy range is shown.

width is plotted versus the number of the resonance (1 = lowest; 38 =highest). The error increases slightly with energy since there is still a residual coupling between the resonances, and more and more thresholds are present. The corresponding effects accumulate over the whole energy range to an additional phase shift of about 8"tr in the case of 6 = 0.5, and about ~ in the other two cases. Since the convergence with respect to the matrix and box size, and the interpolation of Ej(L) is good to 10 -7, the

~ lo-~°'"~".-/t.~At,.~t o.ol

r

,

', ~ ' o i ~ i i o n "~

.

. .

.

.

.

.

.

.

.

.

,

,

, - ., . , . ,.I , ~. . . .,

'" . . . ". . . . .- -

~ ~ 1o-, "~ 10 ~ 0

10

20

30

R........

1o~ ~. . . .~. . ./. . .A. . .~. . . I. . .~. . . .' . .~. . . I. . .I. . .I. . .I. . i ~I, " / k - ~ / ~ / \ l t v k / ~ \ l ' k J \ / ~

]

.~ 01 001 o

1o

zo R ........

~o #

Fig. 2. Relative error of resonance parameters for 6 = 1 / 7 . Triangles are data with subtraction of q~0, and circles are data without that. Note that at resonance No. 25 the fit for the phase with

background exceeds the range of the graph because it did not converge to a meaningful value.

241

R.F. Salzgeber et a l . / Chemical Physics Letters 249 (1996) 237-243

1, , ~ -

. . . . . . . . . . . . . ....

: ' : ~

135 13 125 f~

.... ----

.... -.

.~ ' ~

__I --~

1.2

~

115

.

,- . . . . . . . . . . 0.86 .... gy [eVl

0.85

0.87

range Lmax-Lmin has been increased with respect to the lowest one. The fitting procedure in Eq. (4), which is intricate due to the sudden changes caused by the generally small widths, reacts sensitively to these variations, causing the observed errors, whose mean values are shown for three values of t~ in Table 1. In this table we show results for three molecular models having a weak, an intermediate and a strong coupling parameter (t5 = 1 / 7 , 0.35 and 0.5, respectively). As can be seen in Fig. 2 resonance No. 25 could not be fitted to the t5 = 1 / 7 phase shift with ~b(E) alone. In other words this resonance would not

Fig. 3. An enlarged part of a typical arctan profile is shown here to illustratethe effects of the finite value of Lrain and Lmax in Eq. (2). The lowest curve was stabilized in the smallest range 4.5 < x~'ax < 7.1 ao, the middle and upper curve in larger ranges of 4.5 < x~ax < 12.6 a o and 1.2 < x~ax < 7.1 a o, respectively,

have been found if one had not known its presence. The improvement gained by the subtraction of q50 is approximately a factor of 1.2 in the energies of all cases studied, and in the widths for the weak coupiing case. There is no improvement for the average

dominating error is the result of finite box sizes (Lrnin , Lmax) in Eq. (6). This approximation causes a non-differentiable point in th(E), every time a new curve enters or leaves the box at Lmax or Lrnin, respectively. Note that this effect will be basis dependent if one does not use an explicit box Hamiltonian. One can put Lmin = 0 so that no level vanishes at the left side of the stabilization graph, but a finite Lmax will remain, which dominates the error of the method. The effect can be seen in Fig. 3 where a small enlarged window (of the actual points) of r/(E) is shown. The sequence of dots looks piecewise smooth and differentiable, but often sharp, non-differentiable bends occur. To show the effect of smaller or bigger box sizes, the corresponding plots are also shown. It can be seen that the vertical variation (the distance between the two envelopes of the curve) is smaller in the upper curves, where the

width of the intermediate and the strong coupling case. This can be understood if one realizes that some gain of the new method is counteracted by the additional noise from the numerical subtraction of two noisy sets. It is remarkable, though, that the accuracy of the method seems to be rather constant over a range of coupling parameters. Computations with different basis sizes n = 820, 1275 and 1830 also show that the resonance parameters do not vary significantly within the above errors, a further implication of the fact that the finite box size error is much larger than that from other error sources e.g. the limited basis size. The overall result is that the modified stabilization method with background reduction by the subtraction of th0 is able to find the same resonances as in the complex scaling calculations of this rather general system, even the ones with a large width, an accuracy of = 1 meV in the energy. This is an improvement to the previous version of the stabiliza-

Table 1 Mean relative error of resonanceparameters 8= 1/7 a with ~bo

(AE/E)

0.452×

( A I'll" )

1.093

10 - 3

a ResonanceNo. 25 not included.

tS= 0.35

~5= 0.5

without 4'0

with qbo

without ~bo

with tho

without ~bo

0 . 3 8 2 X 10 - 3

0 . 6 7 0 X 10 - 3

0 . 5 8 0 X 10 - 3

0 . 3 9 4 X 10 - 3

0 . 3 2 1 X 10 - 3

0.846

0.507

0.603

1.512

1.904

242

R.F. Salzgeber et al./ Chemical Physics Letters 249 (1996) 237-243

tion method, since there some resonances will easily escape detection,

prize granted to H.S. Taylor and J. Briggs. RFS thanks H.S. Taylor for his hospitality and V. Mandelshtam for helpful discussions.

4. Discussion References

We have shown that it is possible to reduce a complicated stabilization graph to a n a l m o s t p u r e resonance phase shift, once one can separate the interaction term in a given Hamiltonian. If o n e c a n

[1] W.F. Polik, D.R. Guyer and C.B. Moore, J. Chem. Phys. 92 (1990)3453. [2] S.I. Ionov, C.J. Brucker, C. Jaques and Y. Chen, J. Chem.

separate a given potential into V(r t) + V(r 2) + V ( r 1, r 2) with not too big V(r~, r2) , this method will also work for a fitted ab initio potential surface. The major improvement of the background subtraction is that the resulting arctan forms of the phase

[3] R. Schinke, Photodissociation dynamics (Cambridge Univ. Press, Cambridge, 1993). [4] W.H. Miller, R. Hemandez and C.B. Moore, J. Chem. Phys. 93 (1990) 5657. [5] M. Desouter-Lecomte and F. Remacle, Chem. Phys. 164

shift turn out much clearer and are, therefore, easier to identify and to fit. An extraction of the exact number of resonances in a given energy range at no further computational cost s e e m s possible, but must still be investigated

[6] M. Desouter-Leeomte and F. Culot, J. Chem. Phys. 98 (1993) 7819. [7] W.H. Miller, R. Hernandez and C.B. Moore, J. Chem. Phys. 99 (1993) 10078. [8] K. Someda, N.H. and F.H. Mies, J. Chem. Phys. 187

further. The main source of errors in the modified stabilization method is the primitive approximation for the phase shift ~b(E), which depends on the largest b o x size Lmax, since the number of eigenvalues below a given energy depends on it. If another kind

(1994) 195. [9] U. Peskin, W.H. Miller and H. Reisler, J. Chem. Phys. 102 (1995) 8874. [101 A. Isele, C. Meier, V. Engel, N. Fahrer and Ch. Schlier, J. Chem. Phys. 101 (1994) 5919. [11] w.P. Reinhardt, Ann. Rev. Phys. Chem. 33 (1982) 223. [12] Y.K. Ho, Phys. Rept. 99 (1983) 1. [13] N. Lipkin, N. Moiseyev and C. Leforestier, J. Chem. Phys.

of stabilization procedure is applied, e.g. scaling of the basis functions, this error may also depend on the functional form of the basis, which determines this number of eigenvalues. A n o t h e r advantage of this m e t h o d is the possibility of automating the process of finding a large number of resonances (this would not be easy, if possible at all, from the unsubtracted phase shift). Such an automated search for resonances has not yet been performed in connection with stabilization the-

ory, though it has been used with complex scaling in our group [34], but is quite important if one is interested in the statistical behavior of the r e s o n a n c e s in a system. Acknowledgement This

work

was

supported

by

the

Deutsche

Forschungsgemeinschaft through the SFB 276 " K o r relierte Dynamik hochangeregter atomarer und molekularer Systeme" and by the G e r h a r d H e s s

Phys. 99 (1993) 3420.

(1992) 11.

98 (1993)1888.

[14] V. Ryaboy, N. Moiseyev, V.A. Mandelshtam and H.S. Taylor, J. Chem. Phys. 101 (1994) 5677. [15] C. Leforestier, K. Yamashita and N. Moiseyev, submitted for publication. [16] A.U. Hazi and H.S. Taylor, Phys. Rev. A 1 (1970) 1109. [17] V.A. Mandelshtam, T.R. Ravuri and H.S. Taylor, Phys. Rev. A 48 (1993)818. [18] V.A. Mandelshtam, T.R. Ravuri and H.S. Taylor, Phys. Rev. Letters 70 (1993) 1932. [19] V.A. Mandelshtam and H.S. Taylor, J. Chem. Phys. 99 (1994) 222. [201 V.A. Mandelshtam, H.S. Taylor, V. Ryaboy and N. Moi-

seyev, Phys. Rev. A 50 (1994)2764. [21] V.A. Mandelshtam, T.R. Ravuri and H.S. Taylor, J. Chem. Phys. 101 (1994)8792. [22] A. Dobbyn, M. Stumpf, H.M. Keller, W. Hase and R. Schinke, J. Chem. Phys. 102 (1995) 5867. [23] H.J. Wemer, C. Bauer, P. Rosmus, H.M. Keller, M. Stumpf and R. Schinke, J. Chem. Phys. 102 (1995) 3593. [24] K. Richter, J.S. Briggs, D. Wintgen and E.A. Solov'ev, J. Phys. B 25 (1992) 3929; note that the scaling properties of the Coulomb potential simplify the numerics appreciably. [25] K. Richter and D. Wintgen, J. Phys. B 26 (1993) 3719. [26] R.D. Levine, Ber. Bunsenges. Physik. Chem. 99 (1988) 222. [27] R.D. Levine, Advan. Chem. Phys. 70 (1988) 53.

R.F. Salzgeber et al. / Chemical Physics Letters 249 (1996) 237-243 [28] T. Terasaka and T. Matsushita, Phys. Rev. A 32 (1985) 538. [29] J. Burghardt and P. Gaspard, J. Chem. Phys. 100 (1994) 6395. [30] R.D. Levine, Quantum mechanics of molecular rate processes (Clarendon Press, Oxford, 1969) p. 102. [31] R.H. Bisseling, R. Kosloff, J. Manz and F. Mrugala, J. Chem. Phys. 86 (1987) 2626.

243

[32] D.O. Harris, G.G. Engerholm and W.D. Gwinn, J. Chem. Phys. 43(1965)1515. [33] A.S. Dickinson and P.R. Certain, J. Chem. Phys. 49 (1968) 4208. [34] S. Dallwig, I. Weese, T. Weiss and C. Schlier, to be published.