Numerical Error while extracting pure state from density matrix.

[ADITYA DEV]

Hi all!

I'm trying to simulate a two level system coupled with a High Spin Environment. The numerical code for same is given below.

The `psi_dm` (i.e. the density matrix in psi) comes out as a pure state for most time steps (which is what I expect). But for time, t = (2n + 1)pi/, psi is no more a pure state. This shouldn't be a case, as by definition it's a pure state.

I can't see where the problem is. I tried looking in the documentation and tried figuring out if the problem is in eigenvale solver or any other function.

Any insight or comments are appreciated.

PS. The global state is chosen at random i.e for the example code it's 20, you may chose any other.

Thanks for help in advance.

import qutip as qt

import numpy as np

import matplotlib.pyplot as plt

omega = 1

V0 = 1

J = 50

a = 20

Hs = qt.Qobj(np.zeros([2, 2]), dims=qt.sigmaz().dims, )

Hc = qt.Qobj(omega* qt.jmat(J, "z"))

Vs = V0 *(qt.sigmaz() + qt.sigmax())

Vc = qt.jmat(J, "x")

dimSys = Hs.shape[0]

dimEnv = Hc.shape[0]

def Xi_gauss(t, J, a):

state_ = np.exp(-J**2/a**2)*qt.spin_state(J, -J)*np.exp(-1j*J*t*omega)

for m in range(-J+1, J+1):

state_ += np.exp(-m**2/a**2)*np.exp(-1j*m*t*omega)*qt.spin_state(J, m)

return state_.unit()

H = qt.tensor(Hs, qt.qeye(dimEnv))+ qt.tensor(qt.qeye(dimSys), Hc) + qt.tensor(Vs, Vc)

eigenenergies, eigenstates = H.eigenstates()

glob_ket = eigenstates[20]

def psi_lambda(t, J, a):

"""

makes a density matrix of the form |psi(t)> <psi(t)|

then traces out the env state

then finds the eigenstate corresponding to the maximum eigenvalue

that's the state of the system

"""

psi_dm = qt.tensor(qt.qeye(dimSys), Xi_gauss(t, J, a).proj()) * glob_ket.proj().tidyup()

psi_dm = (psi_dm.unit()).ptrace(0).tidyup().unit() ## we need to keep system 0 after tracing

valeig_, valstate_ = psi_dm.eigenstates()

index_of_max_eigenvalue = np.abs(np.round(valeig_, 5)).argmax()

print(f"Eigenvalue of the system at t = {t}: ", np.abs(np.round(valeig_, 5)))

return valstate_[index_of_max_eigenvalue].unit()

time = np.linspace(0, 4*np.pi, 1000)

oper = qt.sigmax() + qt.sigmax()

psi_lambda_list = [qt.expect(oper, psi_lambda(t, J, a = a)) for t in time]

plt.figure(figsize=(12, 8))

plt.plot(time, psi_lambda_list, "-.", label = "Numerical")

for x in [np.pi/2, 3*np.pi/2, 5*np.pi/2, 7*np.pi/2]:

plt.axvline(x, color='r', linestyle='--', label=f'x={x}')

plt.axvline(x, color='r', linestyle='--', label=f'x={x}')

plt.xlabel("Time")

plt.ylabel("<Operator>")

plt.grid(True)

plt.legend()

plt.show()

[ADITYA DEV]

Just another comment. I know there is a problem with the numeric code, but I'm not sure where the problem is. If it's either in .ptace() or ,proj() or the eigensolver etc. 

[Éric Giguère]

Removing the tidyup() in psi_lambda help. 

This improve speed but increase numerical error.

glob_ket has a lot of small values and tidyup just round them to zeros.

You could also change qutip.settings.auto_tidyup_atol to an appropriate value for your problem.

Adding a `psi_dm = (psi_dm + psi_dm.dag())/2` before the eigenvalues call also help.

There are something numerical error of the order of 1e-8 that makes psi_dm looks non-hermitian.