# Runge Kutta Method for a Lindblad Equation

I am solving a Lindblad equation for a dissipative Harmonic Oscillator. My Hamiltonian is time dependent,

My Lindblad Equation can be written as $$\frac{d\rho}{dt}=\frac{[H(t),\rho]}{i\hbar}+D(\rho,a)$$ where the last part $D(\rho,a)$ represents the Linbald operator due to the dissipation.

Considering Hamiltonian is constant over the interval $t_{1}$ to $t_{1}+dt$. I used the Runge-Kutta fourth order (RK4) method in each time interval and integrated the above equation. I would like to know whether my procedure could be justified? I have seen methods which refer to the interaction picture then apply the RK4 method. I think the method that I mentioned is even more simple, but I would like to know its validity. My mathematical justification is the following \begin{eqnarray} U(t)=\exp{\left(-\frac{i}{\hbar}\int_{0}^{t}H(\tau)d\tau\right)}\\ U(t,t+dt)=\exp{\left(-\frac{i}{\hbar}\int_{t}^{t+dt}H(\tau)d\tau\right)}\\ \end{eqnarray} Now I suppose, in the given tiny interval my Hamiltonian is constant, hence I can move it outside the integral and I get, \begin{eqnarray} U(t,t+dt)=\exp{\left(-\frac{i}{\hbar}H(t)\int_{t}^{t+dt}d\tau\right)}\\ U(t,t+dt)=\exp{\left(-\frac{i}{\hbar}H(t)dt\right)}\\ \end{eqnarray} Now I can use the Taylor expansion and use the Runge-Kutta to integrate the function. But every time the $U(t,t+dt)$ operator changes. Is my method legitimate?

• Your method will violate unitarity in the limit of $D \rightarrow 0$, so if you are worried about that you may want to come up with a unitary integrator with a $D$ kick in the middle. At least, it looks like it might. Runge-Kutta is not intrinsically unitary, although there are unitary Runge-Kutta type algorithms. – webb May 16 '14 at 16:14
• Method's wrong. Violates unitarity. Try using the Fritzpack's modified Runge Kutta method. – user47333 May 27 '14 at 11:28
• Hi @JAMES Could you please provide me a link to the Fritzpack's method. I googled it but I didn't find anything related to that. – Sijo Joseph May 27 '14 at 15:18
• Do you get the correct method now? The method you described is not RK4, and it should be of order one only because you haven't evaluated at the middle point U(t,t+dt/2) – unsym Jul 20 '14 at 7:04

Hi guys I found the answer. The answer is Magnus expansion method. My method is not a good approximation. The method that I described is valid only in a very very small interval of time. Hence the numerical algorithm will be slow and inaccurate. There is a truncation error due to the Zassenhaus formula $$e^{-i/\hbar(t_{1}+t_{2})(H(t_{1})+H(t_{2})}= e^{-i/\hbar\,(t_{1}+t_{2})H(t_{1})} e^{-i/\hbar\,(t_{1}+t_{2})H(t_{2})} e^{1/\hbar^2\frac{(t_{1}+t_{2})^2}{2} [H(t_{1}),H(t_{2})]} \cdots$$. The last additional error factor comes since we deal with time dependent hamiltonian and the Hamiltonian of different time doesn't commute each other. Hence we have to use the Magnus expansion, which is aimed to find an equivalent matrix as a linear combinations of Hamiltonians at different time and the linear coefficients are cleverly chosen to eliminate the error due to Zassenhaus formula. I found this article very useful.
That appears to be a very bad thing to do. The LHS measures how a quantity changes with time. The only time-dependent quantity on the RHS is the Hamiltonian. You hold it time-independent, for your small interval, and integrate $d\rho/dt = f(\rho)$. For the next interval, you again arrest the time-dependence. So, in totality, between these two steps, your Hamiltonian hasn't advanced in time at all, whereas a truly time-dependent Hamiltonian should have!
Yeah! That's right. The method's fine for a time independent Hamiltonian, but not for an $H(t)$.