Hi Simone,

The method proposed by Kosloff et al. [see, e.g., Ann. Rev. Phys. Chem. 45, 145 (1994); J. Chem. Phys. 110, 5538 (1999), and references therein], which is based on an expansion of the time evolution super-operator as a convergent series of either Chebychev or Faber polynomials, has been very useful for me.

Cheers,

Alexis

How big is the system you are considering? I had similar problems but I was dealing with the dissipative dynamics with Lindblad terms. But in my case the system was small enough that I was able to diagonalize explicitly the linear evolution of the density matrix.

Andrea, I have to introduce Lindblad terms too later. But the thing is - I already have the diagonalization, for all the possible values of the Hamiltonian. So I was in fact thinking about just using exponential propagators. The system's size varies depending on the accuracy I require - there's an estimator for the grid spacing that goes with the deBroglie wavelength of the maximum energy of interest - but seldom goes beyond a 100x100 matrix.

Ok without the Lindblad term, you should not have problems in working in the basis where the Hamiltonian is diagonal. Then the time evolution operators are also diagonal and you can achieve arbitrary times without problems, by explicitly solving the Liouville equation.

In presence of the Lindblad term, the situation is more complicate. In principle, the density matrix time derivative is still linear in the density matrix:

drho/dt = -i[H,rho] + L rho = S rho

L is the dissipative part, and S is called usually superoperator because it is a linear operator acting on the space of operators, rho in this case. The main problem is that the size of S is quadratic in the size of the system. So if you take H 100x100, the space of operators has dimension 10000 and S is 10000x10000. This is still doable, but some effort is needed.

The alternative is to implement the dynamics. But it depends on the sizes and on the times you want to achieve. For example, if you just need the stationary state (t->inf), it can be done more easily.

Ok, but in that case I can recast in the local basis set of wavefunctions and use those, right? So since I can describe the system with say, the 10 or so lowest states, the superoperator doesn't get bigger than 100x100 and it's still feasible.

Another issue I am dealing with is time dependent Hamiltonians. I suppose I should be able to apply them by just multiplying subsequent steps dt, right?

@Simone

A second order, very stable, method that I used a lot for 1D to 2D quantum dynamical problems with and without time-dependent Hamilton operator is what I call the asynchronous leapfrog method. It works as follows: In addition to rho you use a second variable of the same data type, which (for reasons that become clear immediately) we will call rhoDot. rhoDot has to be initialized:

rhoDot(t_0) = 1/(ih)*[H(t_0), rho(t_0)] where we have assumed a time dependent Hamiltonian and assumed that we know the initial value of rho. From that step on the two quantities rho and rhoDot are in the game and the task for a time stepping algorithm is to propagate them both. The normal second order leapfrog algorithm also works with two state data rho(t_i) and rho(t_(i+1)). Identifying rhoDot(t_i) with (rho(t_(i+1)) - rho(t_(i-1)) )/(t_(i+1) - t_(i-1)) suggests the following simple and symmetric time step for the state (t,rho,rhoDot):

t += dt/2

rho += rhoDot*dt/2

rhoDot += ( 1/(ih)*[H(t), rho] - rhoDot )*2

rho += rhoDot*dt/2

t += dt/2

One of the beauties od the algorithm is that never a component has to be memorized.

The value changes occur in a natural order. To make things completely explicit let

us consider the central step rhoDot += ..... When we evaluate the term which according to this formula has to be added to rhoDot we use the old value of rhoDot on the right-hand side. The value t in H(t) is, due to the first step t += dt/2 the value of t at 'the midle of the step' (midpoint method). In my C++ programs I have the operators += for the state components implemented and the lines above are virtually identical to the actual code. More about the algorithm you will find in the appended

link.

http://www.arxiv.org/abs/1311.6602

That is an interesting idea, but I'm not sure it would change much compared to my previous algorithm - which was basically

t += dt

rhoDot = 1/(ih)*[H(t), rho]

rho += rhoDot*dt

Of course it's not as refined, but I don't think this difference justifies a disastrous outcome such as the one I experienced. And I also did a Runge-Kutta version of this integration, which I think should be of higher order than a leapfrog. However, the fact that you say that this algorithm was good enough for you sparkles my interest - if that is the case, maybe there was something else going wrong in my case. I'll take a look at the paper, thanks!

@Simone,

your algorithm is known to me and is much less stable (in technical terms all nonzero points of the imaginary axis do not belong to the region of absolute stability, whereas for aynchronous leapfrog a whole interval on the imaginary axis belongs to that region)

than asynchronous leapfrog. You will see a drastic improvement of the stability behaviour for your system, I'm sure!

Ok, thanks, I'll try!

The easiest and more direct way to deal with Liouville-Von Neumann equations is the use of the evolution operator method.

The solution can be cast in the form

ro(x,p,t)=e^(Lt) ro(x,p,0)

I suppose that the Liouville Von Neumann operator is not explicitly time dependent (otherwhise time ordering methods should be used) The method is illustrated in a recent paper by Sabia and myselg (Dattoli) in Phys. Rev. A.

The paper refers to previous publications and has been exploited to treat the classical relativistic harmonic oscillator, the inclusion of the quantum operators does not create significant problems

Pino

...Regarding the time dependence of the Hamiltonian or of the Liouville operator, the evolution operator expansion can bedone using a Magnus type expansion.

Pino

Hi Simone,

The paper [M. Ndong et al., J. Chem. Phys. 132, 064105 (2010)] might be useful to you. The suggested idea can be incorporated into the Chebychev or Faber polynomials based propagators for the LvN equation [e.g., Berman & Kosloff, Comp. Phys. Comm. 63, 1 (1991)].

Cheers,

Alexis

Thanks to everyone for their answers, a quick update:

1) while it does improve the situation a bit, the quantum leapfrog method unfortunately ends up failing in a similar way to the simple ones I had mentioned before;

2) the Chebychev expansion method suggested by Alexis looks very promising but is also very complicated for me. I'll keep it in store but I'd like to begin with something more "quick and dirty" for now - unfortunately I don't have the time I would need to implement something like that for now;

3) for now, the method mentioned by Giuseppe and Andrea is the one I am using and works fine without time dependence. I have to introduce, however, explicit time dependence. I have considered a Magnus expansion, but I am worried by the fact that I've heard matrix powers tend to become numerically unstable very quickly.

One thing: while there IS a time dependence, it is of a very simple form, as the only time-dependent elements are on the diagonal of the Hamiltonian. Could this lead to some simplified analytical solution?

Simone,

if the Hamiltonian is does not depend on time and its discretized version has the operator norm ||H|| the time step dt needs to satisfy dt < 1/||H||. Otherwise, sooner or later, the solution will explode. For time-dependent Hamiltonian things are not so simple. Replacing ||H|| by sup_t ||H(t)|| sugests itself and is certainly not far from correct. It should take you not too much time to run a series with timepsteps reduced by a factor 0.5 and to see when the paradise of stability will be reached. The point is that Runge-Kutta and your simpler integrator will never reach stability for all times, wheras any leapfrog becomes absolutely stable for sufficiently small timesteps.

The problem for me is that after trying a bit it seems to be a very short timestep, which will make the calculation too slow for the timescales I'm interested in.

Simone,

having been in a similar situation many times, I propose the following strategy.

1. Test your Hamiltonian by studying the dynamics not with density operators but with pure states psi:

t += dt/2

psi += psiDot*dt/2

psiDot += ( 1/(ih)*H(t)*psi - psiDot )*2

psi += psiDot*dt/2

t += dt/2

Of course, initialize as psiDot(t_0) = 1/(ih)*H(t_0)*psi(t_0).

The norm of psi has to be constant up to some tiny wiggling under that dynamics. For time-non-dependent Hamiltonians also the energy has a constant real part and a tiny wiggling imaginary part. Any coding error is in the Hamiltonian will spoil these properties virtually certainly from the beginning. If the norm of the state is virtually constant for some time and explodes then, this is says that the timestep has to be reduced (I speak always about constant time step, not adaptive time stepping which never worked for me in quantum dynamical problems) If you know the eigenstates of your initial rho you probably will do the test with one of the eigenstates of highest weight. This will give you also a good idea what the acceptable time steps are since whenever your problem is feasible on your computing resources for some mixture rho it will certainly be feasible for pure states and you will not have to wait too long for visual results.

If you have implemented H it is easy to compute ||H|| the role of which was explained in my last contribution.

See http://www.ulrichmutze.de/articles/qmadd11.pdf equation (80).

If you got comfortable with the dynamics of pure states you could implement the dynamic for mixtures by spectral decomposition. Determine the eigenvectors of rho,

propagate those over time, and assemble the propagated rho from the propagated eigenvectors. If there are many eigenvectors this might take a long time too but you are in better control of what happens compared to listening to the dynamics of the huge object rho. Probably it is an acceptable approximation to ignore many eigenvectors if they contribute to rho only with a small weight.

Dear Simone,

I do not know what is the form of your time dependence so I can't be helpful to tguess a specific transformation which can eliminate it. Anyhow the specific parameters may be such that you can keep only few terms in the Magnus expansion. A further possibility is the use of the Fer expansion, I prefer it being more manageable that Magnus.

All the best

Pino

when I say few terms I mean that you keep few term in the commutator chain and leave allways the exponential form of the evolution operator, so you are sure that the unitarity is preserved.

Pino