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
\begin{equation}
\frac{d\rho}{dt}=\frac{[H(t),\rho]}{i\hbar}+D(\rho,a)
\end{equation}
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?
 A: 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.
A: 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!
A: Yeah! That's right. The method's fine for a time independent Hamiltonian, but not for an $H(t)$.  
