Intramolecular quantum dynamics on intersecting potential energy surfaces: a tutorial review

Abstract

Vibronic coupling-i.e. the interaction between molecular electronic states through the nuclear motion – is an ubiquitous phenomenon for molecular electronic spectra and dynamics in excited electronic states. While the Born–Oppenheimer or adiabatic separation of electronic and nuclear motions remains the cornerstone of our thinking in the field, the ensuing limitations are often thought of almost in line with the theorem itself. In this short review paper, we attempt to give a lucid overview over some basics in the field complemented by relevant examples. While written from a personal perspective, it is nevertheless hoped to provide some useful insight for the beginner and expert alike.

---

I. Brief historic introduction

The well-known Born–Oppenheimer approximation (sometimes also termed adiabatic approximation) is fundamental to our thinking of nuclear motion in molecules. Through the large difference between nuclear and electronic masses, the latter are thought to follow the former almost instantaneously and thus provide a potential for the (slowly changing) nuclear positions. The large difference in masses is reflected in normally much larger electronic than nuclear energy spacings.

Nevertheless, already in 1969 and earlier the importance of small electronic energy gaps was pointed out and argued to lead to significant nonradiative transitions between different molecular electronic states.^1,2 Rartly in parallel the theory of the Jahn–Teller effect has been developed, where the consequences of symmetry-induced degeneracies and their implications for structural distortions and excitation of non-totally symmetric vibrational modes have been elucidated.^3,4 The electronic degeneracy at the high-symmetry nuclear configuration, and its lifting through asymmetric distortions, make up for an early example of a conical intersection of electronic potential energy (PE) surfaces where the nuclei may interconvert freely between the different PE surfaces. Nowadays this has emerged as the typical scenario for rapid interconversion between different PE surfaces, and conical intersections are now established as the paradigm for nonadiabatic excited-state dynamics in the literature.^5–11 Quite generally, such intersections are not expected to be rare events and, for real electronic wavefunctions (that is, no spin–orbit coupling), only two conditions need to be satisfied to achieve degeneracy of two different PE surfaces.^1,12 While this leads to the non-crossing rule for diatomic molecules with a single nuclear degree of freedom, for triatomics and larger molecular systems all but two nuclear degrees of freedom can be varied without lifting, in principle, the electronic degeneracy (again, for systems without spin–orbit coupling).

The efforts to identify conical intersections in specific systems and to illustrate their implications can be categorized as focussing on either electronic structure or dynamical aspects. Quite naturally the electronic structure work dealt with general (i.e. not of model-based shape) PE surfaces as determined from increasingly sophisticated correlated wavefunctions.^7–10 Quantum-dynamical studies often relied on the so-called linear vibronic coupling (LVC) model, where the matrix elements of the Hamiltonian in a diabatic electronic basis arise from the potential energy and are expanded in a Taylor series in first order (in the simplest case) of the nuclear displacement coordinates.⁵ Natural extensions consist in the inclusion of higher-order terms (see, for example,^13,14) and in application of the LVC scheme to the so-called adiabatic-to-diabatic (ADT) mixing angle,¹⁵ to be discussed below.

The construction and use of diabatic electronic wave functions is motivated by the singularity of the derivative coupling matrix elements between different adiabatic electronic wave functions at conical intersections. By some textbook manipulations (see, for example, Chapter 1 by L. S. Cederbaum in ref. 6) this matrix element is found to be

is given here for comparison and completeness. In these expressions x and Q denote the sets of electronic and nuclear coordinates, respectively. The symbol ϕ refers to the electronic wave function. The difference V_n–V_m in the denominator rationalizes the earlier statement about significant nonradiative transitions between electronic states with a small energy gap. The ratio diverges for V_n → V_m which makes adiabatic electronic wave functions tedious to deal with in numerical computations on conically intersecting PE surfaces. This circumstance motivated the use of diabatic electronic wave functions in the literature which have only very crudely been characterized above and are treated in more detail in the next section.

The sections to follow below mostly focus on various examples and phenomena addressed in our work on quantum dynamics on conically intersecting PE surfaces.

II. Diabatic basis sets and adiabatic potential energy surfaces

As stated in the Introduction, diabatic electronic basis sets have been introduced into the literature in order to avoid the singularity that appears when expressing the molecular Schrödinger equation in a basis set of adiabatic electronic wave functions (see Appendix A). For the relevant polyatomic case, let us consider the case of three nuclear degrees of freedom and two interacting electronic states, and express the diabatic electronic wave functions ϕ^d_i(i = 1, 2) as suitable linear combinations of the adiabatic counterparts ϕ_i(i = 1, 2) likeIdeally, diabatic electronic wave functions should be chosen such that the derivative coupling between them vanishes, i.e.

〈ϕ^d₁|∇|ϕ^d₂〉 = 0. (3)

Inserting the above expression eqn (2) into this relation and rearringing the various terms leads to (for real orthonormal electronic wave functions ϕ_i(i = 1, 2))

∇γ = 〈ϕ₁|∇|ϕ₂〉. (4)

From elementary calculus it is known that that the curl of ∇γ (grad γ) vanishes and so should the curl of the r.h.s. of the eqn (4). This is known as curl condition in the literature, but not fulfilled in the general polyatomic case.¹⁵ Consequently, strictly diabatic states do not exist in the general case. Extensions to more than two states have also been given^15,16 and lead to the same result.

Subsequently, efforts were directed to solve the pertinent equations approximately and thus define quasidiabatic electronic states. These are nowaday in the focus of theoretical efforts (and the attribute “quasi” is often suppressed for notational simplicity). The significance of the approximation can be highlighted by being aware that the singular derivative couplings, i.e. those that diverge at a conical intersection, indeed fulfil the curl condition and can be rigorouly removed by constructing suitable quasidiabatic states.^17,18 Two examples will be given below where the focus on the singular derivative couplings indeed gives very useful results.

A big practical advantage of using quasidiabatic electronic states is that the smoothness of the electronic wavefunctions allows for a variety of well-known approximations such as power series expansions of the resulting Hamiltonian matrix elements. Normally, curvilinear or (mostly) rectilinear displacement coordinates relative to the minimum energy structure of the electronic ground state are employed. The latter is usually the initial electronic state before the spectral transition under consideration. The simplest case, using only the first-order terms, constitutes the well-known linear vibronic coupling or LVC scheme^5,6 and the Hamiltonian for the simple two-state two-mode case with nondegenerate electronic states and vibrational modes can be written as eqn (5).

Here the interacting electronic states are taken to be of different symmetry, with Q_g (g = gerade) being a symmetry-preserving coordinate and Q_u (u = ungerade) being a symmetry-lowering coordinate leading to an interaction of the two electronic states and to the off-diagonal matrix element in the Hamiltonian matrix eqn (5). κ₁, κ₂ and λ are intra-state and inter-state electron vibrational coupling constants, and T_N and V₀ denote the kinetic and potential energies of the ground state Hamiltonian in the harmonic approximation.

An analogous Hamiltonian has been employed for the Jahn–Teller (JT) effect with a doubly-degenerate electronic state and vibrational mode, the E × e JT effect.^19,20 The ensuing constraints for the coupling constants are related to Racah coefficients for the molecular point groups and can be explicitly rederived by requiring that the Hamiltonian transforms totally symmetric in the high-symmetry (JT-undistorted) molecular point group. This leads to κ = κ₁ = −κ₂ = λ, and the coordinates are to be reinterpreted as the Cartesian displacement coordinates Q_x and Q_y of the doubly degenerate JT-active vibrational mode. Then the Hamiltonian takes on the following appearance^19,20

Similar considerations and results hold for systems with nondegenerate JT-active modes (E × b JT effect) and for triply degenerate vibrational modes and electronic states (T × t JT effect).^19,20 Linear vibronic coupling (LVC) Hamiltonians like eqn (5) or (6) have been used early on to successfully model the photoelectron spectrum of butatriene²¹ (see also ref. 22) and to explore the PE surfaces and dynamics of the prototypical linear E × e Jahn–Teller effect (see, for example.^19,20 For the latter case, one often switches to complex vibrational (or also polar) coordinates

Q_± = Q_x ± iQ_y = ρe^±iϕ (7)

and similarly with electronic basis functions,The E × e JT Hamiltonian then takes the form given in eqn (B1) of Appendix B, more precisely, in the first 2 terms on the r.h.s. of eqn (B1).

In principle the LVC approximation is well suited to treat small-amplitude displacements. This limitation can be partly overcome by using curvilinear rather than rectlinear coordinates. Another possibility is to carry the Taylor-series expansion to higher order such as forth or sixth-order terms.^13,14 In the multimode case a wide variety of treating the intermode coupling terms arises, which is beyond the scope of this short article. Here, we confine ourselves to the simplest extension obtained by adding second-order terms to the Hamiltonian eqn (6). In the complex electronic and vibrational representation of eqn (7) this leads to eqn (B1) and is the general starting point of the considerations in Appendix B. Hamiltonians like the above consists of the kinetic energy T_N and a potential energy matrix W defined by

H = T_N1 + W. (9)

The eigenvalues of W represent the upper and lower adiabatic PE surfaces V_±. In the 2D linear coupling case g = 0 of eqn (B1) the eigenvalues readThe degeneracy at the origin Q_x = Q_y = ρ = 0 is lifted in first order of the nuclear displacements and the adiabatic PE surfaces thus form a conical intersection.^5,6,19,20 This constitutes one of the earliest examples of a conical intersection in the laterature, and the PE surfaces are also known as mexican hat PE surfaces.^19,20 The transformation matrix for the case of eqn (6) is a simple 2 × 2 rotation where the rotation angle or ADT transformaton angle α, for the Hamiltonian eqn (6) readsRe-expressing the quasidiabatic potential energy matrices in terms of the adiabatic PE surfaces and the ADT angle leads toThis allows for an important generalization by combining the ADT matrix of the linear coupling or LVC scheme with the higher-order, or also general, adiabatic PE surfaces V^gen_±. Augmenting the trigonometric functions in eqn (12) by a superscript “lin” to indicate the linear coupling approximation of the Hamilonian of eqn (6), this leads to the following definition of the regularized diabatic states viz PE matrix W^reg,The quasidiabatic electronic states (and PE matrix) thus constructed are termed regularized diabatic states (and PE matrix): the underlying quasidiabatic states are those pertinent to the LVC scheme while the ensuing adiabatic PE surfaces are the “full” PE surfaces without any reference to the LVC scheme.^17,18,23 The prefix “quasi” will be suppressed henceforth for notational symplicity. The crucial point is that by a correct determination of the coupling constants κ and λ,^17,18,23 this allows to treat the singular derivative couplings, i.e. those that diverge at the point of degeneracy (conical intersection), correctly.^17,18,23 In the present simple two-state–two-mode case the coupling constant κ just needs to reproduce the actual gradients of the PE surfaces at the conical point correctly.^17,23 Examples have been set up which demonstrate the success of the idea of addressing only the singular derivative couplings when constructing quasidiabatic electronic states. This is particularly simple when only two electronic states and two vibrational modes play a role like here. In ref. 23, the case of the quadratic E × e JT effect as defined in eqn (B1) has been investigated and the pertinent equations are included in Appendix B. The adiabatic PE matrix of this system was diabatized employing the eigenvectors of the linear E × e JT effect. This leads to the Hamilonian matrix, eqn (B7) within the regularized diabatic states (RDS) scheme, see also ref. 23. It preserves the adiabatic PE surfaces and also the singular derivative couplings because the latter are embodied in the linear E × e JT Hamiltonian.^19,20 The nonsingular derivative couplings are treated differently than in the original second-order Hamiltonian of eqn (B1). In Fig. 1 the absorption spectra for the full calculation and the RDS scheme are compared and are so close that the difference is easily overlooked. A similar comparison is presented in Fig. 2 for different parameter values with similarly good agreement, and the same holds for several other cases investigated by us.²³ Thus, as far as numerical experience can tell, the neglect of the nonsingular derivative couplings works very well for the type of system addressed.

Fig. 1 Comparison of exact and approximate JT spectra for κ = 1 and g = 0.7, reproduced from ref. 23 with permission from AIP, copyright 1999. The inset shows the corresponding spectrum for purely linear coupling i.e., g = 0).

Fig. 2 Comparison of exact and approximate JT spectra for and g = 0.4, reproduced from ref. 23 with permission from AIP, copyright 1999. The inset shows the corresponding spectrum for purely linear coupling i.e., (g = 0).

To further substantiate the conclusion made, we performed a more extended comparison for a three-mode system, consisting of the 1¹A₂ and 1¹B₁ excited states of SO₂.²⁴ Here two symmetric A₁ modes (“tuning” modes) and one asymmetric B₂ mode coupling the two states play a role. For equal bond lengths (C_2v point group) the two electronic states, being of different symmetry, can cross freely and the crossing seam is shown as green line in Fig. 3. Two different ways of determining λ lead to contour lines indicated by the full red and dashed blue lines in the figure. They are obviously quite different but coincide on the seam where the derivative couplings diverge, see eqn (1). The singular derivative couplings are thus identical, and some details of their determination²⁴ are collected in Appendix C. The resulting absorption spectra are compared in the lower panel of Fig. 4 and the difference is again very small indeed.

Fig. 3 Plot of the coupling constant λ in the C_2v subspace. In dashed blue λ_C2v is evaluated in the entire C_2v subspace. In dashed red λ_seam is determined along the seam of intersection and taken to be independent of ϕ. The green line represents the seam of intersection (along this seam both methods give the same value of λ) Unit of λ is hartree/Å. Figure reproduced from ref. 24 with permission from AIP, copyright 2013.

Fig. 4 Comparison of the theoretical and experimental absorption spectra of SO₂. (Top) Experimental (black line) and theoretical spectrum with λ_seam (red). (Bottom) Spectrum computed with λ_seam (red) and with λ_C2v (blue). Figure reproduced from ref. 24 with permission from AIP, copyright 2013.

As already mentioned, numerical examples can never provide complete “proof”, and we conclude this discussion with three remarks. Further support of the concept of regularized diabatic states is provided by an analysis of the residual derivative coupling matrix elements, i.e. those remaining present in the quasidiabatic basis. For the JT system of Fig. 1 and 2 they have been analyzed numerically for quite a few parameter values and indeed found to be much smaller, often by orders of magnitude smaller, than those in the adiabatic basis.²³ Thus, the smallness of electronic matrix elements remains relevant also when integrating over the nuclear coordinates.²³ Second, the complexity in the most general case seems to derive mostly from the complexity of the intersection seam which can be very complex in shape, also be self-intersecting etc. Knowledge of the location of the seam can be the vital step in dealing with the diabatization problem. Finally, a formally independent (but often important) question is the type of arguments to be used in the Taylor series, rectilinear or curvilinear coordinates or (simple) functions of them. Suffice it to say that in our own work often rectilinear (Cartesian) displacement coordinates have been used, excepting our studies of small to medium-size polyenes, where curvilinear (internal) coordinate have been employed to account for large angular (including torsional) displacements.²⁵

III. Remarks on observables

The spectra in the preceding section have been obtained, like in all our work, by quantum dynamical methods, more specifically through wave-packet propagation with the efficient Multi-Configuration Time-Dependent Hartree (MCTDH) method.^26–29 While this scheme does not overcome the exponenial scaling problem inherent to all quantum techniques, by employing the Dirac–Frenkel variational principle the expansion of the many-body wave-function in terms of time-dependent single-particle functions becomes optimally short, in relevant examples shorter by six or more orders of magnitude. This renders practical applications often manageable at all.^27,28 Given the full wave-function time-dependent observables can be directly computed as expectation values of suitable operators, such as (diabatic or adiabatic) projection operators or (auto)correlation functions. The spectral intensity distribution then follows by Fourier transformation. Further insight is obtained from reduced densities which offer a vivid picture of the time-dependent wave-function. Many useful technical details and more are contained in the original references, and for our purposes it seems sufficient to recall them here collectively again.^27–29

IV. SO₂ and other representative triatomic systems

In a series of papers we have investigated the quantum dynamics and electronic spectrocopy of a series of symmetric triatomic species, the ozone radical cation O₃⁺ (ref. 30) and singlet excited states of NO₂^31–33 and SO₂.^24,34–37 The systems are all characterized by a seam of low-energy symmetry-allowed conical intersections occurring for equal X–O bond lengths (X = O, N and S). The latter two, neutral species represented “famous” or better notorious problems in molecular spectrocopy. Suffice it here to mention some details concerning the 1¹A₂ and 1¹B₁ singlet excited states of SO₂. A main concern for a long time was the observation that the spectral transition of SO₂ in the visible range is dipole-allowed, consistent with the 1¹B₁ final electronic state, but the rotational structure points to geometric parameters characteristic for the 1¹A₂ final (dipole–forbidden) electronic state.³⁸ This has been interpreted by invoking an intricate coupling between these states. A further “pecularity” was the obervation of a Zeeman effect in the presence of a magnetic field in the low-energy part of the spectrum,^34,36 which was attributed to a triplet excited state by the analysis of Merer and Douglas.³⁹ Conflicting assignments have partially been resolved by ab initio electronic structure calculations of ref. 24 and further been clarified by our work,³⁷ to be presented in what follows. Starting point were ab initio MRCI and CCSD calculations with the cc-pVTZ basis set for the sulfur and oxygen atoms. This yields PE curves along the symmetric S–O stretching and O–S–O bending coordinates which are displayed in Fig. 5. The different bond angles in the two excited states are evident which gives rise to the low-energy curve crossing in the upper panel (or conical intersection when invoking the asymmetric stretching coordinate). The dynamical calculations rely on (symmetric) Jacobi coordinates (see Fig. 1 of ref. 24), the regularized diabatization approach for constructing the vibronic Hamiltonian and on the MCTDH wavepacket propagation scheme. The computed visible absorption spectrum of SO₂ is reproduced in the upper panel of Fig. 4 and the theoretical result is seen to match experiment very convincingly. Not only the overall features like band centre and width agree very well, but also spectral details agree almost quantitatively in the lower energy part of the spectrum. For higher energies some very minor differences exist, such as a shift nearly constant in energy which may be caused by residual inaccuracies of the PE surfaces which have here been taken without any adjustment.

Fig. 5 Cuts through diabatic 1¹A₂ (red) and 1¹B₁ (blue) PES within C_2v symmetry. (From top to bottom) Cut along the O–S–O angle ϕ with R_s kept at the equilibrium geometry of the ground state (R_s = 1.448 Å) and cut along the symmetric stretch coordinate R_s with ϕ kept at the ground state equilibrium value (ϕ = 119.1°). Arrows indicate the values at the ground state equilibrium. The dashed lines represent the ab initio potential data after diabatization and the full lines the diabatic potential curves after fitting. Figure reproduced from ref. 24 with permission from AIP, copyright 2013.

Analysis of the time-dependent quantitities reveals an intriguing feature.²⁴ The time-dependent electronic wavepacket stays moving mostly on the lower adiabatic PE surfaces. As the arrows in Fig. 5 show, the system is lifted to the lower PE surface initially by dipole selection rules. When approaching the conical intersection at later times, it bifurcates into two branches with unequal S–O bond lengths, thus avoids nuclear configurations subject to strong nonadiabatic interactions and stays confined to the lower (adiabatic) PE surface although the energy is sufficient to access the upper PE surface.²⁴ This counterintuitive behaviour is related to the specific way in which the wave-packet approaches and crosses the intersection seam and is the main reason for the complex structure of the visible absorption spectrum of SO₂.

In the course of the work also doublet and triplet electronic states of SO₂ have been addressed and, for example, the vibronic structure of the relevant photoelectron sprectral bands been related to the similarities and differences of the corresponding PE surfaces in the neutral and cationic states.³⁵ Concerning triplet states, the importance of the ³B₂ state has been predicted by extensive wavepacket simulations.³⁶ More generally, the intensity of transitions to triplet states can be much enhanced by the combined use of polarized light and molecular alignment which can even make it the only absorption process under suitable conditions.³⁴

Before proceeding to the next set of examples, we emphasize that in this short tutorial review paper the literature cannot be cited comprehensively, and even the papers cited are not all meant to be representative. The purpose is instead to give a better overview over our own work in the field. The same will apply to the sections that follow below.

V. Four-atomic and larger linear polyenes: butadiene and related polyenes

In a different line of work we have investigated butadiene and somewhat larger linear polyenes by a basically similar dynamical approach. The main difference to other non-triatomic systems treated below is the use of internal (curvilinear) rather than Cartesian (rectilinear) displacement coordinates. This renders the kinetic energy operator more complex, but can be handled by the use of the G-matrix technique of Wilson, Decius and Cross.^40,41 The systems treated are linear polyenes with mostly up to four (but also more) conjugated double bonds. Polyenes play a key role in photovoltaic devices and also serve as structural motifs in many biomolecules,^42–44 and their importance in the fields can hardly be overestimated. The UV spectra and related photodynamics of these systems are dominated by two excited singlet states, the bright 1B_u state and the dark 2A_g state. While the former is characterized by a “simple” HOMO–LUMO transition, the latter displays complicated multiconfigurational features, consisting of a HOMO−1→ LUMO, a HOMO → LUMO+1 and a HOMO² → LUMO² double excitation. The complexities due to the electronic structure are combined with those arising from nuclear motion due to the energetic proximity of the two electronic states (see, for example ref. 25 and 45–47). While for smaller systems with up to four conjugated double bonds the bright 1B_u state comes energetically below the dark 2A_g state, for systems with more than four conjugated double bonds their energetic ordering is reversed.

Some time ago we have started working on the smallest system s-trans-butadiene based on ab initio CASPT2/cc-pVTZ computations of PE-curves of the six most relevant vibratoinal modes.⁴⁸ Subsequently, 5 in-plane (3 C_2v preserving and 2 asymmetric modes) and 3 torsional degrees of freedom have been considered.⁴⁵ A prominent feature of the PE curves is a (in C_2v) crossing of the 1B_u and 2A_g state PE curves upon symmetric stretching of the terminal double bonds (mode S₃); this occurs near the minimum of the initial (bright) state and is accompanied by a substantial repulsion of the PE-curves along the asymmetric or coupling mode S₁₇ (asymmetric C–C–C bending). The torsional PE surfaces reveal interesting fourfold minma for cis and trans torsional displacements leading to coherent cis and trans excitation in the quantum dynamical treatment⁴⁸ (see also Martinez et al.⁴⁹).

With regard to the conical intersection, one finds an ultrafast population transfer from the bright 1B_u to the dark 2A_g state,⁴⁵ demonstrated vividly by Fig. 7 of ref. 45. This is typical for conically intersecting PE surfaces and has been observed also for many other cases.^5,6 When comparing the various dimensionalities of the dynamical treatment, the populations are seen to become “smoother” and the population transfer gets more efficient with an increasing number of modes.^45,48 This has been related to the different volume of the configuration space available. The spectral intensity distribution seen in Fig. 6 also displays characteristic features well known from other systems,^5,6 namely a very irregular fine structure arising from a heavy mixing of vibrational levels of upper and lower PE surfaces. It has been obtained by wavepacket propagation with a varying number of vibrational degrees of freedom as specified in the inset to the figure. Noteworthy, under low resolution, corresponding to a damping time of 20 fs, the main progression of the terminal C–C stretching mode, which dominates the fine structure of panel a shows up clearly and is also prominent in experimental recordings of related other polyenes.^50,51 For higher resolution (damping time of 1000 fs) the underlying irregularity becomes predominant. In further studies we proceeded towards the larger systems hexatriene^52,53 and octatetraene^46,54 and also methylated derivatives.⁵⁵ This documents the impact of the changing 1B_u–2A_g energy gap on the dynamics. With an increasing number of conjugated double bonds the dark state approaches the bright state from above and, starting with octatetraene, their energetic ordering becomes reversed. In the transitional regime, great care is needed because the results may depend on the details of the theoretical treatment.^52,53 A further relevant trend is the energetic de/increase of the bright state with terminal dimethylation; this amounts to ∼0.15 eV and has been attributed to hyperconjugation⁵⁶ due to the terminal methyl groups, documented in ref. 55. The change in relative 1B_u–2A_g excitation energy impacts on the minimum 1B_u–2A_g conical-intersection energy: for longer polyenic systems, while the 2A_g excitation energy decreases, the minimum 1B_u–2A_g intersection energy increases; the analogy to the Marcus inverted regime in the theory of electron transfer has been emphasized⁴⁶ and also the relevance for the fluorescence dynamics (dual fluorescence) been discussed.⁴⁶

Fig. 6 Comparison of calculated spectra of butadiene with the experimental result by Vaida et al.⁵⁰ The experiment is shown in the upper right panel. The dynamical results including 4 modes, 6 modes and 8 modes are depicted in the upper left, lower left and lower right panel, respectively. In all cases, ΔE = 0.51 eV. Figure reproduced from ref. 45. with permission from Elsevier, copyright 2017.

VI. Five-membered ring systems: furan, pyrrol, thiophene

The five-membered heteroaromatic ring systems C₄H₄X with X = O, NH and S reveal an interesting change of the first two (A₂ and B₁) ionizaton potentials and related curve crossings in the series.⁵⁷ For furan their energetic difference is substantial, the related conical intersection is rather high in energy and the spectral structures and population dynamics appear as one expects within the B–O approximation (this is nicely displayed in Fig. 1 of ref. 57). Proceeding from O through N–H to S (see Fig. 7) the A₂–B₁ curve crossing moves to lower relative energy until it reaches the B₁ state minimum for X = S. Whereas the population transfer is gradual for the first system, it proceeds on a 20 fs timescale and is almost “instantantaneus” for the other cases, see Fig. 8. This vividly demonstrates the importance of a direct access to the conical intersection (X = S, N–H) as compared to tunneling (X = O) and also becomes well visible in the vibronc structure of the photoelectron spectral bands.⁵⁷ The work has been extended to cover also singlet excited states of the neutral species,^57,58 and a more recent elaborate MCTDH study of thiophene has also become available.⁵⁹ In our work a further focus was on reactive processes such as ring opening (using Jacobi-type coordinates⁶⁰) and N–H photodissociation.⁶¹ In the former case, full coherence between the two C–O bond-breaking pathways has been maintained due to the quantum treatment performed. For the latter system, a five-state approach enabled us to analyze multimode internal conversion processes occurring parallel to the photodissociaton and to confirm a time constant (∼100 fs) for the photoreaction.

Fig. 7 Cuts through the potential energy surfaces of the furan (full lines), pyrrole (long dashes) and thiophene (short dashes) radical cations in their ²A₂ ground and ²B₁ first excited electronic states. The cuts are always taken to be straight lines interconnecting the origin Q = 0 in normal coordinate space with the energetic minimum of the ²A₂–²B₁ conical intersection. The ²A₂ minimum is always taken to be the zero of the energy scale to facilitate the comparison. Note the pronounced shift of the ²B₁ potential energy curves in the series furan-pyrrole-thiophene which leads to a marked lowering of the energy of the curve crossing relative to the minimum of the ²B₁ state in the series. Figure reproduced from ref. 57 with permission from Elsevier, copyright 2004.

Fig. 8 Time-dependent electronic populations of the ²A₂ and ²B₁ electronic states of furan, pyrrole and thiophene cations for a broad band (vertical) excitation to the ²B₁ electronic state in the centre of the Franck–Condon zone. Figure reproduced from ref. 57 with permission from Elsevier, copyright 2004.

VII. The benzene radical cation and derivatives

The benzene radical cation Bz⁺ represents the most extensive multi-mode—multi-state system treated by us to date. Up to 8 coupled PE surfaces and 28 nonseparable vibrational degrees of freedom have been included in the ab initio quantum dynamical treatment.^62,63 The X²E_1g ground state and B²E_2g first excited state (in D_6h symmetry) have been investigated also by many other authors (representative references are^64–67) and we proceed here directly by addressing a three-state study comprising the B²E_2g–C²A_2u excited state of Bz⁺ (henceforth symmetry labels will be mentioned only at first occurrence of a state and spin multiplicity will be suppressed). We mention in passing that even in the B (as well as in the X state) alone⁶⁸ the nuclear is nonadiabatic because it proceeds on the intersecting mexican hat PE surfaces mentioned already above. In the B–C manifold this is augmented by B–C intersections with their own signatures in the photoelectron spectrum and population transfer.⁶⁹ Extending the analysis to include the X state leads to the X–B–C three-state system.⁷⁰ With five coupled PE surfaces the diversity of nonadiabatic couplings is correspondingly increased and not elucidated here in detail. Most prominent appears to be the existence of a comparatively slow (compared to earlier examples) population transfer time of ∼100 fs⁷⁰ which is nevertheless responsible for the absence of fluorescence in this prominent species.^71,72 Proceeding to higher energies, one encounters the D²E_1u and E²B_2u electronic states which resemble the B–C set, but with a smaller energy gap between the degenerate and nondegenerate state.⁷³ With an otherwise similar topology of the PE surfaces also the higher photoelectron spectral bands are unstructured and diffuse.⁷³ – The full set of states becomes relevant when addressing the XUV pump – near VIS/IR probe experimental work of ref. 62. The pulses are included in Fig. 9 along with a schematic, but realistic overview over the relevant PE surface cuts along a suitable effective mode. The first HHG (pump) pulse generates radical cations in the enery range of the D and E states of Bz⁺. In the second, two-photon (probe) step the system is further excited to generate fragment ions as a function of the pump–probe time delay. Through a set of considerations and further analysis, which is beyond the scope of this short article, the fragment ion yield can be related to the electronic populations as obtained e.g. from MCTDH calculations. By a suitable weighting procedure one arrives at the composite populations displayed in Fig. 10 along with the experimental fragment yield so obtained.^62,63 Both curves feature an biexponential decay with fast time constants of 8 fs (theory) and 11 fs (experiment) and slow time constants of 170 fs (theory) and 110 fs (experiment). This represents one of the fastest internal conversion processes observed to this date in a time-resolved experiment. The theoretical calculation allows to interpret this result as reflecting the consecutive passage of the two high-energy conical intersection also depicted in Fig. 9: they involve the E/D and D/B electronic states of Bz⁺ and are energetically close or more distant, and involve smaller or rather large nuclear displacements, respectively. A more detailed analysis of the (adiabatic and diabatic) population in the multi-state system is available.⁶³ A further avenue of investigation was opened when considering benzene cation derivatives, epecifically fluoro-derivatives with one to three F-atoms.^74–77 This is expected to stabilize the sigma-orbitals of benzene but leave the pi-orbitals largely unchanged. The correponding (sigma-derived) cationic states should increase in energy, but not so much the pi-type states. This expectation is confirmed and quantified in the survey diagram of Fig. 11. This reflects on one hand the splitting of the X-state (and other degenerate states’) ionization potential (by 0.2 eV) due to the symmetry lowering from D_6h to C_2v/D_2h, and on the other hand the gradual, but continued increase of the sigma-type ionization potential, deriving from the E_2g state of the parent cation Bz⁺, with increasing fluorination. The minimum energy curve crossing goes up in energy accordingly and the B–X/A population transfer consequently slows down and becomes inefficient.^75–77 This is in nice accord with the observed fluorescence behaviour of these systems: whereas the parent cation Bz⁺, as mentioned above, is non-fluorescent just as is the monofluoro derivative,^71,72 one of the difluoro isomers displays weak fluorescence.⁷¹ Derivatives with 3 or more F-atoms emit strongly and their spectroscopic properties hae been explored in depth experimentally long ago.^71,78

Fig. 9 Schematic of the studied dynamics. Schematic overview of the lowest eight electronic component states of the benzene cation, depicted as potential energy V in eV as a function of a dimensionless effective nuclear coordinate Q_eff. This figure is taken from ref. 62.

Fig. 10 The yield from the adiabatic state calculation is convoluted with the Gaussian instrument response function of τ_IRF = 10 fs and compared to the experimental fragment ion yield (red dots). For more details see text. Figure reproduced from ref. 63 with permission from Elsevier, copyright 2019.

Fig. 11 Correlation between the lowest ionization potentials of benzene and its mono-, di-, and trifluoro derivatives according to the adjusted IPs. The Hartree–Fock canonical orbitals are included for the E_1g and the E_1g-derived states. Figure reproduced from ref. 77 with permission from AIP, copyright 2012.

VIII. Low-energy motion and the geometric phase

The examples addressed so far represent situations where the conical intersection is accessible energetically and the singular derivatice couplings render the nuclear motion (strongly) nonadiabatic. For low-energy motion this is no longer the case, but even in this case the conical intersection may leave its fingerprint on the dynamics. When encircling the point of degeneracy in a closed loop, the electronic wavefunction changes sign^79,80 and so does the nuclear wavefunction because the whole molecular wavefunction must remain single valued. This has been termed the molecular Aharonov-Bohm effect and found wide attention also in the physics literature.^3,79,80 For example, it affects in a characteristic way the quantization pattern of the rotational component of the nuclear wavefunction for the E × e Jahn–Teller effect.⁸⁰

The above sign change admits a more general interpretation in terms of geometric phases. Within the adiabatic representation the electronic eigenstates depend parametrically on the nuclear coordinates R. This allows the definition of the Berry connection (or Mead–Truhlar vector potential)⁸⁰

A_n(R) = i〈ψ_n(R)|∇_Rψ_n(R)〉, (14)

which acts as an effective gauge potential on the nuclear motion. When the nuclei evolve along a closed path C in configuration space, the nuclear wavefunction acquires a geometric phase⁷⁹For a loop encircling a conical intersection this phase equals π, which explains the sign change of the adiabatic electronic and nuclear wavefunctions discussed above, as first recognized by Herzberg and Longuet–Higgins.⁸¹ The molecular Aharonov–Bohm effect can therefore be viewed as the consequence of a gauge structure induced by the electronic motion.

The local quantity associated with the Berry connection is the Berry curvature

Ω_n(R) = ∇_R × A_n(R), (16)

which plays the role of an effective magnetic field in nuclear configuration space.^79,82 A conical intersection acts as a singular source of this field, analogous to a magnetic monopole located at the degeneracy point.⁸³ Using Stokes' theorem, the geometric phase can be written as a surface integral of the Berry curvature over a surface S bounded by the loop C,

The total flux of the Berry curvature through a closed surface defines the Chern number

For an isolated conical intersection one obtains C_n = ±1, demonstrating that the associated geometric phase is topological in nature.^84,85 The phase change of π is therefore not model dependent but represents a topological invariant characterizing the degeneracy. Multiple intersections contribute additively to this invariant, which explains the cancellation of the sign change when an even number of conical intersections is enclosed, as discussed above for multi-state situations.

A noteworthy extension arises when encircling the point of degeneracy several (say n) times. Then the sign change acculumulates and the overall change amounts to (−1)ⁿ. A similar phenomenon occurs then several PE surfaces intersect.^3,79,80 Also then, each conically intersecting surface “contributes” a factor −1 and for three PE surfaces there is no such sign change. Likewise, for glancing intersections such as for the Renner–Teller effect⁸⁶ there is no such sign change. An example for the former case is provided by Fig. 13. This is taken from our study of the so-called B excited state of Na₃.⁸⁷ In the D_3h equilateral triangular reference geometry there are an E′ and an excited electronic states of Na₃.⁸⁷ Upon bending (i.e. along one component of the doubly degenerate e′ mode of equilateral Na₃) the degeneracy of the E′ state is lifted and one component repels from the state (see Fig. 4c of ref. 87). Neglecting the energy difference leads to the shifted parabolae of Fig. 13 where the dashed line represents the state. Including the pseudorotational motion leads to rotationally symmetric PE surfaces with cuts as drawn in Fig. 13. However, the dynamics with and without the nondegenerate state is entirely different: for the 2-state scenario the sign change occurring when encircling the conical intersection leads to destructive interference of the two parts of the wavepacket that have encircled the degeneracy point (at the origin) on opposite sides, and a nodal pattern at the other side than that of the initial wave-packet location results (right hand panels of Fig. 12). For the three-state scenario there is no destructive interference and no such nodal line, i.e. more a “usual” behaviour (left-hand panels of Fig. 12). The consequences for observables such as the fs pump–probe ionization spectrum have been explored.^88,89

Fig. 12 Contour plots of the wave-packets in the excited state for the JT and the PJT treatment: (a) t = O fs; (b) t = 209 fs; (c) t = 500 fs; (d) t = 1755 fs. Before the wave-packets interfere with themselves at ϕ = π (negative Q_y, axis), they appear equal. For long times the difference between the E and E + A systems becomes significant. Figure reproduced from ref. 88 with permission from Elsevier, copyright 1994.

Fig. 13 Potential curves for the JT and PJT cases in C_2v symmetry. The dashed line appears only in the PJT case. The horizontallines mark the (central) energies of the wavefunction in the ground, excited and ionized state and the pulse transitions are sketched with arrows. The energy width for the 60 fs pulses is 124 cm⁻¹. Figure reproduced from ref. 88 with permission from Elsevier, copyright 1994.

We finally mention the intriguing phenomenon that in certain situations the accompanying phase factor may not be ±1 but rather a general unimodular complex number. This has been demonstrated early in the literature⁹⁰ and becomes relevant in the presence of spin–orbit coupling. It may also occur for many-particle statistics, which thus represent neither Bosons nor Fermions, but rather Anyons (as coined by Wilczek⁹¹). Examples of this type have indeed been identified for excitations in low-dimensional systems.⁹² From the geometric viewpoint these phases arise from the Berry connection associated with adiabatic transport in parameter space; the familiar ±1 sign change at an isolated conical intersection represents the special case where the integrated Berry curvature yields a Chern number of unity, whereas more general situations allow arbitrary unimodular phases.

IX. Two system types with avoided crossings of PE surfaces

At the beginning we have pointed out the (incomplete) similarity between intersections and weakly avoided crossings of PE surfaces. Thus it seems natural to address also the latter scenario at least briefly. Two types of such systems have been considered by us and are mentioned here.

A. Criegee intermediates

Carbonyl oxides or Criegee intermediates originate from ozonolysis of alkenes and are a significant fraction of volatile organic compounds in the Earth's lower atmosphere.⁹³ Their UV photodissociation typically occurs through breaking the O–O bond when moving on the coupled X¹A′, B¹A′ and C¹A′ electronic states.^46,94–97 The situation is illustrated by the PE drawings of Fig. 14 which also display an avoided crossing occurring near an O–O distance of ∼1.7 Å. For the smallest system (formaldehyde oxide) the avoided crossing has been investigated closer, but always a finite energy gap was obtained.⁹⁵ Although a degeneracy of the two states cannot be rigorously excluded, it would seem to have no major impart on the dynamics, that is, model simulations with a finite gap allowed to reproduce the experimental UV absorption spectra well, see Fig. 15. The intramolecular couplings lead to a rich population dynamics, and an increasing B–C coupling reduces the lifetime of the Criegee intermediate (see Fig. 14 and 18 of ref. 94): the optical excitation leads to the bound part of the B-state PE curve and the dissociation becomes possible by the B–C interaction (Fig. 14). We conclude by mentioning that spin–orbit coupling has also been included in the theoretical treatment and been shown to lead to new, low-energy dissociaton channels.^94,97

Fig. 14 Adiabatic PECs of the Criegee intermediates along the O–O bond length for the three states X¹A′ (black), B¹A′ (green) and C¹A′ (cyan). The vertical excitation energy, minimum energy of the B state, energy at the avoided crossing and the energy at the barrier maximum are shown by red, blue, black and gray circles, respectively. In the insets, the energy gaps at the avoided crossing geometries are shown. The other molecular coordinates are fixed at the equilibrium structure values. This figure is taken from ref. 94.

Fig. 15 The computed excitation spectrum of CH₂OO for the B¹A′–X¹A′ transition compared with the available experimental and theoretical results. This figure is taken from ref. 94.

B. Weakly coupled symmetry-equivalent systems

In a different line of research we have analyzed the electronic spectra and the so-called quenching of the excitonic splitting in several examples of H-bonded molecular homodimers, such as in the aminopyridin and o-cyanophenol dimers and others.^98,99 In these and other cases, the electronic excitation would lead to degenerate pairs of locally excited dimers, but the degeneracy is lifted by the excitonic interaction which leads to delocalized, symmetry adapted excited states of the homodimer. In terms of delocalized, symmetry adapted vibrational modes one deals with symmetric and antisymmetric combinations of monomeric modes which may be excited in the spectrum. For the relevant system parameters⁹⁸ a symmetric double-well potential results (see Fig. 16) where the quenched S₁–S₂ excitonic splitting emerges as the tunneling splitting in the double-well potential.⁹⁸ This is 1–2 orders of magnitude smaller than the vertical splitting of the two PE curves at the origin, which gives a vivid picture of the origin of the reduction or quenching; in more physical language this would be termed an S₁–S₂ splitting “dressed” by phonons (vibrational excitation).

Fig. 16 Schematic drawing of adiabatic potential energy curves (full lines) indicating the vertical and quenched excitonic splitting Δ_vert (=Δ_el) and Δ_vibron, respectively. The dashed lines near the avoided crossing denote the diabatic curves, corresponding to a localized excitation. Figure reproduced from ref. 98 with permission from AIP, copyright 2012.

There is quite a variety of further aspects which is beyond the scope of this short article.⁹⁸ We rather proceed by addressing a seemingly quite different phenomenon which can nevertheless be treated within a similar theoretical framework, namely, dynamic core–hole localization in systems with several equivalent core-hole sites. Examples with two different sites are the C 1s vacancy in C₂H₄ and C₂H₂ and the O 1s vacancy in CO₂ (and also the analogous core–valence excitations).^100–102 The gerade and ungerade linear combinations of localized core-holes differ energetically by ∼0.02 eV or less. Upon displacement along an asymmetric C–H or C–O stretching mode the equivalence of the core-hole sites is lost and the energy gap of the two cores holes correspondingly increases (strongly). The qualitative picture resembles that of Fig. 16, but with an adapted meaning of the vibrational coordinate. Note that for the core–hole problem one would normally start from a delocalized decription, whereas for the exciton coupling problem a localized approach is the more natural starting point. Re-writing the core–hole problem in a localized basis, casts the core–hole energy difference into the off-diagonal Hamiltonian matrix element where it adopts the role of the excitonic coupling for this other type of systems. Concerning dynamics, the smallnes of the core–hole energy gaps leads to motion on either of the two equivalent left-hand or right-hand shifted parabolae indicated by the dotted lines in Fig. 16. The has been termed “dynamic core hole localization”.^100–102 In the localized, diabatic basis the two PE curves are (approximately) decoupled and the “adiabatic” separation of electronic and nuclear motions is restored, although with a re-defined electronic basis. In this sense the dynamics differs from that addressed earlier in this paper. The vibrational structure in the electronic spectra is regular, but with excitation of asymmetric vibrational modes in odd quanta.^100,101

X. Concluding remarks

In this short review paper we have presented a compact overview over our work on quantum molecular dynamics in the past years and few decades. Time-dependent and time-independent aspects have been covered and high-energy and low-energy regimes been briefly contrasted. Quite naturally, perspectives for future work emerge from such studies, but also open issues and unsolved problems remain. Here we mention briefly two of them.

An intriguing phenomenon is the damping of seemingly coherent oscillations in conical intersection situations with at least three vibrational modes. Whereas for two modes the electronic populations display an arratic, fluctuation-type behaviour, for three degrees of freedom this behaviour changes profoundly and is “replaced” by a seemingly regular, but damped oscillatory behaviour. Note that this happens for an isolated quantum system without coupling to a heat bath or other type of environment. The molecule seems to act as its own heat bath. Only three degrees of freedom are needed for this phenomenon to appear. An intuitive picture underlying the sitaution is missing to our knowledge, as is a theory to obtain the timescale of damping.

An even more counterintuitive situation prevails in the multi-mode and multi-state system of Bz⁺ addressed briefly above. When preparing the system in the highest-energy Ẽ state like in the HHG pump–probe experiment of ref. 62, it undegoes stepwise fs population transfer processes to the lowest excited B̃ state. This is in accord with earlier results and general expectation. But the system gets trapped in the B̃ state without energetic barriers or other mechanisms identified so far. Various possibilities are conceivable and have been proposed, but an answer has to await future work.

More generally, future studies will involve aspects like neural networks, artificial intelligence and similar modern techniques. Neural networks have already been involved in the construction of quasidiabatic states by Yarkony and colleagues. The analysis of complex photochemcial reaction network will benefit from an increasing number of multi-mode multi-state approaches that became available recently. It is hoped that all these developments will contribute to a rich future of the field.

Conflicts of interest

There are no conflicts to declare.

Data availability

No new data were generated in this tutorial review.

Appendices

A. The Born–Oppenheimer approximation

Schrödinger equation for coupled electronic and nuclear motions:

H = H_el + T_N (A1)

H_el = T_e + U(x,Q) (A2)

H_elϕ_n(x, Q) = V_n(Q)ϕ_n(x, Q) (assume solved) (A3)

HΨ(x, Q) = EΨ(x, Q) (A4)

x and Q denote the sets of electronic and nuclear coordinates, respectively. Correspondingly ϕ and χ stands for the electronic and nuclear wave functions.

1. Derivation of the coupled equations. For simplicity, put

By putting:

Λ_nm = 0 (A13)

→[T_N + V_n(Q) − E]χ_n(Q) = 0 (A14)

It follows:

• (Electronic) eigenvalues, V_n(Q), of a given state correspond to the potential energy hypersurface for the nuclear motion.

• Total molecular wavefunction becomes a product of a nuclear and electronic wave function:

Ψ(x, Q) = χ_n(Q)ϕ_n(x, Q) (A15)

• Valid, e.g., when ϕ_n(x, Q) ≈ ϕ_n(x − Q).

Electrons follow the nuclear motion instantaneously (adiabatic), due to the large ratio between nuclear and electronic masses (i.e. the large effective mass of a nucleus compared to that of an electron M_i ≫ m_el).

2. Hellmann–Feynman relation. Re-writing the non-adiabatic (derivative) coupling terms:

Multiplying from the left by and integrating over the electronic coordinates, x, leads to:

In the vicinity of a degeneracy the derivative couplings can diverge and the adiabatic approximation is expected to break down.

B. The quadratic E × e Jahn–Teller effect

Consider the well-known Hamiltonian for the E × e Jahn–Teller (JT) effect, expanded to second order in the nuclear displacements.^19,20 Adopting polar coordinates ρ, ϕ and denoting the corresponding 2D harmonic oscillator Hamiltonian by H₀, one hasHere the rows and columns represent the complex components of the doubly-degenerate electronic wavefunction, 1 stands for the 2 × 2 unit matrix, and κ and g are the first-order and second-order JT coupling constants, respectively. H₀ is the Hamiltonian of the two-dimensional harmonic oscillator with kinetic energy T_N and potential energy (ω_e/2)ρ². The PE matrix is diagonalized to yield the adiabatic PE surfaces as follows.²³ Definingthe Hamiltonian can be written asand the potential energy part is easily seen to be diagonalized by the transformationThis yields the diagonal matrix of adiabatic PE surfaces as follows:Within the RDS scheme the diagonal matrix of adiabatic PE surfaces is “diabatized” by transforming back with the eigenvector matrix of the linear Exe JT effect, i.e.

It follows that (see ref. 23 for details)

Here the unperturbed harmonic oscillator Hamiltonian H₀ is assumed to be of the same (harmonic oscillator) form as in the original representation of eqn (B1). This represents the basic dynamical approximation underlying the regularization (RDS) scheme. Further analysis shows that the terms neglected when using the regularized quasidiabatic states remain finite at the origin ρ = 0 (point of CoIn),²³ whereas those appearing in the adiabatic basis diverge.^19,20,23 Thus, using eqn (B7) should yield considerably better results than obtained in the adiabatic approximation, i.e replacing H byIn the numerical tests displayed in Fig. 1 and 2 the full solution of the Hamiltonian, eqn (B1), is compared with that of the RDS result, eqn (B7). The results are indeed so close that they are easily overlooked in an inspection by eye. This underlines the high quality of the results obtained with the regularized diabatic states. The results of employing the widespread linear vibronic coupling scheme are displayed in the insets of the figures. They differ considerably from those of the full (second-order coupling) calculation which illustrates impact of the second-order coupling terms. Without the second-order terms the RDS scheme would be exact. The comparison illustrates that also the impact of the second-order terms on the coupled dynamics is very well covered by the regularization scheme.

C. Diabatization for an in-plane seam

The 1¹A₂ and 1¹B₁ excited states of SO₂ determine the gross features of its visible absorption spctrum.²⁴ The corresponding adiabatic PE surfaces feature a strong relative shift along the bending angle and an only very moderate displacement along the symmetric S–O stretching coordinate, see Fig. 5. Construction of quasidiabatic states is needed for a general set of coordinates R_s = R₁ = R₂ and ϕ with Q_u = 0 (C_2v point group) and can proceed based on two different alternatives. In one case the given C_2v coordinate set (R_s and ϕ) is projected onto the seam and the determination of λ follows the same basic procedure as in the 2D case of Appendix B by determining the gradient of the potential energy difference on the seamIn the other case λ is obtained by a generalization well known from LVC theory leading to ref. 5 and 6In apparent nomenclature the results for λ are labeled as λ_seam and λ_C2v, respectively. We refer the reader to the literature for details and proceed by comparing the results in Fig. 3. Whereas the contour lines are general closed curves for λ_C2v, they are vertical straight lines for λ_seam by the way the projection is performed. Crucially, the contour lines for a given value of λ cross on the seam, that is, while the numerical values of λ are generally quite different, they are idential along the seam! This is most important since along the seam the derivative coupling terms λ diverge while elsewhere they are finite. In other words, while the derivative coupling terms are generally different for λ_C2v and λ_seam, they coincide along the seam, i.e. where they diverge. Both variants have been employed in determining the visible absorption spectrum of SO₂ and the results are compared in the lower panel of Fig. 4. They are almost indistinguishable in the lower-energy part of the spectrum and differ by a very minor, nearly constant amount in the higher-energy part. Thus the focus on the singular derivative couplings underlying the RDS scheme is expected to lead to very reliable results in many cases.

Acknowledgements

We dedicate this paper to Prof. N. Sathyamurthy's 75th birthday with a Festschrift in Phys. Chem. Chem. Phys. on the structure and dynamics of chemical systems. It is a pleasure to acknowledge a fruitful collaboration with many friends and colleagues over the years, especially L. S. Cederbaum, W. Domcke, S. Faraji, E. Gromov and S. Mahapatra. We are indebted to the Deutsche Forschungsgemeinschaft (DFG) for financial support (KO 945/23-1).

References

1. E. Teller, Isr. J. Chem., 1969, 7, 227.
2. C. Zener, Proc. R. Soc. London, 1932, 137, 696.
3. H. C. Longuet-Higgins, U. Öpik, M. H. L. Pryce and R. A. Sack, Proc. R. Soc. A, 1958, 244, 1.
4. W. Moffitt and W. Thorson, Phys. Rev., 1957, 108, 1251.
5. H. Köppel, W. Domcke and L. S. Cederbaum, Adv. Chem. Phys., 1984, 57, 59.
6. W. Domcke, D. R. Yarkony and H. Köppel, Conical Intersections: Electronic Structure, Dynamics and Spectroscopy, World Scientific, New Jersey, 2004.
7. F. Bernardi, M. Olivucci and M. A. Robb, Chem. Soc. Rev., 1996, 25, 321.
8. D. R. Yarkony, Acc. Chem. Res., 1998, 31, 511.
9. B. G. Levine and T. J. Martínez, Annu. Rev. Phys. Chem., 2007, 58, 613.
10. G. A. Worth and L. S. Cederbaum, Annu. Rev. Phys. Chem., 2004, 55, 127.
11. in Conical Intersections: Theory, Computation and Experiment, ed. W. Domcke, D. R. Yarkony, and H. Köppel, Advanced Series in Physical Chemistry, World Scientific Publishing Co., Singapore, 2011, vol. 17, p. 754.
12. J. von Neumann and E. P. Wigner, Phys. Z., 1929, 30, 467.
13. A. Viel and W. Eisfeld, J. Chem. Phys., 2004, 120, 4603.
14. W. Eisfeld and A. Viel, J. Chem. Phys., 2005, 122, 204317.
15. C. A. Mead and D. G. Truhlar, J. Chem. Phys., 1982, 77, 6090.
16. M. Baer, Chem. Phys. Lett., 1975, 35, 112.
17. H. Köppel and B. Schubert, Mol. Phys., 2006, 104, 1069.
18. H. Köppel, J. Gronki and S. Mahapatra, J. Chem. Phys., 2001, 115, 2377.
19. R. Englman, The Jahn-Teller Effect in Molecules and Crystals, Wiley, 1972, wiley Series.
20. I. Bersuker, The Jahn-Teller Effect, Cambridge University Press, 2006, pp. 12–44.
21. L. Cederbaum, W. Domcke, H. Köppel and W. Von Niessen, Chem. Phys., 1977, 26, 169.
22. C. Cattarius, G. A. Worth, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 2001, 115, 2088.
23. A. Thiel and H. Köppel, J. Chem. Phys., 1999, 110, 9371.
24. C. Lévêque, A. Komainda, R. Taeb and H. Köppel, J. Chem. Phys., 2013, 138, 044320.
25. A. Komainda, B. Ostojić and H. Köppel, J. Phys. Chem. A, 2013, 117, 8782.
26. H.-D. Meyer, U. Manthe and L. Cederbaum, Chem. Phys. Lett., 1990, 165, 73.
27. M. Beck, A. Jäckle, G. Worth and H.-D. Meyer, Phys. Rep., 2000, 324, 1.
28. H.-D. Meyer, F. Gatti and G. Worth, Multidimensional Quantum Dynamics: MCTDH Theory and Applications, WILEY-VCH Verlag GmbH and Co. KGaA, Weinheim, 2009.
29. M. Beck, A. Jäckle, G. Worth and H.-D. Meyer, MCTDH, The Heidelberg MCTDH Package Version 8.4.20, For the current version, see https://www.pci.uni-heidelberg.de/tc/usr/mctdh/doc/index.html.
30. H. Müller, H. Köppel and L. S. Cederbaum, J. Chem. Phys., 1994, 101, 10263.
31. S. Mahapatra, H. Köppel, L. S. Cederbaum, P. Stampfuß and W. Wenzel, Chem. Phys., 2000, 259, 211.
32. S. Mahapatra, H. Köppel and L. S. Cederbaum, J. Chem. Phys., 1999, 110, 5691.
33. E. Haller, H. Köppel and L. S. Cederbaum, J. Mol. Spectrosc., 1985, 111, 377.
34. C. Lévêque, D. Pelaez, H. Köppel and R. Taieb, Nat. Commun., 2014, 5, 4126.
35. C. Lévêque, H. Köppel and R. Taieb, J. Chem. Phys., 2014, 140, 204303.
36. C. Lévêque, R. Taieb and H. Köppel, J. Chem. Phys. Commun., 2014, 140, 091101.
37. C. Lévêque, R. Taieb and H. Köppel, Chem. Phys., 2015, 460, 135.
38. Y. Hamada and H. J. Merer, Czech. J. Phys., 1975, 53, 2555.
39. A. F. Douglas, Czech. J. Phys., 1958, 36, 147.
40. E. B. Wilson, J. C. Decius and P. Cross, Molecular vibrations: the theory of infrared and Raman vibrational spectra, McGraw-Hill, New York City, 1955.
41. J. H. Frederick and C. Woywod, J. Chem. Phys., 1999, 111, 7255.
42. A. J. Musser, M. Maiuri, D. Brida, G. Cerullo, R. H. Friend and J. Clark, J. Am. Chem. Soc., 2015, 137, 5130.
43. S. M. Menke and R. J. Holmes, Energy Environ. Sci., 2014, 7, 499.
44. S. Santra, J. Ray and D. Ghosh, J. Phys. Chem. Lett., 2022, 13, 6800.
45. A. Komainda, D. Lefrancois, A. Dreuw and H. Köppel, Chem. Phys., 2017, 482, 27.
46. B. Nikoobakht, R. Hakim, M. F. S. J. Menger and H. Köppel, Mol. Phys., 2022, 121, e2132186.
47. M. F. S. J. Menger and H. Köppel, J. Phys. Chem. A, 2023, 127, 8501.
48. A. Komainda, B. Ostojić and H. Köppel, J. Phys. Chem. A, 2013, 117, 8782.
49. W. J. Glover and T. J. Martinez, et al., J. Chem. Phys., 2018, 148, 164303.
50. D. G. Leopold, R. D. Pendley, J. L. Roebber, R. J. Hemley and V. Vaida, J. Chem. Phys., 1984, 81, 4218.
51. V. Vaida, Acc. Chem. Res., 1986, 19, 114.
52. A. Komainda, I. Lyskov, C. M. Marian and H. Köppel, J. Phys. Chem. A, 2016, 120, 6541.
53. A. Komainda, A. Zech and H. Köppel, J. Mol. Spectrosc., 2015, 311, 25.
54. I. Lyskov, H. Köppel and C. M. Marian, Phys. Chem. Chem. Phys., 2017, 19, 3937.
55. A. Chattopadhyay, M. F. S. J. Menger, D. S. Sisodiya and H. Köppel, J. Phys. Chem. A, 2025, 129, 2442.
56. Y. Mo, H. Jiao and P. R. Schleyer, J. Org. Chem., 2004, 69, 3493.
57. H. Köppel, E. V. Gromov and A. B. Trofimov, Chem. Phys., 2004, 304, 35.
58. E. V. Gromov, A. B. Trofimov, N. M. Vitkovskaya, H. Köppel, J. Schirmer, H.-D. Meyer and L. S. Cederbaum, J. Chem. Phys., 2004, 121, 4585.
59. M. A. Parkes and G. A. Worth, J. Chem. Phys., 2024, 161, 114305.
60. E. V. Gromov, V. S. Reddy, F. Gatti and H. Köppel, J. Chem. Phys., 2013, 139, 234306.
61. S. Faraji, M. Vazdar, V. S. Reddy, M. Eckert-Maksić, H. Lischka and H. Köppel, J. Chem. Phys., 2011, 135, 154310.
62. M. C. E. Galbraith, S. Scheit, N. V. Golubev, G. Reitsma, N. Zhavoronkov, V. Despré, F. Lépine, A. I. Kuleff, M. J. J. Vrakking, O. Kornilov, H. Köppel and J. Mikosch, Nat. Commun., 2017, 8, 1018.
63. S. Scheit, S. Goswami, H.-D. Meyer and H. Köppel, Comput. Theor. Chem., 2019, 1150, 71.
64. R. Lindner, K. Müller-Dethlefs, E. Wedum, K. Haber and E. R. Grant, Science, 1996, 271, 1698.
65. J. G. Goode, J. D. Hofstein and P. M. Johnson, J. Chem. Phys., 1997, 107, 1703.
66. K. Walter, R. Weinkauf, U. Bösl and E. W. Schlag, Chem. Phys. Lett., 1989, 155, 8.
67. A. B. Burrill, Y. K. Chung, H. A. Mann and P. M. Johnson, J. Chem. Phys., 2004, 120, 8587.
68. J. Eiding, R. Schneider, W. Domcke, H. Köppel and W. V. Niessen, Chem. Phys. Lett., 1991, 177, 345.
69. H. Köppel, L. S. Cederbaum and W. Domcke, J. Chem. Phys., 1988, 89, 2023.
70. H. Köppel, Chem. Phys. Lett., 1993, 205, 361.
71. M. Allan, J. P. Maier and O. Marthaler, Chem. Phys., 1977, 26, 131.
72. O. Braitbart, E. Castellucci, G. Dujardin and S. Leach, J. Phys. Chem., 1983, 87, 4799.
73. M. Döscher and H. Köppel, Chem. Phys., 1997, 225, 93.
74. I. Baldea, J. Franz and H. Köppel, J. Mol. Struct., 2007, 838, 94.
75. S. Faraji and H. Köppel, J. Chem. Phys., 2008, 129, 074310.
76. S. Faraji, H.-D. Meyer and H. Köppel, J. Chem. Phys., 2008, 129, 074311.
77. S. Faraji and H. Köppel, J. Chem. Phys., 2012, 137, 22A531.
78. C. Cossart-Magos, D. Cossart and S. Leach, Mol. Phys., 1979, 37, 793.
79. M. V. Berry, Proc. R. Soc. A, 1984, 392, 45.
80. C. A. Mead and D. G. Truhlar, J. Chem. Phys., 1979, 70, 2284.
81. G. Herzberg and H. C. Longuet-Higgins, Discuss. Faraday Soc., 1963, 35, 77.
82. M. Baer, Beyond Born–Oppenheimer: Electronic Nonadiabatic Coupling Terms and Conical Intersections, Wiley, Hoboken, 2006.
83. J. Moody, A. Shapere and F. Wilczek, Phys. Rev. Lett., 1986, 56, 893.
84. L. S. Cederbaum, J. Zobeley and H. Köppel, J. Chem. Phys., 1997, 106, 661.
85. I. G. Ryabinkin, L. Joubert-Doriol and A. F. Izmaylov, Chem. Rev., 2017, 117, 9865.
86. R. Renner, Z. Phys., 1934, 92, 172.
87. F. Cocchini, T. H. Upton and W. Andreoni, J. Chem. Phys., 1988, 88, 6068.
88. J. Schön and H. Köppel, Chem. Phys. Lett., 1994, 231, 55.
89. J. Schön and H. Köppel, J. Phys. Chem. A, 1999, 103, 8579.
90. A. J. Stone, Proc. R. Soc. London, Ser. A, 1976, 351, 141.
91. F. Wilczek, Phys. Rev. Lett., 1982, 49, 957.
92. S. Dhar, B. Wang, M. Horvath, A. Vashisht, Y. Zeng, M. B. Zvonarev, N. Goldman, Y. Guo, M. Landini and H.-C. Nägerl, Nature, 2025, 642, 53.
93. M. A. H. Khan, C. J. Percival, R. L. Caravan, C. A. Taatjes and D. E. Shallcross, Environ. Sci.: Processes Impacts, 2018, 20, 437.
94. B. Nikoobakht, M. F. S. J. Menger and H. Köppel, Phys. Chem. Chem. Phys., 2025, 27, 20433.
95. B. Nikoobakht and H. Köppel, Mol. Phys., 2021, 119, e1958019.
96. B. Nikoobakht and H. Köppel, Phys. Chem. Chem. Phys., 2022, 24, 12433.
97. B. Nikoobakht and H. Köppel, Phys. Chem. Chem. Phys., 2024, 26, 24591.
98. S. Kopec, P. Ottiger, S. Leutwyler and H. Köppel, J. Chem. Phys., 2012, 137, 184312.
99. P. Ottiger, H. Köppel and S. Leutwyler, Chem. Sci., 2015, 6, 6059.
100. B. Kempgens, H. Köppel, A. Kivimäki, M. Neeb, L. S. Cederbaum and A. M. Bradshaw, Phys. Rev. Lett., 1997, 79, 3617.
101. H. Köppel, F. X. Gadea, G. Klatt, J. Schirmer and L. S. Cederbaum, J. Chem. Phys., 1997, 106, 4415.
102. N. V. Dobrodey, H. Köppel and L. S. Cederbaum, Phys. Rev. A: At., Mol., Opt. Phys., 1999, 60, 1988.

---