# Difference between wave vector and density matrix in numerical calculation of Schrödinger equation

I solved Schrödinger equation for a following tow-level time-dependent Hamiltonian numerically in two ways:

import numpy as np

def H(t):
return np.array([[t,0.5],[0.5,-t]],dtype="complex64")

First, the state is treated as wave vector:

from scipy.integrate import solve_ivp
def schrodinger(t, X):
dXdt = -1j * H(t).dot(X)
return dXdt

T=80.
t_list = np.arange(-T, T, 2*T / 1000.0)
init_state= np.array([1.,0.],dtype="complex64")
solution = solve_ivp(fun=lambda t, X: schrodinger(t, X), t_span=[-T,T], y0=init_state, t_eval=t_list, method="RK45", vectorized=True)

Second, the state is treated as (vectorized) density matrix:

def schrodinger_rho(t, X):
unit = np.eye(2, 2, dtype="complex64")
dXdt = -1j * (np.kron(unit, H(t)) - np.kron(H(t).T, unit)).dot(X)
return dXdt
T=80.
t_list = np.arange(-T, T, 2*T / 1000.0)
init_state= np.array([1.,0.,0.,0.],dtype="complex64")
solution_rho = solve_ivp(fun=lambda t, X: schrodinger_rho(t, X), t_span=[-T,T], y0=init_state, t_eval=t_list, method="RK45", vectorized=True)

The numerical results are as follows:

Here is a question : why cannot the state preserve its unitarity, say $$P_0+P_1=1$$, when the state is treated as wave vector? The larger the value of T, the stronger this tendency of violation becomes. Moreover, the numerical results become worse when I use method="BDF".

I want to treat the state as wave vector because of the numerical cost. Is there any way to improve this phenomenon?

• Explicit Runge Kutta methods rarely conserves integrals of motion. Did you try with an implicit Gauss Runge Kutta method? Mar 3 at 22:50
• @G. Fougeron But isn't it true that for a smaller integration time step the integrals of motion are conserved better, even if the method is not conservative? So for a small enough time step you'll get all conservation laws that the underlying PDE has? Mar 4 at 6:28
• This is true for the short time behavior, but not necessarily for long-time behavior. For more precise info on the subject, I recommend the book by Hairer : Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations. If you can't find it, then maybe chackout his list of preprints unige.ch/~hairer/preprints.html Mar 4 at 10:21
• @G.Fougeron I did not try it because there is no choice of implicit Runge Kutta method in complex domain in scipy. Instead, I found odeintw where LSODA can be applied in complex domain. This method works. Mar 4 at 13:38
• There you go. Then it makes perfect sense. If you check out the resources I linked, you will find proofs / explanations for this 😊 Mar 5 at 8:01