I've only recently started using density matrices in my work but I am confused with the following code that I have whether I am getting the right matrix:
def Hamiltonian(alpha,h):
Sx = np.array([[0,1],[1,0]])
Sy = np.array([[0,-1j],[1j,0]])
Sz = np.array([[1,0],[0,-1]])
I = np.array([[1,0],[0,1]])
H = (alpha*np.kron(np.kron(Sx,Sx),I))
H =+ (alpha*np.kron(np.kron(Sy,Sy),I))
H =+ (alpha*np.kron(np.kron(I,Sx),Sx))
H =+ (alpha*np.kron(np.kron(I,Sy),Sy))
H =+ (h*np.kron(np.kron(I,Sz),I))
return H
So the above gives me my Hamiltonian Function, where alpha is a real number and h is a magnetization parameter acting on one of my qubits.
My ultimate goal with this is to plot $\lvert U\rvert$ vs. $h$, where $U = \text{Tr}\left(H\rho\right)$. I know that the density matrix follows the form of the outer product of the wavefunction:
$\rho = \left|ψ\rangle\langle ψ\right|$
But my question is twofold:
- How do I get the density matrix for the above Hamiltonian function? Has it got something to do with the eigenvalues and eigenvectors?
I have tried the following:
H = Hamiltonian(1,0.5)
print(H)
$$\begin{bmatrix} 0.5&0&0&0&0&0&0&0 \\ 0&0.5&0&0&0&0&0&0 \\ 0&0&-0.5&0&0&0&0&0 \\ 0&0&0&-0.5&0&0&0&0 \\ 0&0&0&0&0.5&0&0&0 \\ 0&0&0&0&0&0.5&0&0 \\ 0&0&0&0&0&0&-0.5&0 \\ 0&0&0&0&0&0&0&-0.5 \end{bmatrix}$$
w, v = LA.eig(H)
print(w)
print(v)
$$\begin{bmatrix}0.5&0.5&-0.5&-0.5&0.5&0.5&-0.5&-0.5\end{bmatrix} \\ \begin{bmatrix} 1&0&0&0&0&0&0&0 \\ 0&1&0&0&0&0&0&0 \\ 0&0&1&0&0&0&0&0 \\ 0&0&0&1&0&0&0&0 \\ 0&0&0&0&1&0&0&0 \\ 0&0&0&0&0&1&0&0 \\ 0&0&0&0&0&0&1&0 \\ 0&0&0&0&0&0&0&1 \end{bmatrix}$$
ro = np.outer(w,w)
print(ro)
$$\begin{bmatrix} 0.25&0.25&-0.25&-0.25&0.25&0.25&-0.25&-0.25 \\ 0.25&0.25&-0.25&-0.25&0.25&0.25&-0.25&-0.25 \\ -0.25&-0.25&0.25&0.25&-0.25&-0.25&0.25&0.25 \\ -0.25&-0.25&0.25&0.25&-0.25&-0.25&0.25&0.25 \\ 0.25&0.25&-0.25&-0.25&0.25&0.25&-0.25&-0.25 \\ 0.25&0.25&-0.25&-0.25&0.25&0.25&-0.25&-0.25 \\ -0.25&-0.25&0.25&0.25&-0.25&-0.25&0.25&0.25 \\ -0.25&-0.25&0.25&0.25&-0.25&-0.25&0.25&0.25 \\ \end{bmatrix}$$
ro_square = ro**2
purity = np.trace(ro_square)
print(purity)
Comes out to be 0.5 which makes sense for my mixed state
U = (np.trace(H*ro))
print(U)
Always $U = 0$ for all $h$...
But for every value of $h$ I use I get $U = 0$ so something is going wrong, I feel like I'm not using the correct parameters for the density matrix in the outer product of the eigenvalues
- How do I get the Hamiltonian in terms of an operator or is it already in this form?
Any help will be appreciated and if this Why is the wrong forum to put this question up let me know and I will put it up on the normal physics forum thanks.
I also just want to clarify I can code the answer when I understand what I'm using to get the density matrix I'm not asking anyone for the code I just need an explanation in terms of the physics of it all!diagonal?
