15966
J. Phys. Chem. B 2008, 112, 15966–15972
ARTICLES Semiclassical Nonadiabatic Surface-hopping Wave Function Expansion at Low Energies: Hops in the Forbidden Region† Michael F. Herman Department of Chemistry, Tulane UniVersity, New Orleans, Louisiana 70118 ReceiVed: June 4, 2008
The accuracy of a semiclassical surface-hopping expansion of the time-independent wave function for problems in which the nonadiabatic coupling is peaked in the classically forbidden regions is studied numerically for a one-dimensional curve-crossing problem. This surface-hopping expansion has recently been shown to satisfy the Schrodinger equation to all orders in p and all orders in the nonadiabatic coupling. It has also been found to provide very accurate transition probabilities for problems in which the crossing points of the diabatic energy surfaces are classically allowed. In the numerical study reported here, transition probabilities are evaluated for energies well below the crossing point energy. It is found that the expansion provides accurate results for transition probabilities as small as 10-11. I. Introduction Semiclassical methods often offer attractive alternatives to fully quantum approaches for the calculation of collision processes in atomic and molecular systems.1-14 A variety of semiclassical techniques have been proposed for problems in which more than one adiabatic electronic state plays a significant role.15-47 It has recently been shown33 that a semiclassical surface-hopping expansion for multistate, multidimensional nonadiabatic problems can be developed into a formally exact solution to the time-independent Schrodinger equation (TISE) for the motion of the nuclei by including single-surface correction terms that allow for energy conserving changes in the direction of the momentum along classical trajectories. This expansion for the wave function (Ψ) is formally exact in the sense that all terms have been shown to cancel when it is inserted into (H - E)Ψ, giving zero for the answer. However, the expansion, each term of which has a primitive semiclassical prefactor and phase function, has the usual semiclassical divergence at classical turning points and caustics. As a result, it is not convergent near these points. Nonetheless, the fact that this expansion formally satisfies the TISE demonstrates that it includes all appropriate hopping and momentum reversal events and that the phases for all trajectories are correctly calculated. Furthermore, the semiclassical version of the surface-hopping expansion, which ignores the single-surface correction terms but keeps the hopping terms, has been shown to provide very accurate transition probabilities, even for problems where there is considerable interference between different avoided crossing regions.23,28,33 Previous calculations have only considered transitions where the avoided crossings of the adiabatic electronic energy surfaces are in the classically allowed region.23,28,33 In this work, the application of the formalism to one-dimensional (1D) problems in which the avoided crossing is in the forbidden region is †
Part of the “Karl Freed Festschrift”.
considered. Accurate results are obtained for transition probabilities down to very small values, corresponding to highly forbidden processes. The extension of the trajectories into the forbidden regions and the inclusion of hops in these regions are required to obtain accurate results. It is found that the small transition probability results from significant cancelation between contributions from hops in the allowed region and from hops in the forbidden region. This paper is organized as follows. The surface-hopping expansion is presented for the 1D case in Section II. A 1D model curve-crossing problem is presented in Section III, and the semiclassical surface-hopping transition probabilities are compared with the results from exact quantum calculations. These results are discussed, and an analysis in terms of contour integrations in the complex plane is presented in Section IV. II. Theory This work numerically studies the behavior of the semiclassical surface-hopping expansion of the time-independent wave function23,33 for 1D problems with turning points. Let the incoming flux be on adiabatic surface W1(x). At each point along a trajectory contributing to the surface-hopping wave function, this trajectory can continue to evolve classically or it can undergo a momentum reversal or it can undergo an energyconserving hop to the other adiabatic surface.33 The energy determines the magnitude of the momentum after the hop, but not its sign. If the momentum has the same sign before and after the hop, this hop is referred to as a transmission or T-type hop. If the sign of the momentum changes, it is referred to as a reflection or R-type hop. After the hop, the trajectory continues to evolve classically until it undergoes another hop and/or momentum reversal. Momentum reversals without a hop correct for the use of semiclassical expressions for the contribution to the wave functions between hops.23,33,48 These single-surface reflection terms are neglected in this work, as they have been in previous numerical calculations,23,27,33 resulting in a semiclas-
10.1021/jp804937q CCC: $40.75 2008 American Chemical Society Published on Web 10/31/2008
Hops in the Forbidden Region
J. Phys. Chem. B, Vol. 112, No. 50, 2008 15967
sical approximation for the multistate wave function. Each trajectory from initial point x0 to final point x with any number of hops gives rise to a contribution to the wave function at x. This contribution contains a prefactor A, a phase factor of the form exp(iS/p), and an amplitude for each hop. For instance, the contribution to the wave function on W2 from a trajectory that starts at x0 on surface W1 with energy E, has a T-type hop to surface W2 at x1, and then continues to x, has the following form;
Ψ(x0, x1, x) )
p1(x0) i τ (x )exp p2(x) 21 1 p
[
∫
x1 x0p1(x )
i p
(|p1 | + |p2 |)
dx +
∫ xx p2(x) dx] (1)
2√ |p1 ||p2 |
η21(x)
1
(2)
where p1 and p2 are evaluated at x, η21 ) 〈φ2|dφ1/dx〉 is the nonadiabatic coupling, φj is the jth adiabatic electronic wave function, and 〈...〉 denotes integration over the electronic coordinates. The contribution from all trajectories with a single T-type hop is obtained by integrating Ψ(x0, x1, x) over all possible hopping points x1, that is, for x0 < x1 < x. The amplitude for an R-type hop is23,33
F21(x) ) -
(|p2 | - |p1 |) 2√ |p1 ||p2 |
η21(x)
(3)
In numerical problems, the trajectory is divided into small steps of length ∆x. It is convenient to express the wave function on surface Wj as Aj1χj1, where Aj1 ) [p1(x0)/pj(x)]1/2 is the usual semiclassical prefactor and the initial incoming wave function is on surface W1. After the first ∆x step, χ11 can be approximated by the no-hop term, χ11(x0 + ∆x) ≈ exp[(i/p)∫p1 dx′]χ11(x0), where the integral is evaluated from x0 to x0 + ∆x, and χ12 can be approximated using eq 1 as
[ pi ∫
χ21(x0 + ∆x) ≈ γ12exp
x+∆x/2 p1 x0
i p
( ) ( )(
x1j(xn) a11a12 x1j(xn-1) ) a21a22 x2j(xn-1) x2j(xn)
where pj is the momentum on Wj. One-dimensional expressions are used here. The many dimensional generalization for problems with any number of adiabatic states is provided elsewhere.33 The amplitude for the T-type hop is given by eq 2,23,33
τ21(x) ) -
all trajectories is the same as that appearing in the exponential in eq 4. Summing all of these terms results in the replacement of γ12 in eq 4 with sin(γ12).23 A similar summation of all terms with an even number of T-type hops (that is, all terms ending on W1) gives χ11(x0 + ∆x) ) cos(γ12) exp[(i/p)∫p1 dx′] χ11(x0), where the p1 integral is taken over the ∆x interval.23 Generalizing this analysis to any ∆x interval and initial state, it is found that
dx +
∫ xx+∆+∆x/2p2 dx]χ11(x0) (4) x
0
where γ12 is the integral of τ12 from x0 to x0 + ∆x. As long as ∆x is small, eq 4 provides a good numerical approximation to the integral over all single T-type hop trajectories that move from x0 to x0 + ∆x. The nonadiabatic coupling can become very large near the avoided crossing points of the adiabatic surfaces. Near these points it is useful to approximately sum the contributions from all trajectories with an odd number of T-type hops between x0 and x0 + ∆x.23,28 All these trajectories end on W2 and make contributions to χ21(x0 + ∆x). This summation can be performed by treating all hops as occurring at the midpoint of the interval, so that the phase function for
)
(5)
where xn ) x0 + n∆x, and j can be 1 or 2. The elements of the step matrix An are given by
(
( )
cos γ12eiS11/p sin γ12eiS12/p a11 a12 ) a21 a22 -sin γ12eiS21/p cos γ12eiS22/p
)
(6)
where
Sij )
∫ x+∆ x
x/2
n
pj dx +
∫ xx +∆x/2p1 dx n+1
n
(7)
and η12 ) -η21 has been used to obtain eq 6. For 1D problems of the type considered here, the wave function amplitudes at x0 + n∆x can be expressed in terms of their values at x0 as
( )
( )
x1j(xn) x1j(x0) ) AnAn-1...A1 x2j(xn) x2j(x0)
(8)
Equation 8 sums all trajectories from x0 and xn with energy E and any number of T-type hops. This matrix multiplication method allows for all contributions to the χkj(x) to be summed. In higher dimensional problems, this summation over all hopping trajectories would generally be performed using a Monte Carlo procedure.28-30 Previous calculations indicate that only T-type hops in the classically allowed region need to be considered if E > Ec to get very good transition probabilities, where Ec is the energy at the crossing point in the diabatic energy surfaces.49,50 Defining B(n)(x0, ∆x) ≡ AnA--1...A1, then Bj1(n)(x0, ∆x) is the amplitude for the component of the wave function for adiabatic state j at x0 + n∆x if the initial wave function has an amplitude of unity for state one and an amplitude of zero for state two at x0. Note that Bij(n)(xn, -∆x) ) Bji(n)(x0, ∆x), given that γ12 changes sign when the sign of ∆x is reversed, that γ21 ) -γ12, and that Sij for the trajectory from xm to xm + ∆x equals Sji for the reversed trajectory from xm + ∆x to xm. Now consider the case in which both adiabatic surfaces, W1(x) and W2(x), approach constant values at large x and increase rapidly for small x. At a given energy E, there is a turning point, xtj, for classical motion on each surface Wj. If W1 is the upper adiabatic surface, then xt1 > xt2. Let ∆x ) (x0 - xt1)/n and choose the incoming wave function to be on surface one, Ψ1(in)(x) ) exp[-ip10(x - x0)/p], where x0 is in the large x asymptotic region and p10 ) p1(x0). In this case, Bj1(n)(x0, -∆x) ≡ bj1(C-) is the incoming wave function amplitude on surface Wj at xt1. Throughout this work, the region x < xt2 is referred to as region A, xt2 < x < xt1 is referred to as region B, and x > xt1 is referred to as region C. The superscript C- on bj1(C-) denotes that this amplitude accounts for propagation in the negative direction across region C to xt1. If j )1, this contribution to the wave
15968 J. Phys. Chem. B, Vol. 112, No. 50, 2008
Herman
function picks up the usual e-iπ/2 turning point phase factor51,52 at xt1 and a factor of bk1(C+) ≡ Bk1(n)(xn, ∆x) for the propagation back from xt1 to x0 ending in state k. If j ) 2, the trajectory on W2 at xt1 continues in the negative direction for x < xt1. If hopping in the region between xt2 and xt1 is ignored, as in previous calculations, then the component of the wave function on W2 picks up a semiclassical phase factor of eiR/p as the trajectory travels from xt1 to xt2, where R ) ∫p2 dx′, with p2 < 0, and the integration is from xt1 to xt2. This yields an incoming amplitude at xt2 of eiR/p b21(C-). This contribution to the wave function picks up a factor of e-iπ/2 at the turning point, xt2, another factor of eiR/p for the propagation from xt2 back to xt1, and a factor of bk2(C+) for the propagation from xt1 to x0. Summing these contributions yields bk1(C+)e-iπ/2b11(C-) + bk2(C+)e2iR/ pe-iπ/2b21(C-) as the outgoing wave function amplitude on surface k at x0. This expression contains contributions from all classically allowed hopping and nonhopping trajectories, including only T-type hops. It is the expression used in previous computations.23,28,33 It is typical for the nonadiabatic coupling to be peaked around a crossing point, xc, for the diabatic surfaces.49,50 At high energies, xc is sufficiently larger than xt1 for the type of problem considered here, and the coupling is largely in the classically allowed region. If this is the case, including only hops that occur in the classically allowed region, as described above, is a very good approximation. As the energy is lowered, xt1 increases. At some energy, xt1 ) xc. Below this energy, xt1 > xc, and the coupling is peaked in the forbidden region. If this is the case, it is necessary to include hops for x < xt1 in order to obtain accurate results. If xt2 < x < xt1, then p1 ) (2m[E - W1(x)])1/2 is imaginary, and p2 is still real. The surface-hopping wave function expansion is still valid in this region, and it is useful to think of the trajectory as continuing into this region with x real, time imaginary, and momentum q1 ) -ip1 when traveling in the negative direction and q1 ) ip1 when traveling in the positive direction. The expression for the transition amplitude τ21, eq 2, becomes
τ21(x) ) -
(|p2 | + i|p1 |) 2√ |p1 ||p2 |
η21(x)
(9)
and τ12(x) ) iτ21(x). The factor of i in τ12 and it absence in τ21 can be understood as follows. An analysis of the surface-hopping wave function in the Schrodinger equation23,33 shows that τ21(x) ) [-(p1 + p2)/p2][A1/A2]η21 for a W1 to W2 hop, where Aj ) [ pi(x0)/pj(x) ]1/2 is the semiclassical prefactor at the hopping point x for a trajectory starting on initial surface i. The ratio A1/A2 gives a factor of [p2/p1]1/2, yielding eq 2 if p1 is real and eq 9 if p1 is imaginary. If the hop is from surface 2 to surface 1, then τ12(x) ) [-(p1 + p2)/p1][A2/A1]η12. Using (p1 + p2)/p1 ) (p2 + ip1)/ip1 for imaginary p1, together with η12 ) -η21, it is readily shown that τ12 ) iτ21. In the region between the two turning points, the amplitude for R-type hops is given by F21 ) -τ/21 and F12 ) -iF21, as can be seen by substituting for the imaginary p1 and comparing eqs 2, 3, and 9. Thus, the magnitude of the R-type transition amplitude, F21, is equal to the magnitude for the T-type transition amplitude, τ21, in this region, in contrast to the allowed region where F21 is always less than τ21. Consequently, the R-type hops cannot be ignored if hops in region B are included in the calculations. In this region, the T-type transitions give rise to a step matrix An, which is calculated using eqs 6 and 7 with p1 ) -ip1 and p2 ) -p2. Since τ12 is complex and τ21 ) -iτ12 in
this region, the summation of all even (odd) hop terms does not approximate a real valued cosine (sine) series. For this reason, the sin(γ12) factors in eq 6 are replaced with γ12 and the cos(γ12) factors are dropped in the calculations for this region. These lower order expressions should still result in a good approximation for small ∆x. If only T-type hops were included, B(n)(xt1, -∆x) ≡ AnAn-1...A1 ≡ b(B-), where ∆x ) (xt1 - xt2)/n, accounts for all hopping trajectories traveling in the negative direction between xt1 and xt2. Similarly, b(B+) ) B(n)(xt2, ∆x) accounts for all (T-type only) hopping trajectories traveling in the positive direction between xt2 and xt1. Since τ21 ) -iτ12, and accounting for the change in the sign of ∆x, b21(B+) ) ib12(B-), b12(B+) ) -ib21(B-), b11(B+) ) b11(B-), and b22(B+) ) b22(B-). To simplify the calculation with R-type hops, trajectories with at most one R-type hop between xt1 and xt2 are included in the calculations. The total contribution from all trajectories with a single R-type hop in region B is a sum over all ∆x intervals between xt1 and xt2 of the contribution with an R-type hop in that interval n
RijB- )
∑ [ti1(k-1,+)r(k)12 t2j(k-1,-) + ti2(k-1,+)r(k)21 t1j(k-1,-)]
(10)
k)1
where tij(k-1,+) ) Bij(k-1)(xk-1,∆x), xk-1 ) xt1 - (k - 1)∆x, tij(k-1,-) ) Bij(k-1)(xt1, -∆x),
[ pi ∫
r(k) 12 ) γFexp
xk-1 xk+∆x/2(i|p1 |
]
+ |p2 |) dx
(11)
r21(k) ) ir12(k), and γF is integral of F12 over the ∆x interval. Equation 10 includes trajectories with one R-type hop and up to 2(n - 1) T-type hops in region B. If x < xt2, then a trajectory is classically forbidden on both adiabatic surfaces, both p1 and p2 are imaginary, and the amplitudes for T-type and R-type transitions are given by eqs 2 and 3. The contribution corresponding to a trajectory decays as it moves in the negative direction. A R-type hop reflects the trajectory back, and the contribution again decays as the trajectory moves in the forward direction. The resulting amplitude, Rij(A-), has the form of eq 10 with
[ p1 ∫
r(k) 12 ) γFexp
xk-1 xk+∆x/2(|p1 |
]
+ |p2 |) dx
(12)
and r21(k) ) r12(k). Rij(A-) includes contributions from all trajectories that travel in the negative direction from xt2 on surface j to the point xk-1 ) xt2 - (k - 1)∆x with k > 0 and with up to k - 1 T-type hops, have a R-type hop in the kth ∆x inteval, and then travel in the positive direction from xk-1 to xt2 with up to k - 1 T-type hops. In this region, tij(k-1,-) ) tji(k-1,+), as is the case in the allowed region. The standard analysis for matching semiclassical expressions across turning points give the following results.51,52 An incoming trajectory in the allowed region results in outgoing trajectories in both the allowed and forbidden regions when it encounters the turning point, xt, where incoming (outgoing) means traveling toward (away from) xt when discussing turning points. The amplitude associated with the outgoing trajectory in the forbidden region is equal to the amplitude associated with the incoming trajectory multiplied by a factor of e-iπ/4 at the turning point, and the amplitude associated with the outgoing trajectory
Hops in the Forbidden Region
J. Phys. Chem. B, Vol. 112, No. 50, 2008 15969
in the allowed region is equal to the incoming trajectory’s amplitude multiplied by a factor of e-iπ/2. An incoming trajectory in the forbidden region results in an outgoing trajectory in the allowed region; the amplitude of which is equal to the amplitude for this incoming trajectory at the turning point multiplied by eiπ/4. This incoming trajectory from the forbidden region also gives rise to an outgoing trajectory in forbidden region, but this trajectory, which is reflected back into the forbidden region, is ignored in this work to simplify the calculations. These turning point phase factors (i.e., e-iπ/4, e-iπ/2, and eiπ/4) are used in the multistate calculations presented below. Let χkj(+)(x) and χkj(-)(x) be the state j to state k wave function amplitudes corresponding to propagation in the forward and backward directions at x, respectively. At the turning point, xt1, the matching conditions give χ1j(+)(xt1+) ) χ1j(-)(xt1+)e-iπ/2 + χ1j(+)(xt1-)eiπ/4 and χ1j(-)(xt1-) ) χ1j(-)(xt1+)e-iπ/4, where xt1- (xt1+) is a point x infinitesimally smaller (larger) than xt1. The quantities χ1j(-)(xt1+) and χ1j(+)(xt1-) represent incoming wave functions amplitudes at xt1, and χ1j(-)(xt1-) and χ1j(+)(xt1+) are the outgoing wave function amplitudes. The contribution to the outgoing wave in the forbidden region, χ1j(-)(xt1-), due to the incoming wave from the forbidden region, χ1j(+)(xt1-), is ignored in these expressions. The wave function amplitudes on surface two, χ2j(+) and χ2j(-) are continuous at xt1. The same analysis is applied at the surface two turning point, xt2, with the subscripts 1 and 2 interchanged. The amplitudes for propagation in the negative direction at xt1- (i.e., traveling into region B) are given by χ1j(-)(xt1-) ) e-iπ/4[b11(C-) χ1j(-)(x0) + b12(C-) χ2j(-)(x0)] and χ2j(-)(xt1-) ) b21(C-) χ1j(-)(x0) + b22(C-) χ2j(-)(x0). These give rise to amplitudes for propagation in the negative direction at xt2- (i.e., traveling into region A) of χ1j(-)(xt2-) ) b11(B-) χ1j(-)(xt1-) + b12(B-) χ2j(-)(xt1-) and χ2j(-)(xt2-) ) e-iπ/4[b21(B-) χ1j(-)(xt1-) + b22(B-) χ2j(-)(xt1-)], where the e-iπ/4 is due to the turning point xt2. These lead to forward propagating amplitudes at xt2+ of χ1j(+)(xt2+) ) R11(A-) χ1j(-)(xt2-) + R12(A-) χ2j(-)(xt2-) and χ2j(+)(xt2+) ) e-iπ/2 χ2j(-)(xt2+) + eiπ/4[R21(A-) χ1j(-)(xt2-) + R22(A-) χ2j(-)(xt2-)] at xt2+. The forward propagating amplitudes at xt1+ are then given by χ1j(+)(xt1+) ) e-iπ/2χ1j(-)(xt1+) + eiπ/4[b11(B+) χ1j(+)(xt2+) + b12(B+) χ2j(+)(xt2+) + R11(B-)χ1j(-)(xt1-) + R12(B-) χ2j(-)(xt1-)] and χ2j(+)(xt1+) ) b21(B+) χ1j(+)(xt2+) + b22(B+) χ2j(+)(xt2+) + R21(B-) χ1j(-)(xt1-) + R22(B-) χ2j(-)(xt1-). Finally, the outgoing wave function amplitudes at x0 can be expressed as χ1j(+)(x0) ) b11(C+) χ1j(+)(xt1+) + b12(C+) χ2j(+)(xt1+) and χ2j(+)(x0) ) b21(C+) χ1j(+)(xt1+) + b22(C+) χ2j(+)(xt1+). These expressions can be combined to obtain the outgoing amplitudes at x0, χ1j(+)(x0) and χ2j(+)(x0), in terms of the incoming amplitudes at x0, χ1j(-)(x0) and χ2j(-)(x0). The semiclassical transition probability for a transition from state one to state two is given by the ratio of the outgoing flux to the incoming flux, (-) 2 2 PS)p2(x0)|A21(x0)χ(+) 21 (x0)| /[[p1(x0)|A11(x0)χ11 (x0)| ] ) (-) 2 2 |χ(+) 21 (x0)| /|χ11 (x0)| (13)
It should be noted that the approximation presented here allows for, at most, one R-type hop in each of the regions A and B, although it allows for multiple T-type hops. This significantly simplifies the calculations, compared with calculations that include any number of momentum reversals accompanying hops. Calculations that allow any number of R-type hops have been performed, but this did not qualitatively change the accuracy of the method.
Figure 1. The diabatic potential surfaces V11, V22, and V12 are plotted. V12 is multiplied by a factor of 10.
III. Results In this section, semiclassical transition probabilities are compared to quantum results for the 1D, two state model defined by the diabatic49,50 energy surfaces V11(x) ) A11 exp(-a11x) + V0, V22(x) ) A22 exp(-a22x), and V12(x) ) A12{1 - tanh[a12(x - x12)]} with A11 ) 1.0, A22 ) 5.0, A12 ) 0.01, a11 ) 1.0, a22 ) 2.0, a12 ) 2.0, x12 ) 2.0, and V0 ) 0.1. These are plotted in Figure 1. The diabatic surfaces V11 and V22 cross at xc ) 1.2975, and their value at the crossing point is Ec ) 0.3732. The adiabatic state energies and wave functions are obtained by diagonalizing the diabatic energy matrix, giving W1(x) ) Vav(x) + R(x)/2, and W2(x) ) Vav(x) - R(x)/2, where Vav(x) ) (V11 + V22)/2 and R(x) ) [(V11 - V22)2 + 4V122]12. The nonadiabatic coupling is given by η21(x) ) 〈η2 dφ1/dx〉 ) [(V11 - V22) dV12/ dx - V12 (dV11/dx - dV22/dx)]/R2. Since W1 is the upper adiabatic state, the classical turning point for this surface, xt1, is larger than xt2, the classical turning point for W2, at all energies. The semiclassical surface-hopping calculations are performed using the matrix multiplication method as outlined in the previous section. The R-type terms are neglected in the classically allowed region, region C. This has consistently been shown23,28,33 to provide accurate results in previous calculations for energies at which the crossing point of the diabatic surfaces is in the allowed region. Both T-type and R-type hops are included in the calculations in region B and in region A. If R-type hops were neglected in region A, the trajectories would continue to travel away from the allowed region with a decaying amplitude. The inclusion of the R-type hop reverses the direction of the trajectory, bringing it back toward the allowed region. Without these R-type terms the effect of the nonadiabatic coupling in region A would be completely neglected in the calculations. For energies below Ec, the nonadiabatic coupling is largest in region A. Quantum calculations are performed for comparison with the semiclassical calculations. At each energy, two independent quantum calculations are started at a point xmin in region A, and the two surface Schrodinger equation is integrated in the diabatic representation using a fourth-order Runge-Kutta routine until x is equal to xf ) 10, which is well into the asymptotic regime. One of the two calculations begins with the wave function solely on V11 at xmin, and the other calculation starts with the wave function on V22 at xmin. Care is taken to avoid numerical problems due to overflows and the potential near linear dependence of the two solutions of the two state problem due to the exponential growth of the wave functions in the classically forbidden region. After both solutions have been evaluated for all x e xf at a given energy, the appropriate linear combination is taken of the two independent solutions.
15970 J. Phys. Chem. B, Vol. 112, No. 50, 2008
Herman
TABLE 1: Comparison of Quantum (PQ) and Semiclassical (PS1, PS2, and PS3) Transition Probabilities for E > Eca
TABLE 2: Comparison of Quantum (PQ) and Semiclassical (PS2 and PS3) Transition Probabilities for E < Eca
E
PQ
PS1
PS2
PS3
E
PQ
PS2
PS3
PS2(FO)
0.38 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 1.00 1.20 1.40 1.60 1.80 2.00
0.618 0.951 0.142 0.835 1.02 × 10-2 0.543 0.599 1.42 × 10-2 0.356 0.637 0.132 0.118 0.327 1.87 × 10-2 0.184 0.485 3.44 × 10-2 0.359
0.440 0.819 0.179 0.761 7.76 × 10-3 0.508 0.589 1.77 × 10-2 0.348 0.650 0.137 0.120 0.327 1.77 × 10-2 0.182 0.487 3.41 × 10-2 0.359
0.576 0.918 0.143 0.838 1.00 × 10-2 0.544 0.599 1.41 × 10-2 0.356 0.638 0.132 0.118 0.327 1.86 × 10-2 0.184 0.485 3.43 × 10-2 0.359
0.587 0.918 0.143 0.838 1.00 × 10-2 0.544 0.599 1.41 × 10-2 0.356 0.638 0.132 0.118 0.327 1.86 × 10-2 0.184 0.485 3.43 × 10-2 0.359
0.36 0.34 0.32 0.30 0.28 0.26 0.24 0.22 0.20 0.19 0.18 0.17
0.275 8.65 × 10-2 1.93 × 10-2 3.00 × 10-3 3.16 × 10-4 2.14 × 10-5 8.54 × 10-7 1.77 × 10-8 1.49 × 10-10 8.50 × 10-12 3.06 × 10-13 5.33 × 10-15
0.261 8.54 × 10-2 1.94 × 10-2 3.03 × 10-3 3.19 × 10-4 2.16 × 10-5 8.55 × 10-7 1.76 × 10-8 1.74 × 10-10 7.07 × 10-12 1.3 × 10-13 2.4 × 10-15
0.290 0.105 2.51 × 10-2 4.00 × 10-3 4.26 × 10-4 2.92 × 10-5 1.19 × 10-6 2.60 × 10-8 2.81 × 10-10 1.39 × 10-11 4.1 × 10-13 4.8 × 10-14
0.288 8.89 × 10-2 1.97 × 10-2 3.05 × 10-3 3.20 × 10-4 2.15 × 10-5 8.56 × 10-7 1.81 × 10-8 1.42 × 10-10 7.10 × 10-12 3.4 × 10-13 8.5 × 10-15
a
PS1 only includes T-type hops in region C. PS2 includes T-type hops in region C, and T and R-type hops in region B. PS3 includes T-type hops in region C and T- and R-type hops in regions A and B.
This linear combination corresponds to no incoming component of the wave function on V22 and an incoming component of the wave function on V11 of Ψ1(in)(x) ) exp[-ip1(x - xf)]. The transition probability is evaluated as the ratio of the outgoing flux on V22 and the incoming flux, PQ ) p2Ψ22/p1, where pj is the magnitude of the momentum on surface Vjj corresponding to the energy E, and all quantities are evaluated at xf. Transition probabilities are presented in Table 1 for energies above the energy of the crossing point for the diabatic surfaces. Three levels of approximation are shown for the semiclassical results. PS1 is the semiclassical transition probability when hops in regions A and B are neglected. This is the level of approximation employed in previous calculations. PS2 includes hops in region B, but not in region A (i.e., x < xt2), whereas PS3 includes hops in all regions. Hops in the forbidden regions contribute very little at higher energies, as expected, but they do improve the accuracy of the calculations at lower energies as the turning point xt1 approaches the crossing point. Transition probabilities for energies below the crossing point energy are shown in Table 2. The PS2 transition probabilities, which include hops in regions B and C, but not in region A, are somewhat more accurate than the PS3 results. These results are consistently within 20% of the quantum results for transition probabilities as small as 10-11 and within a factor of a little more than 2 for transition probabilities down to less than 10-15. The results in the last three rows of Table 2 are numerically accurate to only about ( 2 in the last place shown, due to the large cancelation between the contributions from regions B and C. The PS2(FO) results in Table 2 are similar to the PS2, except that they include at most one hop in region B; that is, all contributions to b21(B() and R21(B-) have one hop and contributions to b11(B() and b22(B() have no hops. Any number of hops are still included in the calculation of the b(C+) and b(C+) matrices, since it is known that multiple hops are needed in region C to get accurate result when E is near or above the crossing point energy. In these calculations, all the terms included in the outgoing wave function amplitude χ21(+)(x0) contain one and only one of the following: b21(B-), b21(B+), b21(C-), b21(C+), or R21(B-), whereas all other factors in each term are diagonals element (b11 or b22) from the matrices
a PS2 and PS3 are defined as in Table 1. The numbers in the last three rows are converged only to within (2 in the last place given. PS2(FO) includes the same hops as PS2, but only includes at most one hop for x e xt1.
b(B-), b(B+), b(C-), or b(C+). If the b(C-) or b(C+) were restricted to no-hop and single-hop terms, then this would be a firstorder (FO) PS2 calculation. Restricting the calculation to single-hop trajectories in the forbidden region would simplify the calculation in higher dimensions. The accuracy of the PS2(FO) results is very similar to the full PS2 calculation when E < Ec. At the lowest two energies considered, it provides results that are somewhat more accurate than the full PS2 calculation, although this improvement may be fortuitous given that it is only seen at the two lowest energies. IV. Discussion The results presented in the previous section demonstrate that the semiclassical surface-hopping expansion employed is able to provide accurate transition probabilities for the case where the crossing point and, therefore, the maximum in the coupling are well within the forbidden region. These calculations include both T-type and R-type hops in the forbidden region. The results indicate that the inclusion of hops in the region where the lower adiabatic surface is classically allowed and the upper surface is forbidden is crucial in obtaining good transition probabilities. These contributions also lead to noticeably improved results at energies slightly above the crossing point energy. On the other hand, ignoring hops in the region where both surfaces are classically forbidden generally yields slightly better results at energies below the crossing point energy than when these hops are included. The amplitudes for the T-type and R-type hops, eqs 2 and 3, are singular at turning points in the classical motion, although these singularities are integrable. It seems surprising that a primitive semiclassical method with turning point singularities should so accurately reproduce the quantum results. The small transition probability at low energies results from the considerable cancelation between terms with hops in the allowed regions and those that have hops for xt2 < x < xt1. This is most easily seen in the PS2(FO) calculations in Table 2. There are four terms in these calculations, each of which have one b21 or R21(B-) amplitude. These four terms have the form b21(C+)e-iπ/2b11(C-), b22(C+)b22(B+)e-iπ/2b22(B-)b21(C-),b22(C+)b22(B+)e-iπ/2b21(B-)e-iπ/4b11(C-), and b22(C+)R21(B-)e-iπ/4b11(C-). The first and second of these have their hop in region C, whereas the third and fourth have their hop in region B. The values for these four terms in the PS2(FO) calculation at E ) 0.2 are given in Table 3 along with their sum. The value of the contribution arising from R-type hops
Hops in the Forbidden Region
J. Phys. Chem. B, Vol. 112, No. 50, 2008 15971
TABLE 3: Contributions to Semiclassical Transition Probabilitiesa term (C+) -iπ/2
(C-)
b21 e b11 b22(C+)b22(B+)e-iπ/2b22(B-)b21(C-) b22(C+)b22(B+)e-iπ/2b21(B-)e-iπ/4b11(C-) b22(C+)R21(B-)e-iπ/4b11(C-) sum of terms 1-4 b22(C+)b22(B+)R21(A-)b11(B-)e-iπ/4b11(C-) R21(C-) a
E ) 0.2
E ) 1.0
6.50398 × + 3.05368 × 10 i 6.18492 × 10-3 - 3.65697 × 10-3 i -6.17858 × 10-3 + 3.66455 × 10-3 i -6.49839 × 10-3 - 3.06182 × 10-3 i
-3.45899 × 10 + 3.85245 × 10-1 i 3.80373 × 10-1 + 7.01941 × 10-2 i 6.94549 × 10-4 - 5.31652 × 10-4 i -6.98721 × 10-4 + 5.26158 × 10-4 i
1.19235 × 10-5 - 5.66904 × 10-7 i
3.45779 × 10-1 + 4.55434 × 10-1 i
10-3
-3
2.82697 × 10-6 - 1.24408 × 10-7 i -3.34119 × 10-3 + 3.73145 × 10-3 i
-2
1.87570 × 10-23 + 2.47054 × 10-23 i -1.10420 × 10-4 + 2.90464 × 10-3 i
See text for details.
for x < xt2, which is included in the PS3 calculation, is also given in the table. The summation of the four terms in PS2(FO) results in almost 3 orders of magnitude in cancelation. Because the interaction is relatively weak in the allowed region at this energy, there is only a small difference between the elements of b(C-), or b(C+), and the corresponding elements that would be obtained in a true first-order (i.e., no multiple-hop trajectories) approximation. If the expressions for the second and third terms in Table 3 are summed in the first-order limit, this can be reexpressed as an integral with the integrand of the form of eq 1 and limits of integration of xt2 and x0. Likewise, terms one and four can also be combined into a single integral from xt2 to x0 in the first-order limit. If the integration path is deformed into the lower half of the complex plane, the value of the integral is not altered, if the deformation avoids the square root branch points in the complex plane where the adiabatic surfaces W1 and W2 are equal. Such a deformation of the integration path is shown in Figure 2. The deformed path continues from xt2 a small distance along the real axis to a minimum value xm < xt2 and then follows a contour below the real axis back to x0 on the real axis. This deformed path avoids the singularities at xt1 and xt2 and the large cancelation that occurs when the integration is taken along the real axis. If the integral corresponding to the first-order limit of terms two and three is calculated along this deformed contour, the squared magnitude of the result is 0.95 × 10-10. The value of the
integral is numerically found to be independent of the choice of xm, as it should be, as long as the contour does not go around the branch point for the adiabatic surfaces, xb. The squared magnitude of the sum of terms two and three in Table 3 is 0.98 × 10-10. The slight difference between this result and the one from the integration along the deformed contour is due to the fact that the values in Table 3 include multiple hops in region C. R-type hops are ignored for x > xt1 in these calculations, as has been the case in previous work. The contribution from R-type hops in region C, R21(C-), is presented for E ) 0.2 and E ) 1.0 in the last line of Table 3. For E ) 1.0, the two terms involving T-type hops in region C dominate the transition probability, and the R21(C-) term is orders of magnitude smaller. However, this term increases somewhat as the energy is lowered, and it is no longer negligible at E ) 0.2. If it were included in the calculation at this energy, the transition probability would be on the order of 10-5, whereas the correct answer is on the order of 10-10. Unlike the terms with T-type hops in region C, there is no compensating term from region B to largely cancel R21(C-), which is dominated by the region near xt1 where the semiclassical approximation is inaccurate and diverging. The fact that there is no cancelation between this near-turning-point contribution and a term with hops in region B is the reason why the transition probability would be very poor if this term were included for E < Ec. The surface-hopping expansion33 is a multistate, multidimensional formalism. Because it is necessary to accurately account for the significant cancelation between contributions from allowed and forbidden regions, the accurate application of this formalism to highly forbidden transitions in multidimensional problems remains an open question at this time. On the other hand, this formalism should certainly be able to provide useful corrections and better accuracy for multidimensional problems where the avoided crossing region is near the turning point region. Conclusion
Figure 2. Integration contours in complex plane. The dark line along the real axis runs from xm to x0. (See text for details.) The semicircular sections along this line near the turning points are shown, since these small deviations from the real axis connect the correct branch of pj ) (2m[E - Wj(x)]) for x > xtj with the correct branch of pj for x < xtj. The actual integrations are performed in the limit in which the radius of the semicircle goes to zero (i.e., along the real axis). The dashed line is a contour in the lower complex half-plane that connects xm to x0 and does not enclose the branch point for the adiabatic surfaces, xb.
The semiclassical surface-hopping expansion presented previously has been applied to the calculation of transition probabilities for a 1D problem in which the avoided crossing is in the classically forbidden region. It is found that the surfacehopping calculations are capable of providing very good values of the transition probabilities down to values on the order of 10-11. This seemingly surprising result, given the singular nature of the semiclassical expression at the turning points in the classical motion, is readily understood when the integrations over hopping points are viewed as contour integrations in the complex plane. To obtain accurate results, it is necessary to include both transmission and reflection type hops in the region
15972 J. Phys. Chem. B, Vol. 112, No. 50, 2008 where the upper adiabatic surface is classically forbidden and to ignore reflection type hops in the classically allowed region. Acknowledgment. This work is funded by NSF grant CHE0715333. References and Notes (1) (a) Miller, W. H., W. H. AdV. Chem. Phys. 1974, 25, 69. (b) 1975, 30, 77. (2) (a) Marcus, R. A. J. Chem. Phys. 1971, 54, 3965. (b) 1972, 56, 311.; (c) 1973, 59, 5135. (3) Heller, E. J. Acc. Chem. Res. 1981, 14, 368. (4) Heller, E. J. J. Chem. Phys. 1991, 94, 2723. (5) Herman, M. F.; Kluk, E. Chem. Phys. 1984, 91, 27. (6) Kluk, E.; Herman, M. F.; Davis, H. L. J. Chem. Phys. 1986, 84, 326. (7) Kay, K. G. Annu. ReV. Phys. Chem. 2005, 56, 255. (8) Zor, D.; Kay, K. G. Phys. ReV. Lett. 1996, 76, 1990. (9) Madhusoodanan, M.; Kay, K. G. J. Chem. Phys. 1998, 109, 2644. (10) Sun, X.; Miller, W. H. J. Chem. Phys. 1997, 106, 6346. (11) Sun, X.; Wang, H. B.; Miller, W. H. J. Chem. Phys. 1998, 109, 7064. (12) Wang, H.; Thoss, M.; L Sorge, K.; Gelabert, R.; Gimenez, X.; Miller, W. H. J. Chem. Phys. 2001, 114, 2562. (13) Liu, J.; Miller, W. H. J. Chem. Phys. 2006, 125, 224104. (14) Shao, J.; Makri, N. J. Phys. Chem. A 1999, 103, 7753. (b) 1999, 103, 9479. (15) (a) Delos, J. B.; Thorson, W. R. Phys. ReV. A 1972, 6, 720. (b) 1972, 6, 728. (16) Miller, W. H.; George, T. F. J. Chem. Phys. 1972, 56, 5637. (17) (a) Meyer, H. D.; Miller, W. H. J. Chem. Phys. 1979, 70, 3214. (b) 1979, 71, 2156; 1980, 72, 2272. (18) Stock, G.; Thoss, M. Phys. ReV. Lett. 1997, 78, 578. (19) Stock, G.; Thoss, M. Phys ReV. A 1999, 59, 64. (20) (a) Bonella, S.; Coker, D. F. J. Chem. Phys. 2001, 114, 7778. (b) 2003, 118, 4370. (21) Preston, R. K.; Tully, J. C. J. Chem. Phys. 1971, 54, 4297. (22) Tully, J. C.; Preston, R. K. J. Chem. Phys. 1971, 55, 562. (23) Herman, M. F. J. Chem. Phys. 1982, 76, 2949. (24) (a) Herman, M. F. J. Chem. Phys. 1984, 81, 754. (b) 1984, 81, 764.
Herman (25) (a) Herman, M. F. J. Chem. Phys. 1985, 82, 3666. (b) 1995, 103, 8081. (26) Herman, M. F.; El Akramine, O.; Moody, M. P. J. Chem. Phys. 2004, 120, 7383. (27) Moody, M. P.; Ding, F.; Herman, M. F. J. Chem. Phys. 2003, 119, 11048. (28) Herman, M. F.; Moody, M. P. J. Chem. Phys. 2005, 122, 094104. (29) Wu, Y.; Herman, M. F. J. Chem. Phys. 2005, 123, 144106. (30) Yang, G.; Herman, M. F. J. Phys. Chem. B 2001, 105, 6562. (31) Wu, Y.; Herman, M. F. J. Chem. Phys. 2006, 125, 154116. (32) Wu, Y.; Herman, M. F. J. Chem. Phys. 2007, 127, 044109. (33) Herman, M. F.; Wu, Y. J. Chem. Phys. 2008, 128, 114105. (34) Tully, J. C. J. Chem. Phys. 1990, 93, 1061. (35) Hammes-Shiffer, S.; Tully, J. C. J. Chem. Phys. 1994, 101, 4657. (36) Burant, J. C.; Tully, J. C. J. Chem. Phys. 2000, 112, 6097. (37) Fang, F.-Y.; Hammes-Shiffer, S. J. Chem. Phys. 1997, 107, 8933. 1999, 110, 11166. (38) Webster, F.; Wang, E. T.; Rossky, P. J.; Friesner, R. A. J. Chem. Phys. 1994, 100, 4847. (39) Bittner, E. R.; Rossky, P. J. J. Chem. Phys. 1995, 103, 8130. 1997, 107, 8611. (40) Ben-Nun, M.; Martinez, T. J. J. Chem. Phys. 1998, 108, 7244. (41) L Volovuev, Y.; Hack, M. D.; Topaler, M. S.; Truhlar, D. G. J. Chem. Phys. 2000, 112, 9716. (42) Hack, M. D.; Wensmann, A. M.; Truhlar, D. G.; Ben-Nun, M.; Martinez, T. J. J. Chem. Phys. 2001, 115, 1172. (43) Jasper, A. W.; Hack, M. D.; Truhlar, D. G. J. Chem. Phys. 2001, 115, 1804. (44) Jasper, A. W.; Stechmann, S. N.; Truhlar, D. G. J. Chem. Phys. 2002, 116, 5424. (45) Jasper, A. W.; Truhlar, D. G. Chem. Phys. Lett. 2003, 369, 60. (46) Zhu, C.; Kamisaka, H.; Nakamura, H. J. Chem. Phys. 2002, 116, 2324. (47) Kondorshiy, A.; Nakamura, H. J. Chem. Phys. 2004, 120, 8937. (48) Bremmer, H. Commun. Pure Appl. Math. 1951, 4, 105. (49) Delos, J. B. ReV. Mod. Phys. 1981, 53, 287. (50) Smith, F. T. Phys. ReV. 1969, 179, 111. (51) Maslov V. P.; Fedoriuk, M. V. Semi-Classical Approximation in Quantum Mechanics; D. Reidel Publishing: Dordrecht, 1981. (52) Schiff, L. I. Quantum Mechanics; McGraw-Hill: New York, 1968; pp 272-275.
JP804937Q
