Solving vibronic dynamics in electron continuum

Solving Vibronic Dynamics in Electron Continuum Martina Cosicová, Jan Dvorák, and Martin Cízek* Cite This: J. Chem. Theory Comput. 2024, 20, 2696−2710 Read Online ACCESS Metrics & More Article Recommendations 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 CO2by 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. ■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⟩ + +e M v M e M v( ) ( ) i f (1) mediated by one or several metastable anion states M−. After the process, the molecule is left in the final vibrational state |vf⟩. The total energy Eduring the collision is conserved = + = +E E E i v f v i f (2) where ϵare electron energies and Evthe 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− +*+M M e M v( ) ( ) f (3) with initial energy Eof 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 energyloss 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 Received: November 2, 2023 Revised: January 10, 2024 Accepted: January 12, 2024 Published: February 7, 2024 Articlepubs.acs.org/JCTC © 2024 The Authors. Published by American Chemical Society 2696 https://doi.org/10.1021/acs.jctc.3c01217 J. Chem. Theory Comput. 2024, 20, 2696−2710 This article is licensed under CC-BY 4.0 Downloaded via TECHL UNIV OF OSTRAVA on November 18, 2024 at 07:29:54 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles. 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 CO2. 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 projectionoperator 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, adiabaticnuclei 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 = | |d d d (4) as a sum over a set of discrete states |d⟩and the complementary operator =I (5) projecting on the background continuum. The basis in the background part can be chosen as the states that solve the background scattering problem | = + |V, ( ) , el 0 0 0 (6) 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 Φ0is 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 | =d d dd (7) | =d0 (8) | = ( ) (9) The electronic Hamiltonian in the -space is described by a matrix | | = +d d V U el dd dd0 (10) 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 | | =d V q( ) el d (11) These elements represent the vibronic coupling 1 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 el for the electron scattering from the molecule for each fixed position of the nuclei qby Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c01217 J. Chem. Theory Comput. 2024, 20, 2696−2710 2697 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 | = + | = |H v T V v E v( ) N v0 0 (12) where TNis the kinetic-energy operator for the nuclei and vis 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 = + +H H U F ef 0 (13) which is the matrix in the indices d,d′and the operator in the space of vibrational degrees of freedom. In the equation above, H0is the Hamiltonian operator for the vibrations of the molecule multiplied by unity matrix δdd′in the discrete-state indices, Uis the matrix with the elements Udd′defined above, and the operator Fdescribes the dynamical coupling of the discrete-state space to the electronic continuum = [ + ]F E H V q E H i V q ( ) ( ) ( )d dd d d 000 1 (14) 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 vito final state vfand leaving in continuum state |ϵfμf⟩is given by 33 = | [ ] |T v V E H V v v v dd fdef dd di 1 ffiif f i i (15) and is closely related to the integral cross section for the vibrational excitation event = | |T 2 v v i v v 3 2 f i i f ffii (16) Finally, to simulate the full 2D electron energy-loss spectra, we have to collect vibrational excitation cross sections for all accessible final states =S( , ) ( ) ( ) i v v v i v f f i f (17) 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 Δϵvd f=Evd f−Evd i=ϵi−ϵfin 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 = + + U E q q q E q g u u g 1 1 2 2 i k j j j j j j j j y { z z z z z z z z (18) This is a completely general form of the dependence of matrix Uon 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 discretestate-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 -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 dand partial wave μ. We therefore assume that the discrete state d= 1 has even symmetry as the μ=epartial 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) = = | |F V( ) 2Im 2 e 11 11 1 2 (19) = = | |F V( ) 2Im 2 o 22 22 2 2 (20) by means of the integral transform [ ]( )/2 ii defined as (compare eq 14) [ ] = + ff x x i x( ) ( ) d (21) This transform can be worked out analytically for the assumed form of the widths =+ a b( ) exp( ) dd d l d 1/2 d (22) (see Berman et al. 47 ). To complete the model description, we must give the vibrational Hamiltonian H0of the neutral molecule. The model simply assumes harmonic vibrations = + = + + H T V q q q q 1 2 1 2 1 2 1 2 N g g u u gguu 0 0 2 2 2 2 2 2 (23) 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 = + + +E n n 1 2 1 2 g g u u i k j j jy { z z zi k j j jy { z z z The numerical values of ωiand the parameters defining direct coupling matrix Uand discrete-state-continuum matrix Vefor the model studied in Estrada et al. 33 and used here for Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c01217 J. Chem. Theory Comput. 2024, 20, 2696−2710 2698 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 =V q E( ) R i Rloc 2 , obtained from the pole of the fixednuclei 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. 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 dand column index μ =V V q V q V q V q ( ) ( ) ( ) ( ) e o e o 1 1 2 2 i k j j j j j j y { z z z z z z (24) where =V q( ) /2 e 1 11 , =V q( ) /2 o 2 22 , 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 CO2molecule. 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 = [ + ] [ + ] V f q f q f q f q ( ) 1 ( ) ( ) ( ) 1 e g o u e u o g 111 2 2 2 i k j j j j j j j j j y { z z z z z z z z z (25) 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 | | = + f a b2 exp( ) dd l d 2 1/2 d (26) with the values of the parameters given in Table 2. Using eq 14, we see that the structure of the nonlocal levelshift operator F(E−H0) is much richer = + [ ] + + [ ] = [ ] + + [ ] + = + [ ] + [ ] + = [ ] + + + [ ] F q f f q q f f q F q f f q q f f q F q f f q q f f q F q f f q q f f q (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) (1 ) g e e g u o o u u e e u g o o g g e e u u o o g u e e g g o o u 11 1 1 1 11 1 22 2 2 22 2 2 12 1 1 2 1 2 2 21 1 2 1 2 1 2 (27) Table 1. Values of the Parameters Describing the PseudoJahn−Teller Model of Estrada et al. 33 , a parameter value parameter value ωg0.258 ωu0.091 E12.45 E22.85 κ1−0.212 κ20.254 λ0.318 a10.086 a20.186 b10.833 b20.375 l12l21 a The dimension of adis eV1/2−ld(ldis dimensionless), bdare in eV−1, and the remaining parameters are in eV. Figure 1. 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. Figure 2. Perspective view of the potential energy manifold for the ECD86 model. The width (inverse lifetime) is marked by the color scale. Table 2. Values of the Parameters Describing the Generalization of the Model a parameter value parameter value a1 e0.07 a2 e0.1 b1 e0.25 b2 e0.5 l1 e0l2 e0 a1 o0.186 a2 o0.15 b1 o0.375 b2 o0.8 l1 o1l2 o1 λ10.2 λ20.1 a The dimension of ad μis eV1/2−ld μ(ld μis dimensionless), bd μare in eV−1, and the remaining parameters are in eV. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c01217 J. Chem. Theory Comput. 2024, 20, 2696−2710 2699 where we used the integral transform (21) again. The ordering of the terms that depend on qiwith respect to ( ) must be kept because we substitute the operator ϵ=E−H0, which does not commute with the normal coordinates. The offdiagonal 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 4do 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. 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 = | | | 2 v v i vv 3 2 f i i f f f i i (28) where we defined auxiliary wave functions |Φv μ⟩with components = | |n n V v, d n n v g u d, , g u (29) and with energy according to the conservation law (2), i.e., ϵ= E−Ev. The anion wave function |Ψvd i μi⟩satisfies | = |E H( ) ef v v i i i i (30) In the harmonic oscillator basis, this equation represents a system of linear equations for unknown components of the discrete-state wave function = | |d n n, d n n g u v, , g u i i (31) where we introduced a compound index α≡(d,ng,nu). Using this notation, the matrix of this system Aα,α′reads = | |A n n E H n n, ( ) , g u ef dd g u, (32) and the scalar product in eq 28 can be written as the sum over the components of Φα vfμfand Ψα | = ** vv v d n n d n n v d n n , , , , , f f i iff g u g u g u (33) 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= 2NgNucomponent vectors and Ais a N×Nmatrix. 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 blocktridiagonal 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 Ais large but sparse. From the character of the problem, it is also complex symmetric but not Hermitian. The structure of matrix Adepends 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 bwith 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 Figure 3. 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. Figure 4. Perspective view of the potential energy manifold for the new model, colored by the width of the resonance. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c01217 J. Chem. Theory Comput. 2024, 20, 2696−2710 2700 the two iteration methods shows the convergence of the residuum rn=b−Axn→0 for the approximation xnof 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 CO2. 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−1A(MT)−1(see Appendix), where Mis easily invertible approximation to A. In the following, we discuss the preconditioning by matrix Mconsisting of the block-diagonal part of Afor different ordering of basis vectors. We show small blocks Mdg,Mdu, and Mgu of the sizes Nu,Ng, and 2, respectively, and the larger blocks Md,Mg, and Muof 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 CO2molecule. 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,ein eq 30. 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 Figure 5. 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,oof the linear system; the bottom two parts display just μ=e. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c01217 J. Chem. Theory Comput. 2024, 20, 2696−2710 2701 reduces the number of iterations considerably. The least efficient preconditioning matrices Mdg and Mgu (overlapping curves in Figure 5) include only a diagonal portion of matrix A and are therefore numerically very cheap to implement. The preconditioner Mdu includes also terms proportional to coupling constants κ1and κ2. For the ECD86 model, there are no terms in matrix Aadded by increasing the size of the preconditioner to Mdand Mu, and the convergence curves thus overlap for these three preconditioners. The best results are obtained with the preconditioning matrix Mg, which has blocks of size 2Nu×2Nuand 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 Mg 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 Mgpreconditioner 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 Fin 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 Mgis 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 Mdu and Mupreconditioning. This can be nicely understood from the structure of matrix Adepicted 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. The faster convergence for the new model is even more apparent for the COCG method in Figure 8. Now, all preconditioning schemes except for Mdg and Mgu are faster than direct iterations. 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. Figure 6. Convergence of the COCG method for the model of Estrada et al.; 33 see also the caption of Figure 5 for details. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c01217 J. Chem. Theory Comput. 2024, 20, 2696−2710 2702 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 (CO2) molecule in the context of the present paper. The vibronic coupling model for the e+ CO2system 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 fourdimensional vibrational space in combination with three electronic states (2Σg+virtual state and two components of the 2Πushape 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 Figure 7. Convergence of the GMRES method for the generalized model; see also the caption of Figure 5 for details. Figure 8. Convergence of the COCG method for the generalized model; see also the caption of Figure 5 for details. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c01217 J. Chem. Theory Comput. 2024, 20, 2696−2710 2703 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 ×105iterations were insufficient to reach the convergence; therefore, a suitable preconditioning is essential. 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 twodimensional bending was a natural choice. Such a preconditioner is analogous to the preconditioner Mgthat performs the best for the ECD86 model and its generalization. In the case of CO2, 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 CO2 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 ϵiin 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 CO2 molecule 24,34 in the region of the Πuresonance. 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 Figure 9. Number of iterations needed to solve the Schrodinger equation for the e+ CO2system using the COCG method without and with preconditioning. In the latter case, the curve is multiplied by a factor of 10. Figure 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. Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article https://doi.org/10.1021/acs.jctc.3c01217 J. Chem. Theory Comput. 2024, 20, 2696−2710 2704