Solving Vibronic Dynamics in Electron Continuum

TL;DR

The paper introduces a model for studying electron-molecule collisions involving metastable states and compares methods for solving resonance dynamics.

Contribution

A novel two-dimensional model of conical intersections and comparison of Krylov-subspace solvers for electron continuum dynamics.

Findings

Two Krylov-subspace methods with different preconditioning were compared for resonance dynamics.

One method was successfully applied to a CO2 vibrational excitation model with four vibrational degrees of freedom.

Electron energy-loss spectra from the models were analyzed to assess scattering outcomes.

Abstract

We present a general two-dimensional model of conical
intersection
between metastable states that are vibronically coupled not only directly
but also indirectly through a virtual electron in the autodetachment
continuum. This model is used as a test ground for the design and
comparison of iterative solvers for resonance dynamics in low-energy
electron–molecule collisions. Two Krylov-subspace methods with
various preconditioning schemes are compared. To demonstrate the applicability
of the proposed methods on even larger models, we also test the performance
of one of the methods on a recent model of vibrational excitation
of CO2 by electron impact based on three vibronically coupled
discrete states in continuum (Renner–Teller doublet of shape
resonances coupled to a sigma virtual state) including four vibrational
degrees of freedom. Two-dimensional electron energy-loss spectra resulting
from electron–molecule scattering within the models are briefly
discussed.

Figures

Sections through the potential energy surfaces in the qg = 0 (left) and qu = 0 (right) planes for the ECD86 model. Blue-shaded areas give the position and width of the fixed nuclei electronic resonance. See the text for more details.

Perspective view of the potential energy manifold for the ECD86 model. The width (inverse lifetime) is marked by the color scale.

Sections through the model potentials in the qg = 0 (left) and qu = 0 (right) planes for the new model. See Figure 1 for more details.

Perspective view of the potential energy manifold for the new model, colored by the width of the resonance.

Convergence of the GMRES method for the model of Estrada et al.33 for various preconditioners. Number of iterations needed for each electron energy (top) and convergence of residuum for ϵi = 2 and 4 eV (bottom two panels). The top part shows results for both right-hand sides μ = e, o of the linear system; the bottom two parts display just μ = e.

Convergence of the COCG method for the model of Estrada et al.;33 see also the caption of Figure 5 for details.

Convergence of the GMRES method for the generalized model; see also the caption of Figure 5 for details.

Convergence of the COCG method for the generalized model; see also the caption of Figure 5 for details.

Number of iterations needed to solve the Schrödinger equation for the e + CO2 system using the COCG method without and with preconditioning. In the latter case, the curve is multiplied by a factor of 10.

2D electron energy-loss spectrum for the ECD86 model (top) and its decomposition to the gerade and ungerade symmetries (bottom). The intensity of spectrum is shown in a logarithmic scale.

Same data as in Figure 10 but for the new model. The black dotted lines in the top figure indicate a triangle containing the feature discussed in the text. This feature extends between the horizontal and diagonal lines almost all the way toward zero energy loss.

Integral cross sections for the elastic channel and the excitation of the two vibrational modes. The dashed curves are the data for a model with the off-diagonal coupling though the continuum switched off.

Tables

Table 1. Values of the Parameters Describing the Pseudo-Jahn–Teller Model of Estrada et al.33,a

parameter	value	parameter	value
ωg	0.258	ωu	0.091
E1	2.45	E2	2.85
κ1	–0.212	κ2	0.254
λ	0.318
a1	0.086	a2	0.186
b1	0.833	b2	0.375
l1	2	l2	1

Table 2. Values of the Parameters Describing the Generalization of the Modela

parameter	value	parameter	value
a1e	0.07	a2e	0.1
b1e	0.25	b2e	0.5
l1e	0	l2e	0
a1o	0.186	a2o	0.15
b1o	0.375	b2o	0.8
l1o	1	l2o	1
λ1	0.2	λ2	0.1

Full text

Introduction

Despite a long history of investigations (see for example reviews^1−5^), the collisions of low-energy electrons with molecules still represent a fascinating and challenging field of study. By low energy we mean here the energy below the electronic-excitation threshold, i.e., the energy that does not exceed few units of electronvolts. Even these low energies lead to many interesting phenomena like the appearance of sharp structures in cross sections^5,6^ or the possibility to select dissociation into different anionic fragments by tuning the energy.^7,8^ This topic is both interesting for practical applications^9^ and challenging for the theory even for small polyatomic molecules.^10−16^

In this paper, we will focus on the process of vibrational excitation in collision of an electron e^–^ with a molecule M initially in a vibrational state |vi⟩

$$[eqn]$$

mediated by one or several metastable anion states M^–^. After the process, the molecule is left in the final vibrational state |v_f_⟩. The total energy E during the collision is conserved

$$[eqn]$$

where ϵ are electron energies and E_v_ the energies of vibrational states of the molecule for the initial and final states before and after the collision. This process is closely related to the process of photodetachment of an electron e^–^ from a molecular anion M^–^

$$[eqn]$$

with initial energy E of the system defined now by the energy of the photon γ shone on the anion to excite it to the state (M^–^)*. The dynamics of both of these processes is driven by potential energy of states of the negative molecular ion and their widths for decay into electronic continuum channels.^17^ The energies of the released electrons are sensitive to the relative position of the anion and the neutral molecular states, and the selection rules are different than for the radiative transitions.^18−20^

The goal of this paper is to advance the detailed theory of the dynamics of electron detachment from anions in such processes. In the development of the theoretical methods, we keep in mind the description of experiments that study in detail the energies of released electrons,^17,21,22^ and in particular, we calculate the two-dimensional electron energy-loss spectrum (2D EELS) for our model. The 2D electron loss spectroscopy was pioneered by Currell and Comer^23,24^ and further developed by Allan and collaborators.^25^ To date, a dozen of high-resolution spectra for different molecules have been measured^17,25−32^ but the detailed understanding of such spectra for polyatomic molecules is mostly lacking.

In this paper, we present and test a general scheme for solving the nuclear dynamics of the negative ion formed in the collision of an electron with a polyatomic molecule. The scheme is tailored for a class of models that are inspired by the pseudo-Jahn–Teller model of Estrada, Cederbaum, and Domcke^33^ with modifications meant to make it a more realistic model of real molecules. This approach combines a model of vibronic coupling of several anionic states expanded in low-order polynomials in vibrational coordinates close to equilibrium geometry of the neutral molecule with a projection-operator approach to include the interaction of the anion discrete states with the electronic continuum. The model is rather flexible in adding states and vibrational degrees of freedom and the present scheme has been used to produce the results in our previous works on CO_2_.^34−36^ In these papers, we did not explain the methods and their performance in detail, a gap that is meant to be filled by this work.

We start the “Theory” section by reviewing the projection-operator approach to the dynamics of vibrational excitation in electron collisions with molecules. We then proceed by reminding the model of Estrada et al.^33^ and propose its generalization by including the vibronic coupling through the electron continuum in addition to the direct vibronic coupling present in the original model. This section is concluded by explaining the representation of the wave function components and the Hamiltonian in a basis constructed from neutral vibrational states. In “Krylov-Subspace Iteration Methods“ section, we first briefly introduce (the details are in the Appendix) used iteration methods and preconditioning schemes, and then we discuss their performance for the models. The section “Discussion of Resulting Spectra for Test Models” is devoted to a brief description of the obtained 2D spectra for the models, and we conclude by summarizing the results in the “Conclusions” section.

Theory

The vibrational and resonance dynamics in electron–molecule collisions has been studied theoretically for a long time (see, for example, one of the review papers^2,4,37,38^). The direct brute-force approach is only tractable for small molecules^39^ or for small deformations.^40^ A number of approximate schemes have therefore been developed: Born approximation, adiabatic-nuclei approximation, zero/effective range, or semiclassical approaches. In the present work, we focus on the development of the numerical schemes for the projection-operator approach based on the existence of an intermediate anion state (or states) that is responsible for the coupling of the electronic and vibrational motion. The approach is often used in its approximate form—the local complex potential approximation, but it is known to fail in predicting interesting phenomena like Wigner cusps or vibrationally excited Feshbach resonances. The nonlocal approach is well-developed for diatomic molecules^4,41^ but the attempts to use it for polyatomic molecules are scarce (see for example the paper by Ambalampitiya and Fabrikant^42^). In addition to bringing more degrees of freedom, the polyatomic molecules also exhibit interesting features like vibronic coupling of resonances, conical intersections, and exceptional points.^43,44^ Here, we follow the work of Estrada, Cederbaum, and Domcke^33^ (ECD86) and extend it to more general form of model functions and vibronic coupling. We start by presenting basic formulas resulting from the projection-operator formalism (for comprehensive review of the approach, see the paper by Domcke^4^). Then, we narrow the model to two vibrational degrees of freedom and two vibronically coupled discrete states.

Nonlocal Model for Multiple Discrete States in Continuum

The main idea of the nonlocal discrete state in continuum model is the assumption that the coupling of the electronic and vibrational degrees of freedom in the electron–molecule collision is mediated by one or a few discrete states and after their removal from the electronic continuum using projection-operator formalism of Feshbach,^45^ the electronic basis consisting of the discrete states and the orthogonalized continuum is diabatic. The vibrational excitation or dissociative attachment then proceeds through capture in the discrete state.

We define the projection operator

$$[eqn]$$

as a sum over a set of discrete states |d⟩ and the complementary operator

$$[eqn]$$

projecting on the background continuum. The basis in the background part can be chosen as the states that solve the background scattering problem

$$[eqn]$$

Here, V0(q⃗) is the potential energy surface of the neutral molecule, i.e., the energy of the ground electronic state |Φ_0_⟩ as a function of the positions of the nuclei q⃗. Since we consider only low-energy electron scattering below the threshold for the electronic excitation of the molecule, the state Φ_0_ is fixed, and we will further omit it from the notation. The electron continuum states |ϵμ⟩ are thus uniquely described by the electron energy ϵ and some other quantum numbers collectively denoted by μ (typically angular momentum). All states |d⟩ and |ϵμ⟩ thus form an orthogonal basis

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The electronic Hamiltonian in the $[eqn]$ -space is described by a matrix

$$[eqn]$$

where all matrix elements depend on the molecular geometry, i.e., positions of nuclei q⃗. The diagonal elements Vd(q⃗) = V0(q⃗) + Udd(q⃗) represent the diabatic discrete-state potentials and the off-diagonal part Udd′(q⃗) the direct vibronic coupling among the states. The coupling between the discrete state |d⟩ and the continuum |ϵμ⟩ is described by the coupling elements

$$[eqn]$$

These elements represent the vibronic coupling1 between the discrete state and the continuum, and they also lead to the second-order vibronic coupling among the discrete states mediated by the continuum as described below.

This way, we parametrized the matrix elements given by eqs 6, 10, and 11 of the Hamiltonian $[eqn]$ for the electron scattering from the molecule for each fixed position of the nuclei q⃗ by functions V0(q⃗), Udd′(q⃗), and Vdϵ^μ^(q⃗). To describe the electron scattering from the molecule including the vibronic dynamics, we start from the definition of the vibrational states |v⟩ of the target neutral molecule

$$[eqn]$$

where TN is the kinetic-energy operator for the nuclei and v is a set of quantum numbers that uniquely determine the vibrational states with energy Ev. It can be shown (see for example the review of Domcke^4^) that the vibronic motion of the anion is described by the effective Hamiltonian

$$[eqn]$$

which is the matrix in the indices d, d′ and the operator in the space of vibrational degrees of freedom. In the equation above, H0 is the Hamiltonian operator for the vibrations of the molecule multiplied by unity matrix δ_dd′_ in the discrete-state indices, U is the matrix with the elements Udd′ defined above, and the operator F describes the dynamical coupling of the discrete-state space to the electronic continuum

$$[eqn]$$

where η is a positive infinitesimal. This operator is also a matrix in the discrete-state indices and a nonlocal operator in the nuclear coordinate q⃗.

The discrete-state contribution to the T-matrix for vibrational excitation by electron scattering in a continuum state |ϵ_i_ μ_i⟩ from the initial vibrational state vi_ to final state vf and leaving in continuum state |ϵ_f_ μ_f_⟩ is given by^33^

$$[eqn]$$

and is closely related to the integral cross section for the vibrational excitation event

$$[eqn]$$

Finally, to simulate the full 2D electron energy-loss spectra, we have to collect vibrational excitation cross sections for all accessible final states

$$[eqn]$$

where ρ(ϵ) is the resolution function of the spectrometer (simulated here with a Gaussian function with full width at half maximum equal to 10 meV, which is comparable to the values in the ref (25)). The energy loss Δϵvf__ = Evf__ – Evi__ = ϵ_i_ – ϵ_f_ in each term in eq 17 is fixed by the energy conservation. The function S(ϵ_i_, Δϵ) gives the full experimental information in the electron energy-loss spectroscopy except for the angular resolution that can also be included,^36^ but it is not of interest in the present paper.

Pseudo-Jahn–Teller Model of Estrada et al.

The model by Estrada et al.^33^ assumes a molecule with an Abelian group of symmetry. They consider two discrete states d = 1, 2 that transform according to different irreducible representations of the symmetry group and are coupled vibronically through a nontotally symmetric vibrational mode qu. They also consider excitation of another totally symmetric mode qg, so that the geometry of the molecule within the model is described by a vector q⃗ = (qg,qu). The symmetry then dictates the structure of matrix U

$$[eqn]$$

This is a completely general form of the dependence of matrix U on coordinates when the terms are restricted up to the first order in q⃗ and the symmetry requirements are taken into account. Similarly, we can expand the matrix of the discrete-state-continuum coupling. For simplicity, we consider only two partial waves |ϵμ⟩ for μ = e, o, representing one even and one odd linear combination of partial waves coupled to the discrete-state space. In principle, we could consider more partial waves but they could be decoupled from the problem by a unitary transformation, thus grouping partial waves into effective channels with the number of channels not exceeding the dimension of the $[eqn]$ -space.^46^ Estrada et al.^33^ considered the coupling matrix Vdϵ^μ^ independent of the nuclear coordinates (we are going to lift this restriction in the next section). The symmetry selection rules then forbid the coupling Vdϵ^μ^ between the different symmetry of discrete state d and partial wave μ. We therefore assume that the discrete state d = 1 has even symmetry as the μ = e partial wave and d = 2 has the symmetry of the odd partial wave μ = o. Using eq 14, we then see that only the diagonal matrix elements F11(ϵ) and F22(ϵ) of the level-shift operator are nonzero. They can be generated from their imaginary parts (widths)

$$[eqn]$$ $$[eqn]$$

by means of the integral transform $[eqn]$ defined as (compare eq 14)

$$[eqn]$$

This transform can be worked out analytically for the assumed form of the widths

$$[eqn]$$

(see Berman et al.^47^). To complete the model description, we must give the vibrational Hamiltonian H0 of the neutral molecule. The model simply assumes harmonic vibrations

$$[eqn]$$

The vibrational eigenstates |v⟩ satisfying eq 12 with this harmonic Hamiltonian can be numbered by two quantum numbers ν = (ng, nu) and the vibrational energies are given by standard harmonic oscillator formula

$$[eqn]$$

The numerical values of ω_i_ and the parameters defining direct coupling matrix U and discrete-state-continuum matrix Ve for the model studied in Estrada et al.^33^ and used here for testing are given in Table 1 (for dimensionless lengths and energies in units of eV). Note that several variants of the model were used by Estrada et al.^33^ Here, we study only the most complex form of the model with the values of parameters, as given in the table. To visualize the character of the model, we show the one-dimensional sections through the model potentials in Figure 1. The functions shown are V0(q⃗), Vd(q⃗) = V0(q⃗) + Udd(q⃗), and the local complex potential $[eqn]$ , obtained from the pole of the fixed-nuclei K matrix, which has to be located iteratively.^35^ The perspective view of the local complex potential Re Vloc colored by values of −2 Im Vloc is also shown in Figure 2.

Sections through the potential energy surfaces in the qg = 0 (left) and qu = 0 (right) planes for the ECD86 model. Blue-shaded areas give the position and width of the fixed nuclei electronic resonance. See the text for more details.

Perspective view of the potential energy manifold for the ECD86 model. The width (inverse lifetime) is marked by the color scale.

Table 1: Values of the Parameters Describing the Pseudo-Jahn–Teller Model of Estrada et al.33,a

Generalized Model with Vibronic Coupling with Continuum States

The vibronic model above assumes the most simple structure of the discrete-state-continuum coupling matrix Vϵ = {Vdϵ^μ^} with row index d and column index μ

$$[eqn]$$

where $[eqn]$ , $[eqn]$ , and V1ϵ^o^(q⃗) = V2ϵ^e^(q⃗) = 0. We go one step beyond the approximation of the coupling matrix by constant terms, and we expand the matrix to the first order in the normal vibrational coordinates. This generalization is very useful in the description of the interaction of resonances through the electronic continuum that is switched off in the equilibrium geometry but becomes nonzero with deformation as, for example, in the pyrrole molecule.^48^ This feature was also an important ingredient of the model for the CO_2_ molecule.^34,35^ We will further assume that the dependence of Vdϵ(q⃗) = f(ϵ)g(q⃗) on the electron energy ϵ and the normal coordinates is separable. Taking into account the symmetry of the system, we get

$$[eqn]$$

where the terms that couple a discrete state to the partial wave of different symmetry must be odd functions of qu. We see that half of the total number of 12 terms (up to first order in q⃗) in the coupling matrix are zero due to the symmetry. For the purposes of the testing of the numerical methods, we choose the same form of the energy dependence as in the original model

$$[eqn]$$

with the values of the parameters given in Table 2.

Table 2: Values of the Parameters Describing the Generalization of the Modela

Using eq 14, we see that the structure of the nonlocal level-shift operator F(E – H0) is much richer

$$[eqn]$$

where we used the integral transform (21) again. The ordering of the terms that depend on qi with respect to $[eqn]$ must be kept because we substitute the operator ϵ = E – H0, which does not commute with the normal coordinates. The off-diagonal terms like F12 and F21 that introduce coupling of the discrete states through the continuum have been studied before in the local complex potential approximation (see for example refs (49−51)) among resonances of the same symmetry. A coupling of resonances of different symmetry requires dependence on vibrational coordinates. To the best of our knowledge, our work is the first one to include the full nonlocal form.

The potentials for the new generalized model are visualized in Figures 3 and 4. Note that the structure of such conical intersections in continuum has been investigated by Feuerbacher et al.^43,44^ In accordance with their findings, the potential manifolds shown in Figures 2 and 4 do not intersect in a single point like regular conical intersections but in a line segment bounded by two exceptional points, where not only the real but also the imaginary part of the two potentials is identical. The form of our model as given by eq 27 is more general than the expansion investigated by Feuerbacher et al.^43,44^ because they studied a linear coordinate expansion of the width function Γ whereas we prescribe the linear expansion of the coupling matrix Vϵ, which is more natural for the subsequent treatment of the dynamics. In our case, the linear form of the coupling matrix also produces quadratic terms in widths in eq 27. When the quadratic terms are omitted, we recover the form used in Feuerbacher et al.^43,44^ However, we cannot omit these terms in the dynamics since it would distort the unitarity of the S matrix.

Sections through the model potentials in the qg = 0 (left) and qu = 0 (right) planes for the new model. See Figure 1 for more details.

Perspective view of the potential energy manifold for the new model, colored by the width of the resonance.

Numerical Representation of the Dynamics

For the numerical solution of the dynamics, we expand wave function components in the harmonic oscillator basis |ν⟩ = |ng, nu⟩ associated with the model Hamiltonian of the neutral molecule (23). We first rewrite eqs 15 and 16 for the cross section as

$$[eqn]$$

where we defined auxiliary wave functions |Φ_v_^μ^⟩ with components

$$[eqn]$$

and with energy according to the conservation law (2), i.e., ϵ = E – Ev. The anion wave function |Ψ_vi__^μi_^⟩ satisfies

$$[eqn]$$

In the harmonic oscillator basis, this equation represents a system of linear equations for unknown components of the discrete-state wave function

$$[eqn]$$

where we introduced a compound index α ≡ (d, ng, nu). Using this notation, the matrix of this system Aα,α′ reads

$$[eqn]$$

and the scalar product in eq 28 can be written as the sum over the components of Φ_α_^vfμf^ and Ψ_α_

$$[eqn]$$

We cut off the basis in each dimension, keeping the states |ng⟩ for ng = 0,1, ···, Ng–1 and |nu⟩ for nu = 0,1, ···, Nu–1. The states Ψ_α_ are thus represented by N = 2NgNu component vectors and A is a N × N matrix. Typical values used in the following calculations, Nu = 100, Ng = 50, result in converged cross sections, as we checked by varying these parameters.

For the solution of eq 30 in the original model, Estrada et al.^33^ devised a specially tailored method based on the block-tridiagonal structure of matrix A. Our aim in this paper is to develop a more general method capable of solving a larger class of models and test it both on the original model and on our generalization. Matrix A is large but sparse. From the character of the problem, it is also complex symmetric but not Hermitian. The structure of matrix A depends on the order of the basis vectors as illustrated in Figure A1 of the appendix for the generalized model.

Krylov-Subspace Iteration Methods

The Krylov-subspace iteration methods (see Appendix for details) are ideally suited for the solution of the system of equations with sparse matrices, since they are based on expansion of the solution in the basis produced by repeated multiplication of the right-hand side b with matrix A. This matrix multiplication can be coded very efficiently (see Appendix). In this section, we test two methods: COCG (Conjugate orthogonal conjugate gradient) method and GMRES (generalized minimal residual) method well-known in numerical linear algebra. The following convergence test of the two iteration methods shows the convergence of the residuum rn = b – Axn → 0 for the approximation xn of the solution Ψ_α_ in the nth iteration. We discuss the convergence for different methods in terms of the number of iterations needed to obtain the residuum on the order of 10^–5^. The computational time is not given since it depends on details of both hardware and software tools used in the calculation, but as a rule of thumb, we can say that we need tens of minutes on desktop PC to get one 2D spectrum consisting of hundreds of energies in the case of the ECD86 model and its generalization but a computational cluster was needed for the calculation of the spectrum for CO_2_. In that case, solving the linear system took up to 10 h for one initial electron energy on a single CPU core and, in addition, a large amount of memory (∼20 GB) was necessary to store the matrix of the preconditioner. Note that the calculation can be easily parallelized over the initial energies since the claculations for different initial energies are independent of each other.

Before discussing the individual tests, we would like to point out that the essential ingredient of the iteration methods is the preconditioning. It is based on modification of matrix A → M^–1^A(M^T^)^–1^ (see Appendix), where M is easily invertible approximation to A. In the following, we discuss the preconditioning by matrix M consisting of the block-diagonal part of A for different ordering of basis vectors. We show small blocks M^dg^, M^du^, and M^gu^ of the sizes Nu, Ng, and 2, respectively, and the larger blocks M^d^, M^g^, and M^u^ of sizes NgNu, 2Nu, and 2Ng. For technical details, we refer again to the Appendix.

Numerical Testing

We applied the methods GMRES and COCG described above to solve system (30) with matrix (32) for the original ECD86 model and for our generalization of the model. The performance of each method for different preconditioning is discussed separately for the two models in the next two paragraphs. The last paragraph also discusses the performance of the COCG method for a realistic model that describes inelastic electron scattering from the CO_2_ molecule. Note that the convergence properties (number of iterations) are independent of the size of the basis (Nu, Ng) used once the size is large enough to achieve the convergence, although the computational costs of the single iteration depend on these parameters.

Performance of the Methods for the ECD86 Model

The performance is shown in Figures 5 and 6. Each of the figures is devoted to one of the methods comparing different preconditioning schemes. The top graph summarizes the number of iterations needed for convergence for all energies, and bottom two graphs demonstrate the decrease of the residuum norm for two selected energies ϵ_i_ = 2 and 4 eV. The different preconditioning methods are shown with different colors. The curves of the same color correspond to two different right-hand sides, μ_i_ = o, e in eq 30.

Convergence of the GMRES method for the model of Estrada et al.33 for various preconditioners. Number of iterations needed for each electron energy (top) and convergence of residuum for ϵi = 2 and 4 eV (bottom two panels). The top part shows results for both right-hand sides μ = e, o of the linear system; the bottom two parts display just μ = e.

Convergence of the COCG method for the model of Estrada et al.;33 see also the caption of Figure 5 for details.

Let us first focus on the graphs at the top of Figure 5 showing the performance of the GMRES method. The method converges rather well (less than 700 iterations) even without any preconditioning. The convergence is extremely fast below 2 eV (several dozens of iterations) but gets slower above this energy with maximum around 4 eV. This is related to the spectrum of the anion. The electron with energy below 2 eV does not have enough energy to populate vibrational states of the anionic potential. The process of the electron scattering is therefore almost elastic, which means that the wave function is not much perturbed with respect to the initial state used for starting the iterations. Above this energy, the dynamics is much richer, which is reflected in the increased number of iterations needed to reach the converged wave function. For the most of the energies in the range of interest, the preconditioning reduces the number of iterations considerably. The least efficient preconditioning matrices M^dg^ and M^gu^ (overlapping curves in Figure 5) include only a diagonal portion of matrix A and are therefore numerically very cheap to implement. The preconditioner M^du^ includes also terms proportional to coupling constants κ_1_ and κ_2_. For the ECD86 model, there are no terms in matrix A added by increasing the size of the preconditioner to M^d^ and M^u^, and the convergence curves thus overlap for these three preconditioners. The best results are obtained with the preconditioning matrix M^g^, which has blocks of size 2Nu × 2Nu and includes terms proportional to λ in eq 18. The convergence of the residuum norm in the lower part of Figure 5 shows a difference in the behavior of different preconditioned methods. While the best method with the M^g^ preconditioner converges exponentionally for all energies, there is a kind of plateau in the other methods, and the iterations without preconditioning can even surpass some preconditioned methods.

The behavior of the COCG method (Figure 6) is different in several aspects. The overall number of iterations is approximately three times larger (for unpreconditioned iterations), but we have to keep in mind that the COCG method is much simpler with computational demands constant over the course of the iterations. For GMRES, the computational demands for one iteration grow quadratically with the number of iterations. The COCG method does not have the minimization property (40). This is reflected in the shape of the convergence curves (two bottom graphs in Figure 6). Unlike in similar curves for GMRES, here the residuum can locally grow, although in general it finally converges to zero. The efficiency of the different preconditioning schemes is similar like in GMRES, although for higher energies only the M^g^ preconditioner is useful.

Performance of the Methods for the Generalized Model

The generalized model has a more complicated structure (27) of the level-shift operator F(E), which is reflected in a more complicated structure of matrix A; see Figure A1. Surprisingly, the iteration methods converge more quickly with this matrix. There are no clear criteria rigorously relating the structure of the matrix to the speed of convergence. We believe that the faster convergence here may be related to the fact that operator F in the generalized model increases diagonal elements of matrix A. Apart from a little bit faster convergence, the graphs in Figure 7 for the GMRES method in the new model look qualitatively similar to those for the ECD86 model. The norm of the residuum is monotonously decreasing for all methods, and the preconditioner M^g^ is again the most efficient. The individual preconditioners now lead to different convergence rates because all choices of the diagonal blocks are distinct for the richer structure of A. The exception is the equivalence of M^du^ and M^u^ preconditioning. This can be nicely understood from the structure of matrix A depicted in Figure A1. We see that the large and small black diagonal boxes in the bottom right matrix differ by a blank area of zero matrix elements.

Convergence of the GMRES method for the generalized model; see also the caption of Figure 5 for details.

The faster convergence for the new model is even more apparent for the COCG method in Figure 8. Now, all preconditioning schemes except for M^dg^ and M^gu^ are faster than direct iterations.

Convergence of the COCG method for the generalized model; see also the caption of Figure 5 for details.

To conclude the numerical experiment section, we add a few notes on the implementation. Even without utilizing the structure of matrix A, we have got by 1 order of magnitude faster calculation of the spectra utilizing the Krylov-subspace iteration methods as compared to a direct solver. Optimizing the matrix-vector multiplication using the structure of matrix A explained at the beginning of the Appendix leads to another order of magnitude speed up. From the previous examples, we see that the proper choice of preconditioning leads to the decrease of the number of iterations needed for convergence by another 1 order of magnitude for both models and both methods.

Performance for the Model of e– + CO2

In the final part of this section, we discuss our earlier work^34−36^ on the electron collisions with the carbon dioxide (CO_2_) molecule in the context of the present paper. The vibronic coupling model for the e + CO_2_ system^35^ follows the general approach presented here in Section “Generalized Model with Vibronic Coupling with Continuum States”; however, the model is more complex. We considered the nuclear motion within the full four-dimensional vibrational space in combination with three electronic states (^2^Σ_g^+^ virtual state and two components of the ^2^Πu_ shape resonance), which are coupled upon bending of the molecule. The Hamiltonian is thus a 3 × 3 matrix in the electronic space and we did not restrict its elements only to the first order in the normal coordinates (some of the elements were expanded up to the fourth order). Additionally, the three discrete states were coupled to four electron partial waves. The vibrational dynamics is described analogously to the scheme given in Section “Numerical Representation of the Dynamics” but there are four vibrational indices instead of two. The vibrational basis was constructed from products of eigenfunctions of 1D harmonic oscillators for symmetric and stretching modes and eigenfunctions of the 2D harmonic oscillator expressed in polar coordinates for the two-dimensional bending mode.

Using the COCG method without any preconditioning, the number of iterations needed to reach the convergence with the stopping criterion of 10^–3^ (sufficient to obtain converged cross sections) rapidly grows with the electron energy; see Figure 9. For energies above 3 eV, even 2 × 10^5^ iterations were insufficient to reach the convergence; therefore, a suitable preconditioning is essential.

Number of iterations needed to solve the Schrödinger equation for the e + CO2 system using the COCG method without and with preconditioning. In the latter case, the curve is multiplied by a factor of 10.

The slow rate of convergence or no convergence at all is caused by the coupling of the discrete states through the bending mode. The stretching modes do not affect the convergence much since we found that the COCG method converges badly even for the case where we did not consider the stretching modes.2 Thus, taking a block-diagonal preconditioner where blocks contain discrete states and two-dimensional bending was a natural choice. Such a preconditioner is analogous to the preconditioner M^g^ that performs the best for the ECD86 model and its generalization. In the case of CO_2_, around 200 iterations were sufficient to reach the convergence for an initial electron energy of 3 eV; see Figure 9.

Discussion of Resulting Spectra for Test Models

It is not the purpose of this paper to study in detail the calculated spectra and their interpretation. This will require a detailed analysis of the final-state distribution and shape of the individual components of the wave function in the coordinate representation and its relation to the shape of potentials and also study of the dependence of the results on the model parameters. It is a quite voluminous work that deserves a separate paper. We also identify specific molecules that can be treated with the model of the current setup or a proper generalization. We already published the generalization of the model^35^ needed to describe the resulting spectra for the CO_2_ molecule^34^ and performed the detailed analysis^36^ including the final-state distribution, the wave functions, and decomposition of spectra due to the contribution of components of different symmetry.

In the following, we show and briefly describe the 2D spectra for the ECD86 model (which were not the subject of the original paper) and for our new generalization of the model. We also separate the contribution of the two right-hand sides in eq 30 corresponding to the gerade and ungerade symmetry, and finally we study the energy dependence of the cross section for excitation of the fundamental modes.

2D Spectrum for the ECD86 Model

The calculated 2D spectrum for the ECD86 model is shown in Figure 10. The intensity given by eq 17 is plotted as a function of both energy loss Δϵ and initial electron energy ϵ_i_ in a color logarithmic scale. It is fully converged result, i.e. it is independent of the method used to calculate it. Interestingly enough, the spectrum is qualitatively quite similar to the 2D spectrum for the CO_2_ molecule^24,34^ in the region of the Π_u_ resonance. The bulk of the spectrum is located at energies of the incident electron between 2 and 4 eV. This is a consequence of the shape of the anion potential manifold (Figures 1 and 2) and its location relative to the potential of the neutral molecule. The understanding of the detailed shape is not trivial. For small electron energy losses, the spectrum is discretized by vibrational frequencies whose ratio is approximately 3:1. But since this ratio is not exact, the spectrum becomes quasi-continuous for energies above 1 eV. At the same time, we see that there is some selection mechanism that singles out narrower structures close to the diagonal threshold line. There are also diagonal rays appearing in the structure of the spectrum (better apparent in the decomposition of the spectrum according to symmetries). Both of these features were present in the case of CO_2_, where we performed the detailed analysis.^36^

2D electron energy-loss spectrum for the ECD86 model (top) and its decomposition to the gerade and ungerade symmetries (bottom). The intensity of spectrum is shown in a logarithmic scale.

2D Spectrum for a New Model

We proposed the new model above to consistently introduce the vibronic coupling in the level-shift operator in the ECD86 model and to test the iteration schemes to solve the dynamics in this model. The choice of the model parameters was guided by our experience with the diatomic molecules but apart from that, the choice is completely random. To our surprise, the resulting spectrum (see Figure 11) has a quite interesting intricate structure, which is furthermore similar to experimental data for some molecules, like benzene and its derivatives.^52^ Particularly, we are talking about the wedge-shaped structure marked in Figure 11. The origin of this structure is not clear and since it is quite common in experimental data, we will dedicate the future study to this phenomenon. It indicates some selection mechanism in the dynamics that forces the system to skip through a region with small energy losses to large losses.^52^

Same data as in Figure 10 but for the new model. The black dotted lines in the top figure indicate a triangle containing the feature discussed in the text. This feature extends between the horizontal and diagonal lines almost all the way toward zero energy loss.

Vibrational Excitation Integral Cross Sections

To better understand the role of the conical intersection in the vibrational excitation process, we show the individual cross sections σ_vf ← vi__ for the elastic channel vf_ = (0, 0) and for the excitation of one quantum of the two modes, vf = (1, 0) and (0, 1) as the function of the initial electron energy ϵ_i_ in Figure 12. Further, we tested the effect of the off-diagonal couplings F12 and F21 of the discrete states through the continuum on the dynamics. The dashed curves in Figure 12 are the results of the calculation where these two terms were switched off by setting the parameters a1^o^ and a2^e^ (see Table 2) to zero. First of all, we observe that the nonlocal coupling has a little effect on the elastic cross section except for the energies in the vicinity of the position of the conical intersection (around 2.5 eV), where the presence of this additional coupling terms increases the cross section by ∼30%.

Integral cross sections for the elastic channel and the excitation of the two vibrational modes. The dashed curves are the data for a model with the off-diagonal coupling though the continuum switched off.

The data nicely reflect the fact that the discussed coupling terms are proportional to the coordinate qu of the nontotally symmetric mode and thus have a larger influence on the excitation of the ungerade mode vf = (0,1). The cross section from the fully coupled model for this state peaks around ϵ_i_ = 3.0 eV, while the dashed curve exhibits a triple maximum close to ϵ_i_ = 2.3 eV, which can be understood by examining the width of the adiabatic states. The width of the upper branch of the conical intersection gets wider as the coordinate qu increases (see the left panel of Figure 3), while the width nearly vanishes in the lower branch close to the intersection. This enhances the cross section of the (0,1) excitation above the intersection energy of 2.5 eV and suppresses it below this point. The width on the lower branch gets broader again close to the crossing with the potential energy surface of the neutral state, which is reflected by the fact that the cross section for the full calculation again exceeds the dashed curve at small energies. Also note that the potential energy curves of the limited model (not shown) are very similar to those of the original ECD86 model (see Figures 1 and 2). The peaks close to ϵ_i_ = 2.3 eV in the dashed blue curve can thus be traced to the minima in the potential energy surface of the anion similar to those in the left panel of Figure 1.

In contrast to this behavior, the cross section for the excitation of the symmetric mode (1,0) is not very sensitive to the presence of the additional terms, and both the solid and dashed curves follow the same behavior with similar magnitudes.

We expect that the off-diagonal coupling will be even more important for the behavior of the differential cross sections because the involvement of different partial-wave components of the electron continuum in the matrix of the coupling amplitudes Vϵ is subjected to selection rules, as discussed above. Each partial wave contributes to the angular dependence of the differential cross section in a different way as we discussed in detail for the CO_2_ molecule.^34−36^ In this series of papers, we found that the detailed analysis of the spectra and differential cross sections requires also the investigation of the wave functions and the variation of various parameters of the model. This goes beyond the scope of this paper, which focuses on the derivation and numerical solution of the model. However, it will be the subject of a future study that will focus on the interpretation of 2D EELS spectroscopy in terms of the potential energy curves of the discrete states and their widths and couplings to different continuum channels.

Conclusions

We derived a generalization of the model of conical intersection in electronic continuum proposed originally by Estrada et al.^33^ by including terms linear in the vibrational coordinates also in the term that couples the two discrete states of the original model to two partial waves of the electronic continuum. The generalization thus produces linear and quadratic terms in the nonlocal level-shift operator F(E) that describes the dynamics of the vibrational excitation of the molecule by collision with an electron. The dependence on the nontotally symmetric vibrational coordinate also allows indirect coupling of the two discrete states of different symmetries through the electron continuum.

To solve the vibronic dynamics, we implemented two Krylov-subspace iteration methods, GMRES and COCG, which are ideally suited for this kind of model that produces a sparse-matrix representation of the Schrödinger equation. Essential for the efficiency of these methods is the choice of a suitable preconditioner. We tested both methods in detail, and we provided guidelines for the choice of the preconditioner.

We also briefly analyzed the resulting two-dimensional electron energy-loss spectra and the energy dependence of the integral cross sections. These spectra exhibit a surprisingly complex structure, and their full understanding will require a more detailed study with the inspection of the wave functions and perhaps also time-dependent calculations. Here, we show that the separation of the contribution of different irreducible representations to the dynamics disentangles the complex pattern to some degree. Another insight was gained by inspection of the energy dependence of cross sections for some final states. The elastic cross section and the excitation of the symmetric mode are not very sensitive to the presence of the indirect coupling through the continuum except in the vicinity of the energy of the conical intersection. On the other hand, the excitation of the nontotally symmetric vibrational mode is strongly influenced by this coupling. This aspect may be important for many molecules with a nontrivial symmetry group.

We believe that the methods tested here can be used for more complicated molecules to get a better understanding of the 2D electron energy-loss spectroscopy. We plan a more extensive parameter study to obtain a deeper understanding of the results. The proposed method is conceptually simple and can be further generalized in a straightforward way to include more anion states, more vibrational degrees of freedom, and higher-order polynomial functions. A more challenging generalization will be needed to also include dissociative channels and anharmonicity of the neutral molecule.