Yun-an Yan

Guizhou Provincial Key Laboratory of Computational Nano-Material Science,Guizhou Education University,Guizhou 550018,China

Beijing Computational Science Research Center,Beijing 100094,China

Low-Storage Runge-Kutta M ethod for Sim ulating Tim e-Dependent Quantum Dynam ics

Yun-an Yan∗

Guizhou Provincial Key Laboratory of Computational Nano-Material Science,Guizhou Education University,Guizhou 550018,China

Beijing Computational Science Research Center,Beijing 100094,China

A w ide range of quantum system sare time-invariant and the corresponding dynam ics is dictated by linear diff erential equations w ith constant coeffi cients.A lthough sim p le in m athem atical concept,the integration of these equations is usually com p licated in practice for com plex system s,where both the computational timeand thememory storage become lim iting factors.For this reason,low-storage Runge-Kuttamethods become increasingly popular for the tim e integration.This work suggests a series ofs-stagesth-order exp licit Runge-Kutta methods specific for autonom ous linear equations,which only requires two times of thememory storage for the state vector.We also introduce a 13-stage eighth-order scheme for autonomous linear equations,which has optim ized stability region and is reduced to a fi fth-order m ethod for general equations.These m ethods exhibit significant perform ance im provem ents over the p revious general-purpose low-stage schemes.As an exam p le,we apply the integrator to simulate the non-Markovian exciton dynam ics in a 15-site linear chain consisting of perylene-bisim ide derivatives.

Low-storage Runge-Kutta,Autonomous linear diff erential equation,Timedependent dynam ics,Time-invariant Ham iltonian

I.INTRODUCTION

The physics laws of quantum dynam ics,such as the Schr¨odinger’swave equation[1]and the Heisenbergm atrix equation[2]for isolated system s,take the form of linear differential equations w ith constant coeffi cients. In reality,a quantum system is inevitably coup led to an environm ent and thus becom es open.W hen the coupling isweak and thememory time of the environment isshort compared to the intrinsic time scale of the system,a M arkovian app roxim ation ignoring thememory time leads to a good description[3−5].Otherw ise, the hierarchicalapproach is a powerful tool to simulate the non-Markovian dynam ics by embedding the non-M arkovian system in a larger,possibly in finite M arkovian system[6−9].Once again,the equations ofm otion for the open system s take a linear differential form w ith constant coeffi cients ifexternaldriving fieldsarenot app lied.Form ally,the dynam ics for all of these cases are dictated by the ordinary diff erential equation

together w ith an initial condition x(t0)=x0at timet0. Here,x isa state vector of lengthNandMisanN×Nmatrix which does not depend on either the timetor the state vector x.

The numerical integration of such an equation is simp le in concept from themathematical perspective.For system sw ith small tomoderatedimension,one can benefi t from the extensive studies on the linear diff erentialequations[10−14].Theaccumulated results include the Lanczos iteration in K ryslov subspace[15−17],the Chebyshev scheme[18,19],m atrix exponentialmethod [20−22],and the classical Runge-Kuttam ethod[23].

In practice,however,the com putations are com p licated by the big size of the state vector in app lying the quantum princip les for com p lex system s.As the computer memory required to store the state vector grows exponentially w ith the size of the system,the lengthNof the state vector can easily reach 109and even higher in various simulations.For large system s,the quantum dynam ical simulation becomes expansive or even prohibited,which is called the curse of dimensionality[24]. This cursew ith respect to thememory issue isenhanced by the numerical integration because the integrator itself needs additionalm em ory to store the state vector.

For this reason,low-storagemethods draw a lot ofattention.The sp lit operatormethod,which isa sym p lectic integrator not requiring additionalmemory location for the wavefunction,is effi cient for quantum mechanical systems w ithout space-momentum m ixed term s in the Ham iltonian[25].For other system s,low-storageRunge-Kutta(LSRK)methods,which use lessmemory locations than the classical Runge-Kutta schemes,are useful in reducing the memory use[26−29].Actually, the LSRK scheme was proposed at the early time of com puter simulations[30−32],when thememory storage was small.Later,different groupsmade schemes using only two memory locations for any stage number(the number of derivative evaluations in each tim e step of the numerical integration)[27,29,33].Kennedy and co-workers systematically studied the low-storage schem esandm ade the low-storage idea popular[34,35].

Regardless of how powerful the high performance com puting becom es and how big the m em ory storage is available,we always ask more to simulate bigger systems to reveal new mechanism,to seek new phenomena and to design new materials.A partial,alternative remedy for this problem is to develop effi cient algorithm w ith not only shorter com putational time but also less memory demand so that we can simulate bigger systemsw ith available com putational resources.Thus,the low-storagem ethods becam e increasingly popular during the past two decades although the available computermemory was dramatically improved.

The LSRK scheme assumes the sam e accuracy and robustness as the traditional explicit Runge-Kutta methods but only needs two memory locations for the state vector regard less of the order.By contrast,the classicalexplicits-stage Runge-Kuttamethod normally requiress+1 memory locations for the state vector. But,the available LSRK methods were designed for general purpose problem s,which suff ers from the order barrier condition,namely,there is nos-stagesthorder explicit Runge-Kutta method fors>4 and no (s+2)-stagesth-order m ethod fors>7[36].The order barrier is severer for LSRK schem es because there are additional conditions in order to achieve the lowstorage integration.For example,for general nonlinear multidimensional diff erential equation,we have to use five stages to get a fourth-order LSRK method[35]. To have an optimal stability region,a 12-stage fourthorder LSRK scheme isalso w idely used[37].The larger stage number w ill significantly increase the com putational tim e for the simulations of the dynam ics.

Actually,we can develop m ore effi cient LSRK schem es if we take advantage of the linearity of Eq.(1). Zingg and Chisholm[38]suggested a four-stage fourthorder LSRK method when adding an inhomogeneous time-dependent term to Eq.(1).Indeed,more effi cient low-storage algorithm exists for linear diff erentialequationsw ith constant coeffi cients.

In this report,we show thats-stagesth-order LSRK method exists for Eq.(1),so that the time-dependent dynam ics of time-invariant system s can be propagated very effi ciently in the use of both com putational tim e and com putermemory.The developed m ethodsare app lied to the hierarchical simulations of the exciton diffusion in a 15-site linear aggregate.

In this work,we show thats-stages-order LSRK methodsexist forany positivesand suggest theschemes fors=4,6,8,10 and 12.A 13-stagemethod w ith optim ized stability region,which is of eighth-order for linear equation w ith constant coeffi cientsand of fi fth-order for general systems,is also suggested.We also presents a hierarchical simulation of the exciton diffusion in a chain of aJ-type linear aggregate w ith the 12-stage 12th-order LSRK method.

II.LOW-STORAGE RUNGE-KUTTA M ETHOD FOR AUTONOMOUS LINEAR

A.Classical exp licit Runge-Kutta m ethods

The Runge-Kutta methods are derivative-free schem es to find the num erical solution for the initial value problem of the ordinary differential equations

for a given function f.A generals-stage Runge-Kutta method discretizes the time w ith a grid{t0,t1,···, tn,···}and approximates the solution of Eq.(2)stepby-step w ithsiterations in each step,

where∆t=tn+1−tnis the size of the time step,and xnisan approximation to x(tn).Here,the coeffi cientsajk,bjandcjare the elem ents of the Butcher tab leau[23]. To have a com pact p resentation,wew ill treat a asamatrix which only hasnon-zeromatrix elementsajkin the strict lower triangular part.The consistency condition of the Runge-Kutta method im poses the constraints

Therefore,onlyajkandbjare free parameters,which are chosen so that xn+1is thep-order Taylor expansion of x(t)at timetn,

If this condition holds,the Runge-Kutta scheme is called apth-ordermethod.

B.Low-storage Runge-Kutta methods

Eq.(3)clearly shows that thes-stage exp licit Runge-Kutta schemenormally requiress+1memory locations for integration,which isnot favorable for large systems. By putting further constraints on the coeffi cientsajkandbj,we m ay achieve a low-storage method which requires lessmemory locations,even twomemory locations regardless of the stage number.In the literature, there exist diff erent schem es to achieve the two-register m ethod[27,29,33].Here,wew ill be sticked to the formulation originally proposed by W illiamson[29],which takes the follow ing iteration in each time step

Here the register y points to the thememory location storing the state vector x,which contains the initial condition xnbefore the iteration and directly gives xn+1after the iteration.An exp licit schem e im p liesA1=0,so there are 2s−1 independent param eters for ans-stage scheme,based on which the classical Butcher tableau can be constructed[37]

For the autonomous linear equation Eq.(1),all second-order partial derivatives vanish,namely,

As a consequence,the order condition Eq.(8)is greatly simp lified,

Therefore,there arepconstraints in order to have apth-order method,which can normally be satisfied because there are 2p−1 param eters in ap-stage schem e. This im plies thats-stages-order LSRK methods exist for the autonomous linear diff erential equations.A lthough there aremu ltip le solutions of Eq.(14),particular interest is in the choiceAj=−1 forj≥2.W ith such a choice,these algebraic equations can be convenientlyand analytically solved w ith a G r¨obner basis decom position[39].The coeffi cientsBjfor the schemes w ith 4,6,8,10 and 12 stages are given in Table I.Hereafter,themethod w ill be nam ed as LSRKs-p(q),wheresrefers to the number of stages,pdenotes the order of themethod for autonomous linear equations.W hen it app lies,qmeans the order of themethod for general diff erential equations.

TABLE I The coeffi cientsBjof the low-storage Runge-Kuttam ethods fors=4,6,8,10 and 12 for the hom ogeneous linear diff erential equations w ith constant coeffi cients.

In num erical app lications,another im portant issue is the stability of the solver,which requires that small errors in the initial state do not drive the Runge-Kutta approxim ation off from the real value.The stability analysisof the solver is com p licated and case-dependent for general system s.Thus a paradigm for the stability analysis is to use the one-dimensional equation w ith a com p lex coeffi cientλto obtain useful guidance.The coeffi cientλcorresponds to the eigenvalues in multi-dimensional problem s.For this one-dimensional case,the exp licit Runge-Kutta m ethod approxim atesx(t0+n∆t)w ithxn=G(λ∆t)xn−1,whereG(z)is the polynom ial

TABLE II The coeffi cients of the 13-stage eighth-order low-storage Runge-Kuttamethod for autonomous linear diff erential equations.Thism ethod is of the fi fth-order for tim e-variant or non linear diff erential equations.

w ith the coeffi cientsγjbeing

The iteration Eqs.(4−6)is robust against small perturbations if the characteristic polynom ialG(λ∆t)fulfi lls the criterion|G(λ∆t)|≤1.The stability region of the Runge-Kutta method is then defined as the set of all pointszin the complex plane satisfying|G(z)|≤1.

For schemes w iths=p,allγjare determ ined by the order condition Eq.(14)so that there is no freedom to m anipulate the stability region.Therefore,in order to optim ize the stability,we have to sacrifice the order of themethod to have free parameters.

Follow ing Ref.[37],we adopt a two-stage process to find the LSRK method w ith optim ized stability.For a given stage numbersand a desired orderp,we fi rst findγj(p

In princip le,we have to optim ize the stability w ith the exact area w ithin the contour|G(z)|=1.Here,we w ill adopt a much simpler,yet reasonable,approximation to this area,namely,

wherexmaxis the maximum positive real number satisfying|G(−xmax)|=1,andyjis them inimum positive number solving|G(iyj−jxmax/10)|=1,w ithing the imaginary unit.

Onceγjare obtained,we proceed to find the LSRK parametersAjandBj.In thiswork,we are interested in the13-stage 8-ordermethod to havebalance between the order and the stability.W ith this choice,there are five free parametersγj(9≤j≤13)to op tim ize the stability at the fi rst step.The optim ization is carried out by a random search w ith the Nelder-Mead method[40]to find the localm aximum as im p lem ented in the NLopt library[41].The optim ized results forγjare listed in Table II.At the second step,there are 25AjandBjparameters satisfying 13 order conditions.In order to fix the rest 12 param eters,we require the scheme to be a fi fth-order method for general differential equations. Besides the order conditions for linear equations w ith constant coeffi cients,there are 12 additional constraints

The nonlinear Eq.(14)and Eqs.(19−30)are solved w ith the Levenberg-Marquardt algorithm[42,43]as imp lemented in the levmar package[44]w ith amodified long double p recision.M ultip le solutions are found in determ iningAjandBj.The results forAj,Bjandcjw ith smallest sum of the absolute values ofAjare given in Table II,show ing 17 effective digits.

C.Stability contours

FIG.1 presents the stability contours normalized by the inverse of the stage number for all the methods listed in Tables I and II.As the number of stages coincides the order of the methods,the stability contoursof LSRK 4-4,LSRK 6-6,LSRK 8-8,LSRK 10-10 and LSRK 12-12 schemes are uniquely defined and do not depend on how the parametersAjandBjare chosen. For these methods,the stability region is larger when the stage number increases.However,FIG.1 clearly shows that the renorm alized region shrinks for larger stage numbers.W hen the stability is the lim iting factor in the numerical integration,themaximum step size is controlled by the stability region.Then,the effi ciency of the un-optim ized high order m ethods is lower because the increase of the stability region is not suffi cient to compensate the cost in perform ing the calculations for the additional stages.In these cases,it is appreciated to sacrifice the order of the m ethod to increase the stability region.As shown,normalized stability region of the LSRK 13-8(5)scheme is significantly bigger than that of the un-optim ized methods.A dom inant feature is that its contour is elongated along thexaxis.Correspondingly,LSRK 13-8(5)is suitable to simulate the over-dam ped m otion,in which the real part of the eigenvalues is larger than the im aginary part. To have a comparison,we also p lot the stability contour for the LSRK 13-4(4)method suggested by Niegemannet al.[37].Thism ethod have a stability region elongated alongy-axis and is suitable to simulate the under-dam ped motion,in which the imaginary part of the eigenvalues is dom inant.

FIG.1 Stability contours for them ethods suggested in Tables I and II.The contours are symm etric w ith respect to thex-axis,and the lower parts thus are not disp layed.For com parison,also shown is the 13-stage fourth-order schem e [LSRK 13-4(4)]for generaldiff erentialequationsp roposed by Niegem annet al.[37].

D.Accuracy verifi cation

We check the effi ciency of all the methods w ith a symmetric m atrixM.Here,the adopted eigenvaluesλjare

whereN=256 is the dimension of thematrixM.W ith this settings,λ1=0 is the only zero of the eigenvalues, which im p lies that there is a unique stationary state of the system s.

For the sim p licity of computation,we choose the discrete Fourier transform matrixFjk=exp[πi(j−1)(k−1)/N]as the eigen-matrix.Then the exact solution at timetis calculated w ith discrete Fourier transform, namely,

whereΛ(t)refers to the diagonal matrix taking the eigenvalues exp(−λjt),and the initial vector x(0)att=0 is fi lled w ith standard norm al Gaussian numbers. The errors of the numerical integration are defined asthe Euclidean distancebetween the exact result and the Runge-Kutta approximation x(appr)(t)

In FIG.2,wepresent theerror ofall themethodsw ith the step size of2n×10−3(0≤n≤12).The results clearly show that there exists a region that the accuracy w ill be com parable w ith the machine round-off errors and they cannot be reduced further simply by decreasing the step size.There also exists an unstable region for the each integrator,when the step size is bigger than a critical value,the error increases explosively.The error corresponding to this critical step size is about one,which im p lies that the step size in this p roblem is controlled by the accuracy.

Not shown in FIG.2,wealso carry out the integration w ith the classical fourth-order Runge-Kutta method, the Fehlberg’s 4(5)m ethod[45],the Calvo 6(5)m ethod [46],and the Prince-Dormand 8(7)method[47].The 4th-order LSRK is identical to the classical fourth-order Runge-Kutta method in errors as well as the stability. Higher-order LSRK schem es are faster than the sam e order classical RK ones,but no more than a factor of two.However,a high-order LSRK(p>6)scheme can be two or three times faster than the corresponding classical Runge-Kuttam ethod which assum es the same number of stages.

Provided w ith a pre-selected error tolerance,higher ordermethods normally can take a bigger step size.For exam p le,w ith the error tolerance of 1×10−5,the step size can be as large as 0.31 for LSRK 12-12,but only 0.0083 for LSRK 4-4.If the error tolerance deceases by two orders,the corresponding step size becom es0.21 for LSRK 12-12 and 0.0026 for LSRK 4-4.Roughly speaking,thespeedup w illbe9 to 21 timesifweusea 12-stage 12-orderm ethod instead of a fourth-order Runge-Kutta schem e.O ther highers-stages-orderm ethods also depict an increase of 3−20 times in effi ciency,depending on the order and the error tolerance.W ith the same step size in the stable region,the stability-optim ized m ethod LSRK 13-8(5)always produces smaller errors than the same ordermethod LSRK 8-8.W hen the additional stagesof the LSRK 13-8(5)scheme are taken into account,the LSRK 13-8(5)scheme has the sam e perform ance as the LSRK 8-8 method.

Again,to have comparison,we consider the LSRK 13-4(4)scheme which is designed for general differential equations.W ith a given step size,it system atically produces smaller errors than that of the LSRK 4-4 method.However,provided thesameaccuracy(1×10−5or 1×10−7),the increase in the step size is less than three times.We thus conclude that the LSRK 13-4(4) scheme is less effi cient than the LSRK 4-4method.

III.APPLICAT IONS TO SIMULATION OF QUANTUM COHERENCE IN BIOLOGICAL-RELEVANT SYSTEM S

We w ill app ly the high-order low-storage m ethods outlined in Section IIto amore realistic simulation,that is,the non-Markovian exciton diff usion in a linear aggregate.Theexciton dynam icsp laysa crucial role in organic photovoltaic cell[48],photosynthesis[49],organic light-em itting diodes[50]and many others[51,52].In such systems,energy transfers to the reaction center along w ith the exciton diff usion.The natural photosynthesis system s exhibit surprisingly high effi ciency in the exciton energy transfer[49].The theoretical understanding of the design features and the underline mechanism w ill greatly help the artifi cial solar energy app lication[53].We are still at the very beginning to understand the processand now expect that dissipation caused by the environm ent is the key for effi cient energy transport[49,53].

Here,we study the noise-assisted exciton diff usion in an artificial chain consisting ofN,N′-Di(N-(2-am inoethyl)-benzam ide)-1,6,7,12-tetra(4-tert-butylphenoxy)-3,4:9,10-tetracarboxylic acid bisim ide(PBI-1). Perylene bism ide system s are very flexible in tuning the geometric configuration as well as the electronic properties in a large extent upon proper substitution or diff erent solvation[54−56].Therefore,they serve as idea test system s for theories to investigate the eff ects of various factors on exciton transfer and have been extensively studied in the literature.For examp le, Wolter and co-workers investigated the size-dependent exciton dynam ics in one-dimensional PBI aggregates [57].A maximal diffusion length was estimated to be 96 nm,which is about 10 tim es larger than that in disordered polymers[58].The photophysical properties of PBImonomer and oligomers have been calculatedin Refs.[59,60].

As the time scale of the exciton dynam ics is of the order of picosecond and longer,it is com parable to that of the environment.Thus,the dynam ics is essentially non-Markovian.Here,the non-Markovian dynam ics is investigated w ith the hierarchical approach,which introduces auxiliary density matrices to characterize the eff ect of the bath and is one of a few available methods to study non-Markovianity[6−9,61−64].Usually,the number of auxiliary matrices are large and the requirem ent for the m em ory storage is high at low tem peratures.The low-storage scheme thus is ideal to integrate the hierarchical equations.

A.H ierarchical approach for dissipative exciton dynam ics

In order to focus on the exciton diffusion,the exciton is supposed to be stable and its dissociation as well as the decay is not taken into account.M oreover,we only investigate the dynam ics in the one-exciton manifold, that is,only the diabatic electronic state corresponding to theS0state is included.W ith this restriction,there are onlyNaggone-exciton states|m⟩(only sitemis inS1)for anNagg-site aggregate.

The Ham iltonian for the exciton diff usion under dissipation can be approximated w ith the system-p lus-bath m odel,

whereJmnis the coup ling constant.

W hen the size of the monom er in the aggregate is big, each site may have a set of localized vibrations coup led to the exciton at this site only.Operatorin Eq.(34) stands for the exciton-vibration coup ling at sitemand a linear form is accepted[65]

For linear dissipation,the effect of the vibrations on the exciton dynam ics is fully captured by the bath response function,

which in turn is determ ined by the spectral density functionJ(ω)[62,64].Here,β=1/kBTis the inverse tem perature rescaled by the Boltzmann constantkB.Generally,the spectral density function can be expressed in terms of the coup ling constant,reduced mass,and frequency of each oscillator in the bath.But most often theoretical simu lations adopt special form to model the bath.The localized vibrational modes normally have high harmonic frequencies and can be described by a dam ped Brownian oscillator bath m odel [66]

whereηis the dissipation strength,ω0the centralmode frequency,andΛthe cutoff.Below we w ill refer to the under-dam ped lim it,where the central mode frequencyω0is much larger than the cutoffΛ.

For this spectral density the response function can be expanded into a sum of exponentials as follows[62]

whereΩ1,2=Λ/2±iξw ithandΩk+2=2πk/(~β)(k>0)is thekth M atsubara frequency of the bath.In w riting this equation,we have separated the short memory contribution of the response function w ith the smallest integerNesatisfyingΩNe≫Λ.

W ith the exponential decom position of the response functions,the hierarchicalapproach is effi cient to propagate the reduced density matrix for the description of the non-Markovian exciton dynam ics.The hierarchical approach is a schem e to attack the non-perturbative non-M arkovian dynam ics for the open quantum system by converting the memory of the bath into a series of auxiliary density matrices.Here we w ill follow the stochastic decoup ling procedure of Shao and coworkers [9,67].Derivation based on fundamental Itˆo calculus can be also found somewhere else[68,69].The resulting equations are as follows

αj;±are defined as:

w ithj′being the index satisfying.Them atrixAhas the dimension ofNagg×Neand itsmatrix elementAm,jrefers to thej-th exponential in the response function of the bath for them-th site.The operator“+/−”is defined as.In Eq.(40),Ris the Red field superoperator which app roxim ates the effect of the shortmemory part of the response function in Eq.(39)[62].

B.Results and discussions

As far as the present m odel is concerned,we w ill adopt the same parameters as those in Ref.[62].Specifically,in the excitonic part we w ill only consider the nearest neighbor coup ling and use the calculated valueJmn=−514 cm−1[60].For the bath,we assume all of them are independent and identical and take the theoretical calculated valuesω0=1415 cm−1andη=0.44 in the spectral density function Eq.(32)[59].The simulations are carried out at a tem perature of 300 K w ith an initial conditionρ(0)=|1⟩⟨1|.

FIG.3 Exciton population evolution at site one.

W ith these settings,there are totally 1,674,992 auxiliary density matrices in the hierarchy and it takes about 5.6 GB memory to store all comp lex matrices.A LSRK method w ith any order w ill need 11.2 GB memory storage for the num erical propagation.By contrast, the classical fourth-order Runge-Kuttamethod requires 16.8 GB space;the fi fth-order Runge-Kutta-Fehlberg method needs 39.2 GB[45].

The numerical simulations are done w ith the Hyshe package[70],which em p loys a shared memory parallelization.Specifically,we use a 24-core Intel(R) Xeon(R)E5-2650v3 CPU w ith a clock speed of 2.30 GHz for the calculations.The effi ciency is estimated as about 0.75 for this simulation.W ith the LSRK 12-12 scheme,the dynam ics can be propagated using a step size of 100 a.u.(∼2.42 fs).As such,the numerical integration takes 6.5 h for 1000 steps.

Theexciton population at siteone isshown in FIG.3. The population follows a steady equilibration after an initial quick drop.Note that there exist three clear spikesaround 0.24,0.47,0.69 ps,andm inor onesaround 0.92,1.2,and 1.5 ps.The peak heights decrease as time elapses,but the peak separation is almost equal in tim e.Onem ay ask why thisquasi-periodic phenom enon emerges?

To answer this question,we scrutinize the exciton population at all sites and present the data in FIG.4. The results reveal that the exciton transfers in the aggregatew ith a constant rate and is bounced by the two ends of the aggregate.Together w ith the dissipation induced by the vibrations,these two factors lead to an under-dam ped exciton m otion.Based on this observation,the exciton hopping rate is about 123.4 ps−1, which is significantly larger than 7.8±2.7 ps−1,a value measured w ith ultrafast transient absorption spectroscopy[56].The reason of this gap is because we are using the theoretically calculated bare Coulomb coup lingJmnof−514 cm−1,which is larger than the experim ental value(120 cm−1)[56].A closer exam ination of the population at each time reveals that the exciton population spread in theaggregate and thus theexciton transfer as a wave.O f course,the environment causes dissipation and the exciton diff uses along the aggregate irreversibly.

FIG.4 Exciton population evolution of the aggregate.The region[n−1/2,n+1/2]along they-axis refers to the population at Siten.

FIG.5 Coherence beating in the aggregate.The region [n,n+1]along they-axis refers to the absolute value of the off-diagonal m atrix elem entsρn,n+1(t)of the reduced density matrix.

Thewave-likemotion of the exciton transfer can also be m anifested by the coherence beating.FIG.5 p lots the absolute value of the off-diagonalmatrix elementsρn,n+1(t)(1≤n≤14)of the reduced density matrix at different times.The coherence exhibits a sim ilar pattern as the popu lation.As time elapses,the coherence doesnot vanish and approachesa constant instead, which reflects the exciton entanglement in equilibrium.

IV.SUMM ARY

Previously,low-storage Runge-Kutta methods were designed for integration of general nonlinear or tim evariant systems,which is effi cient in the use of the memory storage.To generate a low-storage scheme for generalequations,we have to take additional stages because there aremore restrictions for themethod.However,the dynam ics ofmany quantum system s is governed by linear differential equationsw ith constant coeffi cients.For these autonomous linear equations,there are on lyprestrictions versus 2p−1 param eters for ap-stage low-storage m ethod.Thus,the order barrier does not app ly for Runge-Kuttamethods in solving autonomous linear equations.The system ofequations for the parameters in the current scheme can be solved analytically and it is straightforward to design a new LSRK method.Thisallowsus to obtainp-stagepth-order lowstoragemethodswhich require lessstages than the classicalpth-order Runge-Kutta schem es whenpis larger than four.W hile the general LSRK method consumes longer com putational time than the same order classical Runge-Kutta scheme,the currentmethod takes less com putational tim e for autonomous linear equations. In thiswork,we have introduced a series of low-storage schemes specific for autonomous linear equations.

Numerics shows that the 12-stage 12th-order LSRK method could be an order ofmagnitude faster than the fourth-order Runge-Kuttam ethod.

As an exam p le,we have app lied the suggested 12-stage12th-orderm ethod for thehierarchicalsimulations of the exciton diff usion in a 15-site linear aggregate, where it takes 5.6GBmemory to store the reduced and auxiliary density m atrices.Thanks to the effi cient lowstage schem e,we can propagate the non-M arkovian dynam icsup to 2.4 psw ithin 6.5 h w ith a 24-core Intel(R) Xeon(R)E5-2650v3 CPU.

V.ACKNOW LEDGM ENTS

This work is supported by the National Natural Science Foundation of China(No.21373064),the Program for Innovative Research Team of Guizhou Province(No.QKTD[2014]4021),and the Natural Scientific Foundation from Guizhou Provincial Department of Education(No.ZDXK[2014]18).A ll the calculations were performed at Guizhou Provincial High-Perform ance Com puting Center of Condensed M aterials and Molecular Simulation in Guizhou Education University.

[1]E.Schr¨odinger,Ann.Phys.384,361(1926).

[2]M.Born,W.Heisenberg,and P.Jordan,Zeit.Phys.35, 557(1925).

[3]F.B loch,Phys.Rev.70,460(1946).

[4]A.G.Red field,IBM J.Res.Dev.1,19(1957).

[5]G.Lindblad,Commun.M ath.Phys.48,119(1976).

[6]V.Shapiro and V.Loginov,Physica A 91,563(1978).

[7]Y.Tanimura and R.Kubo,J.Phys.Soc.Jpn.58,101 (1989).

[8]Y.Tanimura,Phys.Rev.A 41,6676(1990).

[9]Y.A.Yan,F.Yang,Y.Liu,and J.Shao,Chem.Phys. Lett.395,216(2004).

[10]L.Adrianova,Introduction to Linear Systems ofDiffi erential Equations,Translations of Mathematical M onographs,Vol.146,Providence,Rhode Island:AMS, (1995).

[11]A.G reenbaum,Iterative M ethods for Solving Linear Systems,Frontiers in Applied Mathematics,Vol.17, Philadelphia:SIAM,(1997).

[12]E.Coddington and R.Carlson,Linear Ordinary Diffi erential Equations.Philadelphia:SIAM,(1997).

[13]A.Forsyth,Part 3.Ordinary Linear Equations,Theory of Diffi erential Equations Vol.4,Cambridge:Cambridge University Press(1902).

[14]Y.Saad,Iterative M ethods for Sparse Linear System s 2nd Edn.,Philadelphia:SIAM,(2003).

[15]C.Lanczos,J.Res.Nat.Bur.Standards 45,255(1950).

[16]R.Haydock,Com put.Phys.Commun.20,11(1980).

[17]V.D ruskin and L.Knizhnerm an,J.Com put.App l. Math.50,255(1994).

[18]H.Tal-Ezer and R.Kosloff,J.Chem.Phys.81,3967 (1984).

[19]C.Leforestier,R.Bisseling,C.Cerjan,M.Feit,R.Friesner,A.Guldberg,A.Hamm erich,G.Jolicard,W.Karrlein,H.Meyer,N.Lipkin,O.Roncero,and R.Kosloff, J.Com put.Phys.94,59(1991).

[20]R.C.W ard,SIAM J.Num er.Anal.14,600(1977).

[21]M.Hochbruck,C.Lubich,and H.Selhofer,SIAM J. Sci.Com put.19,1552(1998).

[22]R.B.Sid je,ACM Trans.Math.Softw.24,130(1998).

[23]J.C.Butcher,Num ericalM ethods for Ordinary Diffi erential Equations,2nd Edn.Chichester:W iley,(2008).

[24]R.E.Bellm an,Dynam ic Programm ing,Princeton: Princeton University Press,(1957).

[25]M.Feit,J.Fleck,and A.Steiger,J.Com put.Phys.47, 412(1982).

[26]P.Houwen,Numer.Math.20,149(1972).

[27]P.J.Van Der Houwen,Construction of Integration Form ulas for Initial Value Problem s,Vol.19,Am sterdam: North-Holland,(1977).

[28]A.W ray,NASA Am es Research Center,California: Moff ett Field,202(1990).

[29]J.W illiam son,J.Com put.Phys.35,48(1980).

[30]S.Gill,M ath.Proc.Cambridge Phil.Soc.47,96(1951).

[31]E.K.B lum,M ath.Com put.16,176(1962).

[32]D.J.Fy fe,M ath.Com put.20,392(1966).

[33]D.I.Ketcheson,SIAM J.Sci.Com put.30,2113(2008).

[34]M.H.Carpenter and C.A.Kennedy,Fourth-order 2Nstorage Runge-Kutta schem es,Tech.Rep.NASATM-109112,VA:National Aeronautics and Space Adm inistration,NASA Langley Research Center,(1994).

[35]C.A.Kennedy,M.H.Carpenter,and R.Lew is,App l. Num er.M ath.35,177(2000).

[36]J.Butcher,BIT Numer.Math.25,521(1985).

[37]J.Niegem ann,R.Diehl,and K.Busch,J.Com put. Phys.231,364(2012).

[38]T.Becker,H.K redel,and V.Weisp fenning,Gr¨obner Bases:a Com putational Approach to Commutative A lgebra,Corrected Edn.,Graduate Texts in Mathematics, Vol.141,London,UK:Springer-Verlag,(1993).

[39]D.Zingg and T.Chisholm,App l.Numer.M ath.31, 227(1999).

[40]J.A.Nelder and R.M ead,Com put.J.7,308(1965).

[41]S.Johnson,The Nlopt Non linear-Optim ization Package.

[42]K.Levenberg,Quart.J.App l.M ath.2,164(1944).

[43]D.W.Marquard t,J.Soc.Indust.App l.Math.11,431 (1963).

[44]M.Lourakis,Levm ar:Levenberg-Marquardt Nonlinear Least Squares A lgorithm s in C/C++,http://www.ics. forth.gr/∼lourakis/levm ar/July(2004).

[45]E.Fehlberg,Com puting 6,61(1970).

[46]M.Calvo,J.M ontijano,and L.Randez,Com put.M ath. App l.20,15(1990).

[47]P.Prince and J.Dorm and,J.Com put.App l.M ath.7, 67(1981).

[48]C.W.Tang,App l.Phys.Lett.48,183(1986).

[49]H.Lee,Y.C.Cheng,and G.R.Flem ing,Science 316, 1462(2007).

[50]J.H.Burroughes,D.D.C.Brad ley,A.R.B rown,R. N.M arks,K.M ackay,R.H.Friend,P.L.Burns,and A.B.Holmes,Nature 347,539(1990).

[51]G.Trinkunas,J.L.Herek,T.Pol´ıvka,V.Sundstr¨om, and T.Pullerits,Phys.Rev.Lett.86,4167(2001).

[52]V.M ay and O.K¨uhn,Charge and Energy Transfer Dynam ics in M olecular System s,3rd Edn.Weinheim: W ILEY-VCH,(2010).

[53]M.delRey,A.W.Chin,S.F.Huelga,and M.B.Plenio, J.Phys.Chem.Lett.4,903(2013).

[54]Z.Chen,V.Stepanenko,V.Dehm,P.Prins,L.Siebbeles,J.Seibt,P.Marquetand,V.Engel,and F.W rthner, Chem.Eur.J.13,436(2007).

[55]X.Q.Li,X.Zhang,S.Ghosh,and F.W¨urthner,Chem. Eur.J.14,8074(2008).

[56]M.R.Wasielewski,Acc.Chem.Res.42,1910(2009).

[57]S.Wolter,J.Aizezers,F.Fennel,M.Seidel,F.Wurthner, O.K¨uhn,and S.Lochbrunner,New J.Phys.14,105027 (2012).

[58]H.Marciniak,X.Q.Li,F.W¨urthner,and S.Lochbrunner,J.Phys.Chem.A 115,648(2011).

[59]D.Ambrosek,H.Marciniak,S.Lochbrunner,J. Tatchen,X.Li,F.W¨urthner,and O.K¨uhn,Phys. Chem.Chem.Phys.13,17649(2011).

[60]D.Ambrosek,A.K¨ohn,J.Schulze,and O.K¨uhn,J. Phys.Chem.A 116,11451(2012).

[61]Q.Shi,L.Chen,G.Nan,R.X.Xu,and Y.Yan,J. Chem.Phys.130,084105(2009).

[62]Y.A.Yan and O.K¨uhn,New J.Phys.14,105004 (2012).

[63]Y.A.Yan and S.Cai,J.Chem.Phys.141,054105 (2014).

[64]Y.Yan,J.Chem.Phys.144,024305(2016).

[65]O.K¨uhn and V.Sundstr¨om,J.Chem.Phys.107,4154 (1997).

[66]H.G rabert,P.Schramm,and G.L.Ingold,Phys.Rep. 168,115(1988).

[67]J.Shao,J.Chem.Phys.120,5053(2004).

[68]M.Schr¨oter,S.Ivanov,J.Schulze,S.Polyutov,Y.Yan, T.Pu llerits,and O.K¨uhn,Phys.Rep.567,1(2015).

[69]Y.Yan and J.Shao,Front.Phys.11,110309(2016).

[70]Y.A.Yan and Y.Zhou,Hybrid Stochastic-Hierarchical Equations.http://nano.gznc.edu.cn/∼yunan/hyshe. htm l(2012).

ceived on March 3,2017;Accepted on March 17,2017)

∗Author to whom correspondence shou ld be addressed.E-m ail: yunan@gznc.edu.cn