Quantum Dynamics Project (1).py

#!/usr/bin/env python
# coding: utf-8

# In[53]:


from sympy import *
from numpy import sqrt
from math import pi
from qiskit import *
from qiskit.visualization import plot_bloch_multivector


# ## Two Qubit Systems

# When considering sytems with two or more parties, a new mechanical feature arises - entangelement. 

# Seperability - if a state $ \rho $ can be written as the convex combo of product states: 

# $$ \rho = \sum p* \rho (A) \otimes \rho (B) $$

#  then $ \rho $  is called seperable.

# Entanglement is if a state $ \rho $ is not seperable. 

# A major challenge is in determining the degree of entanglement in a system. FOr a two qubit system, the concurrence is a good measure of entanglement.

# If we consider a pure state, we can express it in the complete basis constructed from maximally entangled bell states:

# $$ \left| e1 \right\rangle =  {1\over\sqrt{2}}(\left| 00 \right\rangle + \left| 11 \right\rangle )$$

# $$ \left| e2 \right\rangle =  {i\over\sqrt{2}}(\left| 00 \right\rangle - \left| 11 \right\rangle )$$

# $$ \left| e3 \right\rangle =  {i\over\sqrt{2}}(\left| 01 \right\rangle + \left| 10 \right\rangle )$$

# $$ \left| e4 \right\rangle =  {1\over\sqrt{2}}(\left| 01 \right\rangle - \left| 10 \right\rangle )$$

# Thus, we can write

# $$ \left| \psi \right\rangle = \sum \alpha \left| e \right\rangle $$

# The concurrence is defined as 

# $$ C(\left| \psi \right\rangle) = \sum \alpha^{2}$$

# This is zero for a seperable state and 1 for a maximally entangled state.

# In a mixed state, we can perform a spin-flip operation like

# $$ \rho ' = (\sigma 2 \otimes \sigma2) \rho^{*} (\sigma 2 \otimes \sigma 2) $$

# where $\sigma 2$ is a pauli matrix. Now we can define $$ R^{2} = \rho * \rho ' $$ where the lambda terms are the squra roots of the eigenvalues of the matrix R^2. The concurrence of the mixed state is then given by:

# $$ C(\rho) = max(0, \lambda 1 - \lambda 2 - \lambda 3 - \lambda 4) $$

# ### Bell States:

# All pure and maximally entangled states are called Bell states. For a 2 qubit system, a set of four Bell states
# build a basis for all states. This basis is:
#     $$ \left| 00 \right\rangle , \left| 01 \right\rangle , \left| 10 \right\rangle ,  \left| 11 \right\rangle $$

# These can be expressed in vector form as:
# $$ \left| e1 \right\rangle = 
# \begin{bmatrix}
# 1 \\
# 0 \\
# 0 \\
# 0 \\
# \end{bmatrix}
# $$

# $$ \left| e2 \right\rangle = 
# \begin{bmatrix}
# 0 \\
# 1 \\
# 0 \\
# 0 \\
# \end{bmatrix} $$

# $$ \left| e3 \right\rangle = 
# \begin{bmatrix}
# 0 \\
# 0 \\
# 1 \\
# 0 \\
# \end{bmatrix} $$

# $$ \left| e4 \right\rangle = 
# \begin{bmatrix}
# 0 \\
# 0 \\
# 0 \\
# 1 \\
# \end{bmatrix} $$

# The bell state we will use is 
# $$ \left| \psi \right\rangle =  {i\over\sqrt{2}}(\left| 01 \right\rangle + \left| 10 \right\rangle ) $$

# The corresponding density matrix is
# 
# $$ \left| \psi \right\rangle \left\langle \psi \right| = {1 \over 2} \begin{equation}
# \begin{bmatrix}
# 0 & 0 & 0 & 0\\
# 0 & 1 & -1 & 0\\
# 0 & -1 & 1 & 0\\
# 0 & 0 & 0 & 0
# \end{bmatrix}
# \end{equation} 
# $$

# ### The Geometry of Qubits

# We will talk about the simplest case - a qubit with two degrees of freedom, spin up and down (1/2 & -1/2).

# The qubit is given by:
#     $$ \left| \psi \right\rangle = \alpha \left| 0 \right\rangle + \beta \left| 1 \right\rangle $$

# Using the normalization condition, we can reparameretize the qubit state vector to
#  $$ \left| \psi \right\rangle = \cos \theta /2 \left| 0 \right\rangle + \sin \theta / 2 e^{i \phi} \left| 1 \right\rangle $$

# The angular parameters define a point on the unit sphere in 3 dimensional space. This space is known as Bloch sphere. Here is an example for the $ \left| 1 \right\rangle $ state.

# In[22]:


psi = [0,1]
plot_bloch_multivector(psi)


# ## Exercise 
# ### Decoherence of a Two Qubit System

# In order to tackle this problem, we are going to write a Markovian quantum master equation to describe how the system evolves in time. We will use a specific form of this master equation, the Linblad master equation (below), of a quantum system weakly coupled to a Markovian reservoir. Getting to this point is not necessarily straightforward so we'll not delve into that discussion.

# $$\frac{d\rho}{dt} = -i[H, \rho] - \mathcal{D}(\rho)$$

# To get to the functional form of the operator used here, we futher assume that Linblad operators (which make up $\mathcal{D}$) are Hermitian and commuting with the Hamiltonian. We can then replace the Linblad operators with projection operators and arrive at 

# $\mathcal{D}(\rho) = \frac{1}{2} \sum_{n=k}^{d} \lambda_k (P_k \rho + \rho P_k -2P_k \rho P_k)$

# Here the $P_k$ project onto the different k's of the eigenbasis $\left| e_k \right\rangle$ as below:

# $$P_k = \left| e_k \right\rangle \left\langle e_k\right|$$

# Taking all this into account, and that the projeciton operators must sum to 1, $\mathcal{D}(\rho)$ becomes

# $$\mathcal{D}(\rho) = \lambda (\rho - \sum_{n=k}^{d} P_k \rho P_k)$$

# The advantage of using this dissipator is that the positive and real parameters $\lambda_k$ are the decoherence parameters and represent the strength the the system coupled to its surroundings, which is a process that reduces the purity of a microstate. The master equation is now: 

# $$\frac{d\rho}{dt} = -\mathcal{D}(\rho)$$

# We will now try to solve the master equation with this decoherence formalism for a 2 qubit system. Before we write down any specific density matrices, we will first determine how the dissipator acts on different bases. We will call the first decoherence mode simply $E{\otimes}E$, where both qubits are written in the Hamiltonian eigebasis. Therefore, we find that the dissipator acted on a general density matrix looks like

# $$ \begin{pmatrix}
# \rho_{11} & 0 & 0 & 0 \\
# 0 & \rho_{22} & 0 & 0 \\
# 0 & 0 & \rho_{33} & 0 \\
# 0 & 0 & 0 & \rho_{44} 
# \end{pmatrix} $$

# In this case, differential equations decouple and we are left with the solutions:

# $$\frac{d\rho_{ij}}{dt} = -\lambda, 
# \frac{d\rho_{ii}}{dt} = 0 $$

# In the general rotated basis, we find a different story. Consider the following rotation operator:

# $$ \left| \alpha \right \rangle = 
# \begin{pmatrix}
# \cos\frac{\theta}{2} & \sin\frac{\theta}{2} \\
# -\sin\frac{\theta}{2} & \cos\frac{\theta}{2}
# \end{pmatrix} $$

# Which operates on either of the following:

# $$ \begin{pmatrix}
# \left| e_1 \right\rangle \\
# \left| e_3 \right\rangle
# \end{pmatrix} 
# \begin{pmatrix}
# \left| e_2 \right\rangle \\
# \left| e_4 \right\rangle
# \end{pmatrix} $$

# Now we have a situation where each of the four basis vectors, under this general rotation operator, will give four new basis vectors $\left| \beta_k \right\rangle$ which have projection operators: $P_k = \left| \beta_k \right\rangle \left\langle \beta_k\right|$. We call this new space $GR \otimes E$. Now we play the same game again and assemble the dissipator operator. Below is the right side (summation part):

# $$ \begin{pmatrix}
# (\cos^{4}\frac{\theta}{2}+\sin^{2}\frac{\theta}{4})(\rho_{11}+\rho_{33}) & 0 & \sin\frac{\theta}{2}\cos\frac{\theta}{2}((\cos^{2}\frac{\theta}{2}-\sin^{2}\frac{\theta}{2})(\rho_{11}-\rho_{33})-2\cos\frac{\theta}{2}\cos\frac{\theta}{2}(\rho_{13}+\rho_{31})) & 0 \\
# 0 & (\cos^{4}\frac{\theta}{2}+\sin^{4}\frac{\theta}{4})(\rho_{22}+\rho_{44}) & 0 & \sin\frac{\theta}{2}\cos\frac{\theta}{2}((\cos^{2}\frac{\theta}{2}-\sin^{2}\frac{\theta}{2})(\rho_{22}-\rho_{44})-2\cos\frac{\theta}{2}\cos\frac{\theta}{2}(\rho_{24}+\rho_{42})) \\
# \sin\frac{\theta}{2}\cos\frac{\theta}{2}((\cos^{2}\frac{\theta}{2}-\sin^{2}\frac{\theta}{2})(\rho_{11}-\rho_{33})-2\cos\frac{\theta}{2}\cos\frac{\theta}{2}(\rho_{13}+\rho_{31})) & 0 & (\cos^{4}\frac{\theta}{2}+\sin^{4}\frac{\theta}{4})(\rho_{11}+\rho_{33}) & 0 \\
# 0 & \sin\frac{\theta}{2}\cos\frac{\theta}{2}((\cos^{2}\frac{\theta}{2}-\sin^{2}\frac{\theta}{2})(\rho_{22}-\rho_{44})-2\cos\frac{\theta}{2}\cos\frac{\theta}{2}(\rho_{24}+\rho_{42})) & 0 & (\cos^{4}\frac{\theta}{2}+\sin^{4}\frac{\theta}{4})(\rho_{22}+\rho_{44})
# \end{pmatrix} $$

# This matrix now gets subtracted from $\rho$, which is messy algebra. Anyways, we will compute this problem numerically. So let's now turn our attention to the python script.

# In[25]:


import math
import numpy as np
from numpy import linalg as LA
#import scipy.integrate as integrate
import matplotlib.pyplot as plt

rho = [[1,1,1,1], [1,1,1,1], [1,1,2,1], [1,1,1,1]] # trial rho
l = 1 # lambda
x = (math.pi)/2 # angle

change_rho = np.zeros((4,4))

def densityMatrix(v):
    return np.outer(v, v)
    
def DecoherenceA(change_rho, rho, l):    
    for i in range(4):
        for j in range(4):
            if i == j:
                change_rho[i][j] = 0
            else:
                change_rho[i][j] = -l*rho[i][j]
    
    return(change_rho)

def oneRotated(change_rho, rho, l):
    
    for i in range(4):
        for j in range(4):
            if i == j:
                if i < 2:
                    change_rho[i][j] = ((math.cos(x/2)**4)+(math.sin(x/2)**4))*(rho[i][i]+rho[i+2][i+2])  
                else:
                    change_rho[i][j] = ((math.cos(x/2)**4)+(math.sin(x/2)**4))*(rho[i][i]+rho[i-2][i-2])
            elif j == i+2:
                change_rho[i][j] = (math.sin(x/2)*math.cos(x/2))*(((math.cos(x/2)**4)-(math.sin(x/2)**4))*(rho[i][i]-rho[i+2][i+2])+(2*math.cos(x/2)*math.sin(x/2))*(rho[i][i+2] + rho[i+2][i]))
            elif i == j+2:
                change_rho[i][j] = (math.sin(x/2)*math.cos(x/2))*(((math.cos(x/2)**4)-(math.sin(x/2)**4))*(rho[i][i]-rho[i-2][i-2])+(2*math.cos(x/2)*math.sin(x/2))*(rho[i][i-2] + rho[i-2][i]))
            else:
                change_rho[i][j] = 0
                
    new_rho = -l*np.subtract(rho, change_rho)
    return(new_rho)

def PlotRho(f, rho, change_rho, i, j):
    drho = np.zeros((4,4))
    rho_list = []
    time = []
    
    steps = np.arange(0, 50, 1)

    for n in steps:
        if n == 0:
            rho_list.append(rho[i][j])
            time.append(0)
            
        drho = f(change_rho, rho, 1)
        rho = drho*0.1 + rho
        t = n*0.1
        
        rho_list.append(rho[i][j])
        time.append(t)
    
    plt.title('Density matrix')
    plt.xlabel('time / a.u.')
    plt.ylabel('Density matrix terms')
    plt.plot(time, rho_list, 'b')

def purity(rho):
    return np.trace(np.matmul(rho,rho))

def concurrence(rho):
    y_mat = np.array([[0, -1j], [1j, 0]], dtype=np.complex)
    y = np.tensordot(y_mat, y_mat)
    rho_tilda = y*rho*y
    R2 = np.sqrt(rho)*rho_tilda*np.sqrt(rho)
    R = np.sqrt(R2)
    
    evals = LA.eigvals(R)
    evals_asc = np.sort(evals)
    evals_desc = evals_asc[::-1]
    #print(evals_desc)
    
    c = (evals_desc[0] - evals_desc[1] - evals_desc[2] - evals_desc[3])/2
    
    if c > 0:
        return c
    else:
        return 0

def PlotPurity(rho, f, purity):
    drho = np.zeros((4,4))
    rho2_list = []
    time = []
    
    steps = np.arange(0, 50, 1)

    for n in steps:
        if n == 0:
            rho2_list.append(purity(rho))
            time.append(0)
            
        drho = f(change_rho, rho, 1)
        rho = drho*0.1 + rho
        t = n*0.1
        
        rho2_list.append(purity(rho))
        time.append(t)
    
    plt.title('Density matrix')
    plt.xlabel('time / a.u.')
    plt.ylabel('Density matrix terms')
    plt.plot(time, rho2_list, 'g')

def PlotConcurrence(rho, f, concurrence):
    drho = np.zeros((4,4))
    rho2_list = []
    time = []
    
    steps = np.arange(0, 50, 1)

    for n in steps:
        if n == 0:
            rho2_list.append(concurrence(rho))
            time.append(0)
            
        drho = f(change_rho, rho, 1)
        rho = drho*0.1 + rho
        #print(rho)
        t = n*0.1
        
        rho2_list.append(concurrence(rho))
        time.append(t)
    
    plt.title('Density matrix')
    plt.xlabel('time / a.u.')
    plt.ylabel('Density matrix terms')
    plt.plot(time, rho2_list, 'b')
    #print(rho2_list)

singlet =np.array([0,1/math.sqrt(2),-1/math.sqrt(2),0], dtype=np.complex)
singlet_densmat = densityMatrix(singlet)
PlotConcurrence(singlet_densmat, DecoherenceA, concurrence)
PlotConcurrence(singlet_densmat, oneRotated, concurrence)


# ## Questions

# 1. Using the qubit state vector equation, rotate the (0,1) vector by three different angles (thetas) and plot them on a bloch sphere. What do you observe is happening in geometric space? What value will give you the pure (1,0) state? 
# Now, we are going to visual this using a quantum circuit. Create a quantum circuit using qc = QuantumCircuit(1). We will then apply a Hadmard gate using qc.h(0). This will create a superposition state of spin up and down. Next we will apply a rotation gate using qc.rz(angle1,angle2). You can then draw the gate using qc.draw(). To visualize your final rotated superposition state in vector and geomertric form, we will use a quantum simulator: 

# In[ ]:


backend = Aer.get_backend('statevector_simulator')
final_state = execute(qc,backend).result().get_statevector()
print(final_state)
plot_bloch_multivector(final_state)


# Simulate the rotations that you applied numerically on your quantum circuit. 

# 2. Now let's plot the $GR \otimes E$ case for the singlet initial state. First plot a test case of $\theta = \frac{\theta}{2}$ and show the time dependence of the $\rho_{11}$, $\rho_{12}$, and $\rho_{13}$ terms. What happens? Then, choosing an arbitrary time, modify the code to plot the $\theta$ dependence of each of these terms. 

# 3. A useful way for us to describe these density matrices is by computing the purity and concurrence of the matrix, as defined below. A pure and maximally entangled state will have P=1 and C=1. What happens to these values for a singlet initial state that evolves by decoherence mode A? What about the $GR \otimes E$ case?

# $$P = Tr[\rho^{2}], C = max[0, \lambda_1 - \lambda_2 - \lambda_3 - \lambda_4]$$

# Where $\lambda$ are the eigenvalues of the matrix $\sqrt{\rho \widetilde{\rho} \rho}$. Here $\rho$ is the density matrix and $\widetilde{\rho}$ is the matrix $(\sigma_y \otimes \sigma_y)\rho^{*}(\sigma_y \otimes \sigma_y)$. We again compute these two numerically at every chosen time point.

# In[ ]:


1. Answer:


# In[55]:


qc = QuantumCircuit(1)
qc.h(0)
qc.rz(pi/4, 0)
# See the circuit:
qc.draw()


# In[56]:


# Let's see the result
backend = Aer.get_backend('statevector_simulator')
final_state = execute(qc,backend).result().get_statevector()
print(final_state)
plot_bloch_multivector(final_state)