Quantum dynamics of vibrational energy flow in oscillator chains driven by anharmonic interactions

A new model of vibrational energy transfer in molecular systems taking into account anharmonic (third order) interactions of localized vibrations with a chain of harmonic oscillators is developed. The role of the energy spectrum of the chain and of the magnitude of the non-linear coupling is discussed in detail by an exact numerical solution of the quantum dynamical problem based on the tensor-train (matrix product state) representation of the vibrational wave function. Results show that the type of wave packet motion is determined by the eigen-spectrum of the chain and by its excitation time. It is found that when the excitation of the chain takes place on a much shorter timescale than the energy transfer along the chain the vibrational wave packet moves in a ballistic way independently of the length of the chain. On the other hand when the excitation of the chain takes place on the timescale of the energy transfer along the chain the overall motion becomes superballistic. These findings shed new light on recent observations of ballistic energy transfer along polymethylene chains.


Introduction
Vibrational energy redistribution in molecular systems is one of the fundamental topics of modern physical chemistry. Bonds are broken and formed as a result of energy transfer between nuclear degrees of freedom moving on a certain potential energy landscape: it is this flow that ultimately controls the outcome of all chemical reactions. The main features of the energy flow depend on the topography of the potential energy surface (PES) and on the initial conditions of the system, yet the determination of general laws is impossible if not for the simplest cases.
A wide range of theoretical studies investigating models consisting primarily of chains of bonded atoms demonstrated how a variety of phenomena can emerge even in such simple systems. Most of these studies are based on second order (local) potentials and nearest neighbour couplings [1][2][3][4][5][6][7][8]. However, when second order models of vibrational energy flow are employed to compute relaxation and dephasing times a very poor agreement with available experimental data is obtained [9].
On the one hand, the role of the anharmonicity of the PES-mainly Fermi resonances-on the vibrational energy relaxation and dephasing has been recognized since the early 70's [9,10], while the first approximate rate theory of vibrational energy flow in molecules induced by cubic couplings was developed by Logan and Wolynes in the 90's [11][12][13]. Rate theory based approaches have been widely applied during the last decades [14][15][16][17], while, due to the inherent mathematical and numerical complexity, full quantum dynamical treatments of vibrational energy flow processes in anharmonic systems are still restricted to a few cases [18][19][20][21][22].Furthermore, investigations on model systems have been focussed almost exclusively on the properties of the chain which carries the vibrational energy [2], and little attention has been paid to the way the initial vibrational excitation is created, and to the extent of the delocalization of the wave function during the energy flow. On the other hand, ultrafast IR and 2D-IR measurements have shed new light on energy relaxation paths in molecular chains [23,24,[24][25][26][27][28][29]. In these experiments vibrational energy redistribution is observed after an excess energy has been deposited in a spatially localized region of a system by photo-excitation. The region of the molecule being excited is associated with a localized optically active vibration such as C = O, C ≡ N, amide or C-H, stretching. After the initial excitation the energy content of other localized vibrations can be monitored by some type of time-resolved spectroscopic technique [26]. A prototype of these systems, recently studied by Rutbtsov and co-workers, is sketched in figure 1. Normal mode analysis reveals that the N ≡ N and the C = O stretching vibrations are localized on a very small group of atoms. Upon excitation of the N ≡ N stretching, the energy is transferred to the C = O mode via a two step mechanism: first the excess energy is transferred to the chain by a subpicosecond process, and then it travels 'ballistically' reaching the acceptor carbonyl group after a time that is proportional to the length of the polymethylenic chain. Theoretical studies based on approximate rate theory have highlighted the importance of the vibrational band structure of the polymethylenic chain on the overall energy transfer process suggesting that the ballistic energy transfer is obtained because the energy is distributed among high frequency vibrations with little or no friction at all [28].
However, rate theories cannot provide information about the structure of the vibrational wave function during the process, and cannot disentangle the interplay between coherence and delocalization in the efficiency of the energy transfer. While it might be possible to tackle the problem by a direct computation of the multi-dimensional PES of the system depicted in figure 1, here we follow a different route which aims at disentangling the role of anharmonic and harmonic terms of the PES in the transport mechanism.
The remainder of the paper is organized as follows: in section 2 we develop a realistic microscopic model of vibrational energy transfer in oscillator chains by taking into account the initial process of photo-excitation of localized vibrations and the subsequent transfer of energy to a chain via cubic anharmonic interactions; the tensor-train (TT) methodology is briefly introduced in section 3; then in section 4 an accurate numerical study of the model is performed using recent developments of the time-dependent TT formalism [30][31][32][33][34][35][36]. The main conclusions are summarized in section 5.

Model of energy flow in a chain
The physico-mathematical formalism required to describe the process of vibrational energy transfer in molecules is rooted in the well known GF method by Wilson [37], and in the extensions that have been developed to include anharmonic interactions (see for example [38]). Once a proper set of internal coordinates, comprising bond distances, bond angles and dihedrals, has been chosen the prescriptions to write the kinetic energy operator are readily implemented, while the computation of the PES remains a formidable task but for the simplest cases.
Here we wish to construct a model system that mimics a molecular structure which is characterized by two localized vibrational modes, Q A and Q F , interacting with the ends of a chain of harmonic oscillators. The vibration Q A interacts with a radiation field via a dipole operator enabling a direct photo-excitation of one side of the chain (see figure 2).
The most generic structure of the Hamiltonian operator H implementing the above model can be written in the form ( = 1) The localized vibrations A and F are represented as harmonic oscillators with frequencies ω A,F , respectively; H c is the Hamiltonian of the chain connecting A and F and V(Q A , Q F , {q i }) represents the interaction of the two modes Q A and Q F with the chain modes {q i }; finally the term μ A Q A E(t) describes the dipole interaction of the oscillator A with the external electric field (the dipole moment is assumed to be linear on the oscillator coordinate Q A , E(t) is the dimensionless electric field envelope, and μ A is the coupling where the actual values of the elements of the G and F matrices define the eigen-spectrum of the chain. We refer the readers to reference [37] for a detailed derivation of the G and F matrices and for a comprehensive explanation of the GF methodology. Here q i represent local possibly internal coordinates, and p i are their conjugate momenta. The overall number of oscillators in the model, including the Q A and Q F modes, is N = N c + 2. Rather than trying to compute H c for a specific molecular structure we model the chain by defining a central frequency ω 0 and the bandwidth of its spectrum. Therefore H c can be written in the final form Here D is a coupling constant which determines the bandwidth of the energy eigenvalues of the chain normal modes. While it might seem an oversimplification, direct computations of the spectrum of polymethylenic chains indeed show this type of band structure for several types of vibrations like C-C stretchings and CH 2 waggings [39,40], which, as we shall see in the next sections, represent a most relevant route of energy transport. The mechanism of energy flow from A to F depends both on the chain properties, and on the couplings between the localized vibrations and the chain. We assume that these latter are Fermi type interactions, i.e. cubic anharmonicities, which involve the vibrations Q A , Q F and the first and last oscillator of the chain, q 1 and q N c , respectively This type of interaction can induce transitions in which one quantum is destroyed in the localized vibrations and two quanta are created in the chain (and viceversa). In order for the cubic interactions to effectively drive energy transfer to the chain, the frequencies ω A,F of the two vibrations Q A , Q F must be close to the double of the characteristic chain frequency ω 0 (ω A,F ≈ 2ω 0 ). Here we neglect any other cubic coupling between the modes Q A , Q F and the remaining modes of the chain. This choice is motivated by the observation that the third order coupling coefficients are determined mostly by the geometrical overlap of the normal modes [41], therefore it is equivalent to state that the harmonic mode Q A (Q F ) has a significant overlap only with the harmonic local mode q 1 (q N c ). The relevance of third order coupling terms is also discussed in [18,19] although within a quite different model. We further mention that a similar type of energy transfer mechanism has also been described very recently in a Mo 2 -N 2 complex and has been suggested as a possible way to activate specific chemical bonds [42,43]. Finally, we remark that while most experimental data on IVR are obtained from condensed phase experiments, the model described above does not include any dissipative process. This is justified by the different timescales of inter-and intra-molecular vibrational energy transfer, which are also experimentally observed in the systems described above [26]. In order to clarify the relation of the above model with other approaches in which IVR processes are triggered by anharmonic interactions it is useful to transform the Hamiltonian of the chain into its diagonal form by a unitary transformation (see appendix A) where the frequencies of the eigen-modes, Q k , satisfy the dispersion relation ω k =ω 0 +2D cos(πk/(N c + 1)). Therefore the chain represents a collection of harmonic oscillators with frequencies centred around ω 0 and having a bandwidth 4D. The eigen-modes of the chain are related to the site modes by the linear transformation Q k = l T kl q l , where the explicit structure of the T matrix is given in the appendix A. The eigen-mode transformation also induces a transformation of the anharmonic interaction The explicit expression for the coupling constants β Xmn is given in the appendix A and depends both on the parameters β A,F and on the eigenvector matrix T. Therefore, the initial model with the Hamiltonian of equation 1 is completely equivalent to a system in which the harmonic oscillators Q A , Q F are coupled, via cubic interactions, to a collection of harmonic oscillators with eigen-energies falling in the range This model can be easily extended to include local anharmonicities of the chain but for the moment we skip this complication. After having defined the equivalence of the two models it is clear that the energy transport does not depend exclusively on the chain properties but rather on the interplay between the chain structure, i.e. its spectrum and eigenvectors, and the third order coupling constants β A , β F . In the next sections we address the problem of the role played by these parameters in the mechanism of energy flow by tackling directly the quantum dynamical problem using numerical techniques based on TTs.

Numerical methodology
The quantum dynamical analysis of the above IVR model requires efficient numerical methods for the solution of the associated time-dependent Schrödinger equation. Here we employ a special type of representation of wave functions known as TT format (or matrix product states, MPS, in the physics literature) which has turned out to be a promising approximation in multi-dimensional problems [34,[44][45][46][47][48][49]. Below we sketch the basic principles of the TT decomposition, and show how it can be applied to efficiently solve the time-dependent Schrödinger equation. The reader is referred to the original papers [34,46,49] for a detailed analysis of the TT decomposition.
Let us consider a generic expression of a state of a d dimensional quantum system in the form where |i k labels the basis states of the kth dynamical variable, and the elements C(i 1 , . . . , i d ) are complex numbers labelled by d indices. If we truncate the summation of each index i k the elements C(i 1 , . . . , i d ) represent a tensor of rank d. The evaluation of the summation in equation (8) requires the computation (and storage) of n d terms, where n is the average size of the one-dimensional basis set, which becomes prohibitive for large d. Using the TT format, the tensor C is approximated as where G k (i k ) is a r k−1 × r k complex matrix. In the explicit index notation The matrices G k are three dimensional arrays, called cores of the TT decomposition. The ranks r k are called compression ranks. Using the TT decomposition of equation (9) it is possible, at least in principle, to overcome most of the difficulties caused by the dimensions of the problem. Indeed, the wave function is entirely defined by d arrays of dimensions r k−1 × n k × r k thus the required storage dimension is of the order dnr 2 . The accuracy of the approximation is controlled by the value of the compression ranks, with larger r k providing better approximations.
In a time-dependent theory the cores G k (i k ) are time dependent complex matrices whose equations of motion can be found by applying the time-dependent variational principle (TDVP) to the parametrized form of the wave function The resulting equations of motion can be written in the form and provide an approximate solution of the original equation on the manifold of TT tensors of fixed rank, M TT . In equation (12),P T (G(t)) is the orthogonal projection operator into the tangent space of M TT at |Ψ(G(t)) . Several techniques exist to compute the time evolution of TT/MPS [35,48,50,51]. Here we adopt the so-called adaptive 2-site time-dependent variational integrator (2TDVP) [30,35]. With respect to the one-site TDVP integrator this scheme has the advantage to dynamically increase the ranks r k of the cores of the TT decomposition until convergence is reached or up to a maximum threshold value. We refer the reader to the recent review [30] for numerical details and properties of the 2TDVP integrator. All the computations presented in this paper have been performed using a in-house code based on the itensor software library [52]. We also mention that previous exact treatments of IVR in anharmonic chains reported in reference [19] employ the multi-layer multi configurational time-dependent Hartree method.

Computational results
Taking as starting point the system depicted in figure 1, in the following numerical study we let the frequency of the modes A, F be 2000 cm −1 , which is comparable to the usual frequency of a carbonyl stretching vibration and assume ω 0 = 1000 cm −1 [28]. The quantum dynamics of the model system is analysed as a function of the third-order coupling parameters β A,F , of the bandwidth parameter D and of the length of the chain. In order to keep the system at a size that is physically realistic but still suitable for a numerical quantum dynamical treatment we employ 16, 32 and 64 site chains. Although the two Hamiltonian operators of equations (1) and (7) are equivalent we prefer to adopt the local mode representation [equation (1)] for the actual numerical solution of the problem since it is more suitable for TT methodologies. Furthermore, the local mode representation enables a direct visualization of the energy flow along the chain in real time (see below).
Since only high-frequency vibrational modes are considered, thermal effects are negligible and temperature can safely be set to zero. Hence the initial state of the system is obtained first by determining the ground state of the full anharmonic Hamiltonian and then exciting it with an ultrashort laser pulse. In order to simplify the numerical procedures we assume a δ(t) envelope of the electric field. We point out that the ground state of the system |G already contains small correlations between the vibrational modes Q A , Q F and the chain modes due to the third order couplings. However, upon the instantaneous excitation only the mode Q A is allowed to absorb energy, thus energy is deposited only at one end of the chain. This way of preparing the initial state is very close to the real physical process involved in ultrafast time-resolved IR measurements.
Since we use the TT/MPS representation of the wave function the ground state |G of the Hamiltonian operator in equation (1) is obtained by using the numerical density matrix renormalization group theory (DMRG) [53]. Once we have obtained |G , the dipole operator Q A is applied to simulate the instantaneous photo-excitation process |ψ(0) is not an eigenstate of H and its evolution in time is obtained by solving the associated time-dependent Schrödinger equation. In all our numerical simulations the threshold in the singular-value decompositions of the tensors is set to 10 −10 , ensuring energy conservation during the evolution which is fundamental to correctly describe long-time hydrodynamic behaviour of quantum systems.
In figure 3(a) the populations of the chain sites and of the oscillators A and F are reported as a function of time for the case of D = 50 cm −1 and β A,F = 150 cm −1 . The resulting bandwidth of the harmonic oscillator chain is 4D = 200 cm −1 which is typical for the C-C stretching and CH 2 wagging bands of polymethylene chains. The first thing to notice is that the Fermi resonance effectively removes the single quantum that has been created by photo-excitation from mode Q A enabling a very fast energy transfer to the chain which takes place in about 200 fs. This initial depopulation of the mode Q A is Gaussian, n A (t) = exp{−(t/τ A ) 2 }, where τ A = 71 fs is the so-called Zeno time which can be evaluated as described in [54,55]. After 800 fs the excitation of the last vibration of the chain reaches its maximum, and then the occupation number of the mode Q F starts to increase reaching a maximum at 1000 fs. Interestingly, the maximum value of the occupation numbers of the chain sites decreases along the chain but for the last two sites. We notice that the occupation number n 1 , of the chain oscillator q 1 which is directly coupled to Q A reaches the value of 1.3 at t ≈ 100 fs, which reflects the Fermi-resonant nature of the coupling. In fact, the complete annihilation of one quantum on mode Q A and the creation of two quanta on the first site of the chain is never observed due to quantum interference effects caused by the remaining sites of the chain. After 200 fs, when n A (t) ≈ 0, there are at least 6 vibrations which share the excess energy. The extent of delocalization of the wave packet can be easily determined by looking at the inverse participation ratio defined as [33,56,57] IPR(t) = 1 when there is a maximum of localization in a single site of system, while IPR(t) ≈ N when the energy in fully delocalized on the entire system (remind that N = N c + 2 is the total number of oscillators in the system). In figure 3(b) IPR(t) is reported as a function of time. Since the initial wave packet is fully localized on site A, IPR starts from its minimum value of 1 and then increases with time reaching a value of 8 around 800 fs. This means that not all sites of the chain have the same excitation level. From 800 to 1200 fs the IPR decreases. Comparing the IPR value with the occupation numbers of the chain of figure 3(a), this sudden decrease of the IPR value can be clearly associated with a localization of the excess energy on the last two sites of the chain and on the mode Q F . After this transitory localization, the IPR increases again reaching a maximum value of 10 after 1.5 ps. This corresponds to a significant delocalization of the energy between the degrees of freedom of the system.
The computed occupation numbers clearly indicate that only a small fraction of excess energy becomes actually localized on mode Q F .
A global description of the dynamical evolution of the vibrational wave packet is provided by its second moment The square root of this quantity follows a general power law of the type √ m 2 (t) ∝ t ν where ν = 0 corresponds to complete localization, ν = 1/2 to a classical small-step diffusion, ν = 1 to ballistic motion, ν ∈ (0, 1) to anomalous diffusion, and ν > 1 to superballistic motion.
In a complex dynamical process the parameter ν will not be a constant, and the type of dynamics can be described as piecewise ballistic or diffusive according to the approximate value of ν in a specific time  interval. From figure 3 we can observe that up to 800 fs ν ≈ 1, suggesting that the motion is almost perfectly ballistic. After 800 fs the curve flattens and fluctuates corresponding first to a partial localization of the energy at one end of the chain and then to a diffusive-like behaviour. Figure 4 shows the results of dynamics for a system with 34 harmonic oscillators, i.e. with a chain of 32 sites, and a Fermi coupling constant β A = β F = 150 cm −1 . The linear dependence of √ m 2 (t) up to 1800 fs suggests a ballistic type of motion which is followed by a localization up to about 2.5 ps. The time for the observation of the excitation of the Q F mode is about 2000 fs, exactly the double of the time observed in the 16 oscillator chain model. This is clearly the result of the ballistic motion of the system which follows the excitation of Q A . The qualitative behaviour of the IPR is also similar to that described for the 16 site model, with an initial increase up to 1800 fs followed by a decrease as a result of the energy localization at the end of the chain and then to an increase associated with a spread of the wave packet along the chain.
As can be clearly seen from figure 5, doubling again the length of the chain to 64 sites gives a qualitatively similar picture. Indeed, the energy transfer time, τ ET , defined as the time required to observe the maximum occupation number on the last site of the chain doubles, see figures 4 and 5. In the models treated so far the ballistic behaviour is predominant causing a linear increase of τ ET with increasing the length of the chain. This behaviour is in line with experimental findings obtained by time-resolved 2D-IR spectroscopy on molecules belonging to the family described in figure 1 [25,28].
The vibrational dynamics described so far is characterized by two distinct timescales: the ultrafast de-excitation of the mode Q A with the time τ A which corresponds to an ultrafast energy transfer to the harmonic chain, followed by ballistic motion along the chain, characterized by the time τ ET (the characteristic energy transfer time) that is proportional to the length of the chain.
However, if the coupling between the excited oscillator Q A and the chain is lowered to β A = 50 cm −1 , then τ A = 213 fs, the separation of timescales starts to fade at least for the shortest 16 site chain, and vibrational wave packet dynamics becomes more complex. This behaviour can be clearly observed in figure 6(a) where the three bell-shaped curves represent the occupation number, n F (t) , of mode Q F for N c = 16, 32, 64 oscillators, and the decaying dashed line represents the population of mode Q A . Indeed, with increasing N c from 16 to 32 the transition time does not double, while changing from 32 to 64 insures that τ ET τ A , and a clear doubling of τ ET is observed. Figure 7 shows a direct comparison of evolution of the occupation numbers and of the root mean square deviation, √ m 2 (t), of the 16 site chain for two different values of the Fermi resonance coupling parameter, β A,F , while keeping D = 100 cm −1 . When the coupling parameter is fairly large, β A,F = 150 cm −1 , the mode Q A transfers completely the excess energy which travels along the chain to reach the mode Q F in a ballistic type of motion ( √ m 2 (t) ∝ t). This two stage process is associated with a quite large population of the harmonic oscillators of the chain (see dotted curves in figure 7(a)). We also notice that around 1.3 ps a revival of the occupation number of mode Q A (black full line) is evident: energy is reflected back from Q F to Q A . During this localization process √ m 2 (t) reaches a minimum. Conversely, when the coupling parameter is small, β A,F = 50 cm −1 , the chain keeps a very small level of excitation during the entire process. This can be clearly seen from the occupation numbers of the chain oscillators which remain much smaller than in the previous case. In this sense the chain is not an efficient transport layer. When this type of quantum transition becomes predominant the energy transfer proceeds in a superballistic way. Indeed, examination of figure 7(b)) reveals that √ m 2 (t) ∝ t ν with ν > 1 at least up to 600 fs.
In order to identify the main features of the present model that rule the transition from ballistic to superballistic motion it is convenient to analyse the root mean square deviation √ m 2 (t) for different values of the bandwidth 4D and couplings β A,F as reported in figures 8-10. The system with a 16 oscillator chain, see figure 8, shows two distinct behaviours. For D = 50 cm −1 the motion is ballistic up to almost 1 ps, with a minor deviation emerging, as already discussed, for the vale β A,F = 50 cm −1 . With doubling D the motion looses its ballistic behaviour but for the largest value of β A,F = 150 cm −1 . For small values of the cubic coupling constants a higher order component appears in the behaviour of √ m 2 (t) and the dynamics can be described as superballistic (ν > 1). A similar analysis holds true for the model with a 32 oscillator chain  (see figure 9). The main difference with the previous case is that the ballistic behaviour here is more pronounced even for small values of β A,F . On the other hand, the 64 site chain exhibits a ballistic behaviour in almost all conditions, with small superballistic features observable at short times and only for the smallest value of β A,F .
An alternative point of view on the dynamics of the energy flow can be derived by an examination of the transformed Hamiltonian operator of equation (7). For very large Fermi resonance couplings β X , X = A, F the transformed parameters β Xkl remain large enough to couple effectively the modes Q A,F to a significant fraction of the chain eigen-modes Q k . The large density of coupled states makes the energy flow from Q A to the chain very close to an irreversible rate process. As a rule of thumb this happens when ρ(E)β Akl 1, where ρ(E) ≈ N c /4D is the average density of states of the chain, and the brackets mean some average value. In the opposite limit, i.e. ρ(E)β Akl 1, resonant tunnelling through a few energy levels of the chain becomes important (see right panel of figure 2), causing deviation from the ballistic regime.
We further notice that in standard IVR theories Fermi resonances are known to be effective in the energy transport mechanism when the ratio between the coupling constant and the resonance offset becomes large, i.e. β Akl /(ω A − ω k − ω l ) 1, whereas in the current model their efficacy is mostly  determined by the bandwidth parameter of the chain, D. However, the two criteria are not in contrast with each other and the former is recovered when the chain Hamiltonian is transformed into the eigen-mode representation.

Conclusions
In this paper we have developed a new model for energy transfer in molecular systems. The model comprises two local oscillators connected via cubic anharmonic couplings to the two ends of a chain of harmonic oscillators. We have numerically studied the evolution of a vibrational wave packet after the initial excitation of one local oscillator using TT (matrix product state) techniques combined with a 2TDVP integrator.
The results show a remarkable variety of wave packet dynamics depending on the specific choice of a few system parameters, namely the magnitude of the cubic coupling constants, the energy bandwidth of the chain of harmonic oscillators, and the length of the chain. It has been found that two regimes exist for the energy transport process. The first is a two step mechanism in which the energy is first transferred completely from the local oscillator Q A to the chain, and then travels ballistically along the chain to reach the local oscillator Q F . This mechanism is active when the cubic coupling is fairly large and the bandwidth of the chain is relatively small. The second mechanism appears when the time required for the mode Q A to lose its excess energy is comparable to the time required for the energy to be transported along the chain itself. Under these conditions a non-linear dependence of the square root of the mean standard deviation m 2 (t), which is associated with a non ballistic spread of the vibrational wave packet, is observed.
The first regime is very likely active in several molecular systems described recently in the scientific literature [26]. The second regime might not be stable under ordinary conditions, as it requires a quasi-resonant tunnelling and very small cubic coupling constants. Although this condition is very likely fulfiled, small anharmonicities of the PES of the chain, which have been neglected in the present treatment, might strongly affect the wave packet dynamics. Work is currently in progress to assess this point.
which, upon expansion gives Using equation (A.4) we obtain The maximum value of the transformed third-order parameters is β Akl = 2β A /(N c + 1) for k = l = (N c + 1)/2. This result implies that the values of the couplings β Akl depend on the length of the chain, and, in our case, they are one to two orders of magnitudes smaller that the actual value of β A,F .