Compacting the density matrix in quantum dynamics: Singular value decomposition of the surprisal and the dominant constraints for anharmonic systems
We introduce a practical method for compacting the time evolution of the quantum state of a closed physical system. The density matrix is specified as a function of a few timeindependent observables where their coefficients are time-dependent. The key mathematical step is the vectorization of the surprisal, the logarithm of the density matrix, at each time point of interest. The time span used depends on the required spectral
<!>I. INTRODUCTION<!>( )<!>II. COMPACTING THE TIME EVOLUTION OF A QUANTUM STATE<!>A. The surprisal<!>B. Vectorization<!>C. Compacting the surprisal<!>D. Summary<!>III. DYNAMICS ON THE MORSE POTENTIAL<!>A. The algebra for a Morse potential<!>B. The coherent state as an initial state<!>C. Analytical expression for the time-dependent surprisal<!>D. Dominant constraints from SVD<!>IV. NON-ADIABATIC DYNAMICS IN PYRAZINE<!>V. CONCLUDING REMARKS<!>SUPPLEMENTAL MATERIAL
<p>resolution. The entire course of the system evolution is represented as a matrix where each column is the vectorized surprisal at the given time point. Using singular value decomposition, SVD, of this matrix we generate realistic approximations for the timeindependent observables and their respective time dependent coefficients. This allows a simplification of the algebraic procedure for determining the dominant constraints (the time-independent observables) in the sense of the maximal entropy approach. A nonstationary coherent initial state of a Morse oscillator is used to introduce the approach. We derive analytical exact expression for the surprisal as a function of time and this offers a benchmark for comparison with the accurate but approximate SVD results. We discuss two examples of a Morse potential of different anharmonicities, the H 2 and I 2 molecules. We further demonstrate the approach for a two coupled electronic states problem, the well studied non radiative decay of pyrazine from its bright state. Five constraints are found to be enough to capture the ultrafast electronic population exchange and to recover the dynamics of the wave packet in both electronic states.</p><!><p>Physics tells us how to propagate the state of the system in time. For a closed system these differential equations are fully agreed on. In classical mechanics it is Newton's law or Hamilton's equations. In quantum mechanics it is the time-dependent Schrödinger equation. What physics does not tell is how to specify the initial state of the system which is a necessary input to the differential equations of motion. Of course, the initial state is ours to specify. But the question is, in general, how to tell what it is. There are persuasive and power arguments why to use an initial state of maximal entropy. [1][2][3][4][5][6] In general the initial state is not as 'uniform' as possible but it is subject to constraints that it must satisfy.</p><p>In quantum mechanics these constraints are taken to be operators, Âk , k = 0,1,.. whose mean values, Âk = Tr Âk ρ ( ) restrict the choice of the unknown density matrix ρ . Often one takes the zeroth constraint to be the identity, Â0 = I ! , to insure the normalization of the unknown density matrix ρ . Subject to reasonable mathematical limitations on the operators that constrain ρ 3, 6 one shows that the density of maximal entropy has the exponential form:</p><p>where the λ k ' s are the Lagrange multipliers that are used in the standard procedure of There are as many equations as there are unknowns but even in the classical limit solving the implicit equations is not fully trivial. 7 In quantum mechanics the problem is far more demanding because the operators need not commute so the exponential form needs to be re-expressed as a product form. 8 Our purpose in this paper is to demonstrate a practical route to the computations of the Lagrange multipliers in the quantum mechanical case of a density matrix that is evolving in time and we know it at different points in time. We here use a direct time propagation 9 of the surprisal, the logarithm of the density matrix, which for the state of maximal entropy is a linear function of the constraints:</p><p>. This enables us to simplify determining the values of the Lagrange multipliers.</p><p>One can seek to determine a density matrix of maximal entropy that is explicitly a continuous function of time. Jaynes 10 and Zubarev 11 advocated this approach with special reference to macroscopic open systems. In particular Zubarev emphasized the quantum mechanical case and proceeded to derive kinetic equations and other relaxation processes and to relate it to non-equilibrium thermodynamics. Our aim is different. We seek to compact the density matrix for a closed physical system, a density that is propagated according to the Liouville-von Neumann equations of motion. In principle we can approach or even reach an exact description but first we want to identify a set of dominant constraints, constraints that capture the essence of the dynamics. Time-independent constraints that are valid throughout the time evolution and that are derived from the dynamics.</p><p>For a system that is closed one can show that if the initial state is a state of maximal entropy it will remain a state of maximal entropy throughout its time evolution. [12][13][14] Even more, if the initial state is a state of maximal entropy and the system is closed it is immediate to write down a set of time-dependent constraints that exactly specify the system. 12 To determine explicit, time-independent, constraints it is very advantageous if one can generate a set of operators that remains finite under the action of the Hamiltonian of the system. Explicitly, the commutator of any constraint with the Hamiltonian Âk , Ĥ ⎡ ⎣ ⎤ ⎦ must then be expressed as a linear combination of the other known constraints. Examples include systems of coupled harmonic oscillators 15 including for several electronic states. 9,16 For harmonic systems the needed constraints are derived from the well-known algebra of creation/annihilation operators. 17 There are four operators per oscillator â , â † , â † â ,I ! or six operators â , â † , â † â , â2 , â †2 ,I ! when the frequency of the oscillator can change. 12 When the set of constraints does not remain finite under the operation of commutation with the Hamiltonian we make do with approximate constraints. We will put these in a hierarchy starting with those with the largest Lagrange multipliers. These we call the dominant constraints. Nature helps in that it is often the case that there is a hierarchy of scales of the Lagrange multipliers. In such cases the algebraic route to identify the dominant constraints becomes harder to implement as one needs to guess/identify which constraints will dominate during the dynamics. In this paper we aim to develop a straightforward numerical route to the dominant constraints for any given surprisal matrix.</p><p>Strictly speaking our approach does not require that the density matrix is an exponential function as in equation (1). All that we strictly need is that the density matrix, our initial state, can be approximated as a function of the constraints. As an example, in his review Fano 18 expands the density matrix as a linear function of observables. In such a case one can apply the methodology presented in this work directly to the density matrix.</p><p>It is however compelling to argue that amongst all possible density matrices that are functions of the constraints one should chose the one whose entropy is maximal.</p><p>Realistic oscillators in molecules are typically anharmonic so in this paper we explore an algebra for a dynamics of a Morse-like oscillator. That one can use the SU(2) algebra for the energy levels of the Morse oscillator is well studied. [19][20][21][22][23][24] Here however we determine the algebra for a non-stationary coherent state 25 and derive an analytical expression for the time-dependence of the surprisal. This expression is then used to compare with the results of the here proposed numerical procedure for the determination of the Lagrange multipliers. In addition this form of the surprisal offers an analytical basis for computing recurrences and revivals in the anharmonic dynamics. [26][27][28] The density matrix of a closed system is a Hermitian square matrix. At any given time point one can diagonalize the density matrix. One can diagonalize a density matrix of maximal entropy even when the individual operators do not commute. 6,29 Then one can retain only a small set of leading terms in the spectral resolution of this matrix. The eigenvalues of a normalized density matrix will be in descending order below the unity.</p><p>Retain those eigenvectors with the largest eigenvalues meaning those eigenvalues nearest to unity. Convergence is achieved when the sum of the eigenvalues is acceptably near one.</p><p>The density matrix is exactly reproduced when the sum of the eigenvalues is exactly one.</p><p>We can repeat this for the next time point of interest but this approach does not relate the results at different time points, the eigenvectors of the density will depend on time. It is not an analogue of the maximal entropy approach where each Lagrange multiplier is determined as a function of time and the different constraints are the same for different time points. Hence we need to proceed in a different way, a way when all the time points of interest are treated together. The time span needs to be such that it suffices to provide the spectral resolution required. We seek to express the quantum time evolution of a closed system for all times of interest by one matrix, a matrix that can then be compacted. The matrix need not be square so to 'diagonalize' it we will use the procedure of Singular</p><p>Value Decomposition, SVD.</p><p>For some time we have been using SVD to compact big biological data. [30][31][32][33] Also for biology we meant to represent the data in a manner suggested by the maximal entropy approach identifying on the classical level constraints dominant for a specific process of interest. At a given time point data such as gene expressions are represented as a vector. In quantum dynamics we use a finite basis representation for the operators, 9 hence the density, surprisal or constraints are matrices. To get the density matrix and/or its surprisal as a vector we use the operation of vectorization, 34 Scheme 1. This retains the complete data of the matrix in a vector format. A set of surprisal matrices represented at a series of time points becomes thereby a single 'unified' matrix where each column of the unified matrix is a matrix of the surprisal at a particular point in time.</p><p>Scheme 1. Vectorization of the surprisal matrices, ( )</p><p>, computed at L time points, n t , in an N by N finite basis representation. The columns of each matrix are stacked sequentially in one vector, and then all the vectors are concatenated in one 2 N by L rectangular matrix T .</p><p>Remarkably, the route to determining the Lagrange multipliers as a function of time is rather similar in biology and in the quantum mechanical case. This is a technical point that is discussed below and it is a direct result of a mathematical property of vectorization. If</p><!><p>vec A is the vectorization of matrix A then ( ) ( ) † † ( )</p><p>where the trace is the scalar product of two operators in the Hilbert Schmidt sense. Thereby the constraints identified via SVD route in a vectorized form as orthogonal right eigenvectors,</p><p>= δ kl , are also orthogonal in the Hilbert Schmidt scalar product of their respective operators,</p><p>In coupled systems the number of constraints can be quite large but they are typically far from being all equally important. The advantage of SVD is that it effectively identifies the 'dominant' constraints meaning those constraints whose weight, (their Lagrange multiplier), is high. Already in our very first application of SVD 30 the importance of constraints varied in time. Some constraints can be important throughout. Others can contribute more at earlier time, etc. In practice SVD is applied to data covering a finite time span, see scheme 1. The range in time determines the available spectral resolution which in an anharmonic system can vary. Therefore the constraints determined for a small time span in SVD can be different from those determined for longer times up to the largest period in the system.</p><p>We discuss two examples of application of the proposed methodology. In the first example we contrast the non-equilibrium vibrational dynamics in H 2 and I 2 molecules.</p><p>Morse potentials of these two molecules have different parameters of anharmonicity, different number of bound states. We compare the dynamics given by the analytical expression with the results obtained by considering only the dominant constraints provided in the SVD procedure. We derive an analytical form for the surprisal. It shows that the algebra is not closed under commutation with the Hamiltonian. For an exact representation of the surprisal there are as many constraints as the size of the vibrational basis. Still we show by example that for early time dynamics the number of dominant constraints can be much smaller than this maximum.</p><p>In the second example we describe the non-equilibrium dynamics both in electronic and vibrational degrees of freedom. We use a well-studied system in nonadiabatic dynamics the radiationless decay of the bright state of pyrazine. [35][36][37][38] Here the degree of compaction is higher, out of 800 constraints in finite basis representation we identify 5 constraints that describe the population dynamics both in electronic and nuclear Hilbert spaces. The detailed description of coherences between different electronic states can need more constraints.</p><p>In section II we begin with the general description of the SVD approach applied to the surprisal of a quantum state. Section III provides analytical and numerical results for the surprisal in the case of the dynamics on the Morse potential. Section IV presents the numerical results for the compaction of the non-adiabatic dynamics in pyrazine.</p><p>Concluding remarks are discussed in Section V.</p><!><p>The compaction proceeds in three steps. Determine the surprisal, vectorize the surprisal, apply singular value decomposition to a matrix whose columns are the surprisals at different points in time. To make the mathematics simpler we take all matrices to be of finite, if perhaps large, dimension N. In the actual numerical work all operators are to be understood as N by N matrices in the same space. The dimension needs to be large enough not only for accurately describing the initial state but also for the later time dynamics.</p><p>Strictly speaking this can be an approximation but it is quite often numerically accurate enough.</p><!><p>The surprisal Î is the logarithm of the density matrix, Î = − ln ρ . It is the quantum mechanical extension of the classical surprisal. 12,39 The expectation value of the quantum expansions where the operators are also orthogonal). The time-evolution of the surprisal is determined by the time-evolution of these projectors, time-dependent constraints. In terms of a complementary time-independent constraints we aim to write the surprisal, cf.</p><p>equation ( 1), as:</p><p>Here the time dependence is carried by the Lagrange multipliers. The constraints are Hermitian time independent operators and we typically arrange them to be mutually orthogonal in the sense of Hilbert-Schmidt, 9 Tr Âk † Âl ( ) = δ kl . This orthogonality will be shown to be inherently present in the results of the approximation by SVD. It is a very useful property because it gives rise to a practical expression for the Lagrange multipliers:</p><p>Sometimes it is convenient to work with a not Hermitian constraint. Then its Hermitian adjoint, for example â and â † , needs to be included in the sum. The Hermiticity of the density or of the surprisal means that the Lagrange multipliers of two such conjugate constraints must be complex conjugate to one another.</p><p>In the numerical propagation of the surprisal 9 we compute ˆ( ) I t in the finite basis representation via Liouville-von Neumann equation of motion:</p><p>Rather than propagating the surprisal in time using equation ( 4) we can, where convenient, use a basis that diagonalizes the Hamiltonian. This allows propagating the matrix elements of the surprisal in parallel where each one is given by:</p><p>Here i ε is an eigenvalue of the Hamiltonian Ĥ . Note that depending on the initial state the surprisal is not necessarily diagonal in this basis.</p><!><p>At each time point t the surprisal matrix is vectorized which we write as a column</p><p>It is often useful to write the time propagation of the surprisal as:</p><p>where</p><p>) is the time-evolution operator from the initial time to the present. There is a corresponding result for the vectorized surprisal:</p><p>The 2 N by 2 N matrix ( ) t Λ propagates the surprisal vector in time.</p><!><p>We use singular value decomposition 40 as a mathematical tool to provide a compact representation of the surprisal that is valid over a range of points in time. There is a subtle mathematical point about the meaning of convergence that we will point out. However when the compacted surprisal provides a good fit the reservation is not of practical importance.</p><p>We consider a sequence of L time points { } 0 1 2 1 , , ,.., L t t t t − and define a rectangular matrix T of dimensions 2</p><p>N by L. Each column of this matrix is a vectorized surprisal at the time point with the column index. Since the number of states N is typically quite large L is the smaller dimension of the rectangular matrix T . The range of L depends on the required resolution. For high spectral resolution of the dynamics one needs to span a wide range. This can generate a set of dominant constraints that are valid throughout. Otherwise the constraints can be limited to describing a short time span where low frequency details may not be resolved which may be important in anharmonic systems.</p><p>To determine the singular values of T we construct the L by L symmetric non-negative matrix T † T and diagonalize it:</p><p>The matrix T † T is a real positive matrix. In the basis of eigenstates of the Hamiltonian, Eq. ( 5), its diagonal elements are constant. This and other key properties of the matrix are discussed in section S1 of the Supplemental Material. The non-negative eigenvalues are written as ω α 2 and we arrange them in descending order. Some of the eigenvalues can be zero or zero to within the numerical noise. The normalized real eigenvectors L α have L components and for different eigenvalues α are mutually orthogonal. If there is a degeneracy the eigenvectors can be arranged to be orthogonal.</p><p>The matrix TT † is 2 N by 2 N dimensional. It is however not necessary to diagonalize this typically enormous matrix because as shown in general 40 or in particular by equation (10) the non-zero eigenvalues of TT † are the same as the L eigenvalues of the matrix T † T . We do however need the eigenvectors of TT † that correspond to these L eigenvalues. A practical route 30 is to multiply equation ( 9) from the left by T :</p><p>TT TL TL (10)   The complex valued L vectors ( )</p><p>To each vector G α of 2 N components one can relate a square N by N matrix G α . From the definition of vectorization one can show that the orthogonality of the different eigenvectors G α of TT † can be written as a Hilbert-Schmidt orthogonalization of the different matrices G α :</p><p>The compact singular value decomposition is a numerically exact spectral-like representation of T using the right and left eigenvectors G α and L α : 40</p><p>At the particular time t and for the component n of the vectorized surprisal this yields:</p><p>The surprisal is compacted by using fewer terms α in the expansion equation ( 14). To make a clear connection with maximal entropy we write, using ( ( ))</p><p>nt n vec t = I T meaning the n'th component of the vector vec(I(t)) and defining the real valued ( )</p><p>M in equation ( 15) is the number of constraints that are used to approximate the surprisal.</p><p>By our convention on the labeling of the eigenvalues, ω 0 > ω 1 ≥ ω 2 .. The λ α (t)'s are also in decreasing order since the eigenvectors L α are normalized. Technically 'the best' is meant here in the sense of a norm of the difference between nt T and the left hand side of equation (15). This is almost but not quite the solution of maximal entropy. At maximal entropy we aim to minimize the difference but where each point is weighted by its term in the density ρ n (t) . In other words, maximal entropy aims to do a better job for the larger elements ρ n (t) .</p><p>The orthonormality of the constraints enables us to write that for each value of α : † ( )</p><p>Recall that the index n is a pair of indices so that the equation reads in more practical detail as in equation ( 6):</p><!><p>Consider the quantum mechanical surprisal as a matrix in a finite, N dimensional, space. To describe the dynamics the matrix is given for each time point in a sequence of L points. For each time point the N by N surprisal matrix is vectorized. The sequence of matrices at time points is thereby expressed as a sequence of vectors of typically high, 2 N , dimensions. These vectors are arranged as columns of a matrix T . We diagonalize the L by L dimensional matrix T † T and arrange the L eigenvalues in decreasing order. These eigenvalues and respective eigenvectors determine the Lagrange multipliers,</p><p>, in order of decreasing importance. The large dimensional matrix TT † has the same non-zero eigenvalues as T † T and its high dimensional eigenvectors with nonzero eigenvalues can be determined without an actual diagonalization. These are the time- (10). The matrix elements of the compacted surprisal with n ≡ k,l and for M constraints are given by:</p><p>The density matrix is determined from the surprisal as noted in section II.A. The new feature of quantum surprisal analysis is that the surprisal is a not necessarily diagonal matrix. In SVD the matrix elements of the surprisal are given in equation (18). To have an explicit expression for the density matrix it is therefore necessary to exponentiate the surprisal matrix, ρ = exp(−I) . This can be efficiently done by diagonalizing the surprisal matrix and using the spectral representation. Explicitly:</p><p>where s labels the eigenstates of the surprisal matrix and ε s are the corresponding eigenvalues. It is often the case that the surprisal matrix is not diagonal in the basis of physical interest and then the populations, the diagonal elements ρ kk , are determined also by the off diagonal elements of the surprisal. In other words, the populations are controlled by the coherences. In the basis of the eigenstates of the Hamiltonian the diagonal part of the surprisal is steady in time. The corresponding constraint thereby defines a density that is a steady state.</p><!><p>When there is a closed algebra that can specify the initial state and the elements of this algebra are closed under commutation with the Hamiltonian there is an exact expansion of the surprisal using a finite number of time independent constraints. When the algebra is not closed under commutation with the Hamiltonian we can only expand the surprisal in time dependent constraints. In practice this often means that approximations are needed. To benchmark the proposed approach for this general case and to illustrate the formalism we here use a wave packet moving on a Morse potential, a coherent superposition of the eigenstates of the Morse potential. Therefore the populations of the different vibrational states will be constant in time and only the coherences will vary. The more states there are in the wave packet, the quadratically more coherences need to be described. The higher the anharmonicity, the more the periods of the coherences will differ. In the coordinate representation the density will vary and it will do so to reflect the coherences between the vibrational states. The analytical results for the time dependence of the surprisal clearly show the periodicity and time revival of the wave packet.</p><p>Surprisal analysis using algebras based on harmonic oscillators have been extensively discussed. 9,12,14,41,42 Here we wanted to emphasize that neither the algebraic approach nor surprisal analysis are limited to coupled but harmonic systems. To analyze the effects of the anharmonicity on the dominant constraints and their Lagrange multipliers we contrast the results of the dynamics that unfold on the Morse potentials of H 2 and I 2 molecules.</p><!><p>Using an algebraic approach for the bound states and also the continuous spectrum of a</p><p>Morse potential has been extensively reported. [19][20][21][22][23] Here we use the SU(2) algebra and notation as introduced in Ref. 22. The intention was to make it as similar as possible to the harmonic oscillator case where we use â, â † ,I ! . There is a creation Q+ and an annihilation operator Q− that are adjoints on one another Q± ( ) † = Q∓ . The replacement for the identity operator is 2 Q0 , so that the commutation relation</p><p>tends to an identity only in the harmonic limit. The Hamiltonian can be written as</p><p>) in a form analogous to that of the harmonic limit. A is the energy scale and equals the anharmonicity, A = ω e χ e . 23 We define the phase such that Q0 , Q±</p><p>with the result that the analogue of the identity commutes with the Hamiltonian Q0 , Ĥ ⎡ ⎣ ⎤ ⎦ = 0 . We derive the analytical expressions for the surprisal and its time change in an algebraic form. To compute we need the matrix elements of the operators in the (finite) basis of bound states of the Morse oscillator. These are given in the sections S2-S3 of the Supplementary Material (SM).</p><!><p>To generate motion we need to start from a not stationary initial state. A convenient choice is a coherent state. There is not a unique definition of a coherent state for an anharmonic motion. 43 We chose one that is analogous to what is used for the harmonic oscillator namely a shifted ground state. 17 It is the definition used by Gilmore et al 25 but note that our Hamiltonian, as expected for an anharmonic case, 20 is bilinear in the operators of the algebra. This leads to a rather complex motion. 26,27,44 Even so we are able to derive an analytical expression for the surprisal as a function of time. This exact expression is used as a benchmark for comparison with a surprisal analysis as described in section II.</p><p>The coherent state, the density matrix at time zero, is defined by a displacement of the vibrational ground state 0 of the Morse potential:</p><p>where D is the displacement operator:</p><p>The complex parameter ζ is the analogue of the complex parameter α in the definition of the coherent state of a harmonic oscillator. 17 One can rewrite the displacement operator as a product of exponentials, see section S2 of the Supplemental Material (SM) but this is not essential for us because we will work with the surprisal which is the logarithm of the density. To take the advantage of the algebra we write the density 0 0 as the very low temperature limit, of the thermal density matrix Z −1 exp(−β Ĥ ) . Here Z = exp(λ 0 ) is the partition function and βω > 1 where ω is the harmonic frequency, so β A = βω / (2k +1) is finite (k is the total number of bound states). Thereby the surprisal at the initial time is defined as:</p><p>The key details of the transformations needed to rewrite the surprisal are spelled out in the SM, section S2. The resulting expression is given by:</p><p>H.C. stands for the Hermitian adjoint, 2 Q0 + 1</p><p>) . Therefore, the initial set of constraints consists of a set of seven operators:</p><p>} that are not closed under commutation.</p><!><p>The operator 0 Q is diagonal in the basis of the Morse vibrational states but 0 2 Q tends to a multiple of the identity operator only in the harmonic limit, 2 Q0 → (2k + 1)I ! , where</p><p>is the total number of bound states. 23 Thereby the operators Q+ and Q− are almost but not quite closed under commutation with the Hamiltonian. For example, under the action of the evolution operator Û we have:</p><p>Note how the constraint ( )</p><p>breaks the linearity in the dependence of the surprisal on the time-independent constraints. This will show up in the SVD procedure: the number of dominant time-independent constraints will increase in time. The rate of growth in the number of constraints is proportional to the anharmonicity parameter e e</p><p>A ω χ = and not dependent on the extent of the initial shift ζ .</p><p>The final expression for the time propagated surprisal is:</p><p>where the operator K = ln Z + β Ĥ cos 2 2ζ + β Asin The operators Q+ and Q+ 2 and their Hermitian conjugates are off diagonal in the basis of Morse eigenstates, equation (S6) of the SM. Their matrices are one below the diagonal and one twice below and vice versa. Therefore the route to the elements of the density matrix is not immediate, see equation (19).</p><p>Next we express the surprisal, Eq. ( 25), in the form where it depends linearly on timeindependent operators. We use the finite basis representation, the basis of eigenstates of Morse potential, , k ν , where k is the total number of bound states of Morse potential and ν denotes the vibrational quantum number. In this representation the operator 0 Q is diagonal:</p><p>, and hence also its function:</p><p>The identity A ω e =1 (2k + 1) enables taking the harmonic limit. Inserting this in Eq. (25)   gives the following:  27) is an exact representation of the time-dependent surprisal via five time-independent constraints but with Lagrange multipliers that are state dependent. In the harmonic limit of a large number of bound states, as in equation ( 26), the Lagrange multipliers do become a function of time only.</p><!><p>The exact dynamics of the surprisal, equations ( 25)-( 27 The hydrogen molecule, H 2 , has a high frequency, 4395 cm -1 , and small number, 18, of bound states. The frequency of the iodine molecule, I 2 , is much lower, 215 cm -1 , and the number of bound states is 174. I 2 is therefore almost exactly an order of magnitude more harmonic than H 2 . More details on the computational parameters are given in section S3 of the SM. We define the two initial states using the same parameter ζ for both molecules, see Fig. First we examine the accuracy that is achievable by keeping a finite number of terms in an SVD expansion for the surprisal, Eq. (18). By diagonalizing the surprisal we obtain the density, equation (19). We plot the density of the wave packet as a function of the bond length as shown in Fig. 2. The plot is for the initial states with a larger, 2ζ shift so that more excited vibrational states are accessed, see Fig. 1. The movies S1-S4 available in the supplemental materials show the propagation of the wave packet in time and the approximation by M = 9, 14 and 19 time-independent constraints for smaller and larger shift for each molecule. For all the times of the wave packet rattling between its two turning regions the approximation given by 9 constraints for the smaller shift and by 14 constraints for 2ζ is accurate. Adding five more constraints does not improve the fit to graph reading accuracy showing the convergence of the SVD expansion, Eq. ( 18). The resulting mean values and dispersion of the coordinate are shown in Fig. S3 of the SM. In the vibrational basis of a Morse potential the populations of the states should be constant in time and only the coherences will vary. The corresponding numerical results using SVD are shown in Fig. 3 for different numbers of constraints computed from short time dynamics. The population of vibrational states in hydrogen, a molecule with not many bound states, is well represented with nine constraints for the small shift and 14 for the larger shift, Fig. 3. For the larger shift using only 9 constraints provides not constant populations and the norm deviates from unity. This is because the off diagonal elements of the surprisal increase with the shift ζ and the coherences in the surprisal control the populations as discussed following equation (19). In iodine, where at a higher shift many more states are in a superposition, more frequencies contribute. Hence one needs to increase the number of constraints up to 14, see Fig. S4 in the SM. The coherences are also shown in Fig. 3. The diagonal elements of the surprisal, equation ( 26) are analytically seen to be independent of time. The time dependence of the approximate populations, the diagonal elements of the density matrix comes from the off diagonal matrix elements of the surprisal contributing to the diagonal elements of the density matrix. The magnitude of the analytical off diagonal matrix elements of the surprisal increases with the shift.</p><p>The results shown in Figs. 2 and 3 are computed as follows. The Lagrange multipliers are linearly proportional to the time-independent singular eigenvalues ω α of the T matrix, ( )</p><p>λ ω = L see equations ( 14) and (15). To extract the dominant constraints, we arrange the eigenvalues in the order of decreasing size as seen in Fig. 4. The fast drop in magnitude is clear for all four cases shown. This enables us to keep only the first few terms in the SVD expansion, thereby providing the best approximation for the exact dynamics with that number of terms. Here 'best' means best in the sense of L 2 norm.</p><p>The leading eigenvalue ω 0 is almost the same for the different values of the shift parameter ζ and in H 2 it is an order of magnitude bigger than the next eigenvalue. The shift does significantly affect the eigenvalues for the other constraints. These other constraints reflect the off diagonal elements of the surprisal and the off-diagonal elements are more significant for the larger shift. For both values of the shift Fig. 4 shows that the leading constraints for α > 0 are about degenerate pairs. It is also the case that their conjugate Lagrange multipliers have a phase shift, see Fig. The eigenvalues α ω do not depend on time. The dynamics of the Lagrange multipliers is determined by the real eigenvectors L α of the covariance in time matrix</p><p>, see Fig. 5. In terms of the Lagrange multipliers the time dependence of the surprisal matrix is given as:</p><p>The matrix elements of the surprisal in the basis of vibrational eigenstates of the Morse potential are given in the analytical derivation in equation ( 27). This equation shows that the diagonal elements do not depend on time while off diagonal do. To provide an accurate approximation the SVD procedure should yield the Lagrange multiplier λ 0 to be time independent and the conjugate matrix G 0 that is diagonal. When only a short time span</p><p>Δ is used for SVD it is seen in Fig. 6 that λ 0 is constant but only to better than 1% accuracy. Comparison of the corresponding G 0 diagonal elements with the analytical results is shown on Fig. 6. The changes of the diagonal elements with respect to the vibrational quantum number are accurately captured by SVD and only a small deviation of the SVD results from the analytical terms is observed. However, when a long time span is used, a time span sufficient to include the longest periodicity in the data, the resulting λ 0 (t) is time independent. From equation (27) it follows that the off diagonal elements of the surprisal should be strictly confined to two sub and two super diagonals. This is fully recovered by the SVD expansion, all the constraints with 0 α > have non-zero off-diagonal values on the first two rows above and below the diagonal, see Fig. In conclusion, over a wide dynamic range we find that the SVD procedure provides a practical and accurate description of anharmonic motion on a single potential. In the limit of a long time span, characterized by the lowest frequency in the system, the SVD procedure provides a Gelfand basis representation. We turn next to a well studied multielectronic-state problem.</p><!><p>The non-radiative transfer between the bright and dark states of pyrazine Hamiltonian is defined following Domcke, 36 by taking into consideration two vibrational normal modes for two excited electronic states, bright * π π − (S 2 ), and dark * n π − (S 1 ). The anharmonicity of the vibrations is due to the coupling mode, with frequency 10a ω (Fig. 8a), that represents the out-of-plane motion of hydrogen atoms, which breaks the symmetry and is responsible for the coupling and population exchange between the bright and dark states. [35][36][37][38] The other mode is a tuning mode, with frequency 1 ω (Fig. 8b), that is inducing the nuclear motion towards the conical intersection between the two electronic states. The Hamiltonian is defined in a diabatic picture as: The initial state is localized in the Franck Condon region on the bright electronic state.</p><p>In this Franck Condon region both vibrational modes are undisplaced. Very rapidly dynamics unfolds bringing the wave packet to the conical intersection, Fig. 8(c). Exact computational results on a 2D grid of nuclear coordinates for the population and coherent dynamics are plotted in Fig. 9. To generate the surprisal matrix we use 20 vibrational states in each vibrational normal mode in a diabatic picture of the electronic states as in equation (29). The time duration we use is 25 fs as the transfer is rather rapid. See the details in the section S3 of the supplemental material. For comparing with the exact numerical dynamics shown in Fig. 9 are the SVD results using 5, 9 and 26 constraints. Even 5 constraints capture the population transfer, Fig. 9(a) and the overall motion of the wave packet, Fig.  In the basis of zero order states the constraints have matrix elements that carry three</p><p>indices. An electronic state index n =1,2, a vibrational index, s =1,2,..,20, for the coupling mode and an index, q =1,2,..,20, for the tuning mode. To present the results we need to agree on the order of these indices. We therefore arrange the matrices in blocks. The elements of these matrices are labeled as:</p><!><p>To represent a quantum mechanical time evolution we need to account not only for the populations of the states but also for their coherences. In an anharmonic system the coherences exhibit a range of periods. The demands from an approximation are therefore tighter. The surprisal offers a tool that has a double parentage. On the one hand it can be motivated through the information theory based maximum entropy formalism. This expresses the surprisal as a linear combination of time independent operators with time dependent coefficients. The operators are the constraints that restrain the entropy from reaching its global maximum. Their coefficients are the undetermined Lagrange multipliers that prescribe the form of the constrained density matrix. The numerical values of the coefficients are determined so as to enforce the density operator of the system to be consistent with the expectation values of the constraints. Other operators do not constrain the entropy and their Lagrange multipliers have the value zero. In empirical surprisal analysis one determines the numerical values of the coefficients (= the Lagrange multipliers) by a fit to the input data. Such a fit is not fully trivial even in the limit when only the populations are of interest. Here we showed how to perform such a fit in an approximate but accurate manner using singular value decomposition for the quantum dynamics via a vectorization of the surprisal matrix.</p><p>A complementary motivation for the use of surprisal is a Lie algebraic approach to the dynamics. It shows that the surprisal as a linear combination of constraints can be an exact representation of the dynamics also in anharmonic systems and/or when several electronic states are coupled. Indeed, the importance of anharmonicity in surprisal analysis was rather clear from the very earliest applications, see for example, figure 6</p><!><p>Diverse technical aspects and detailed analytical derivations are described in the supplemental material. This includes the discussion of the properties of the timecovariance matrix T † T for a vectorized surprisal of a general form, section S1. The derivation of the time evolution of the surprisal for a Morse coherent state is shown in section S2. Section S3 contains computational details for the SVD procedures and benchmark calculations. Additional figures are listed in the section S4.</p>