CoCalc Public FilesQuantum Dynamics Project (1).py
Authors: Supraja Chittari, Cullen Walsh
Views : 30
Compute Environment: Ubuntu 18.04 (Stable)
1#!/usr/bin/env python
2# coding: utf-8
3
4# In[53]:
5
6
7from sympy import *
8from numpy import sqrt
9from math import pi
10from qiskit import *
11from qiskit.visualization import plot_bloch_multivector
12
13
14# ## Two Qubit Systems
15
16# When considering sytems with two or more parties, a new mechanical feature arises - entangelement.
17
18# Seperability - if a state $\rho$ can be written as the convex combo of product states:
19
20# $$\rho = \sum p* \rho (A) \otimes \rho (B)$$
21
22#  then $\rho$  is called seperable.
23
24# Entanglement is if a state $\rho$ is not seperable.
25
26# 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.
27
28# If we consider a pure state, we can express it in the complete basis constructed from maximally entangled bell states:
29
30# $$\left| e1 \right\rangle = {1\over\sqrt{2}}(\left| 00 \right\rangle + \left| 11 \right\rangle )$$
31
32# $$\left| e2 \right\rangle = {i\over\sqrt{2}}(\left| 00 \right\rangle - \left| 11 \right\rangle )$$
33
34# $$\left| e3 \right\rangle = {i\over\sqrt{2}}(\left| 01 \right\rangle + \left| 10 \right\rangle )$$
35
36# $$\left| e4 \right\rangle = {1\over\sqrt{2}}(\left| 01 \right\rangle - \left| 10 \right\rangle )$$
37
38# Thus, we can write
39
40# $$\left| \psi \right\rangle = \sum \alpha \left| e \right\rangle$$
41
42# The concurrence is defined as
43
44# $$C(\left| \psi \right\rangle) = \sum \alpha^{2}$$
45
46# This is zero for a seperable state and 1 for a maximally entangled state.
47
48# In a mixed state, we can perform a spin-flip operation like
49
50# $$\rho ' = (\sigma 2 \otimes \sigma2) \rho^{*} (\sigma 2 \otimes \sigma 2)$$
51
52# 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:
53
54# $$C(\rho) = max(0, \lambda 1 - \lambda 2 - \lambda 3 - \lambda 4)$$
55
56# ### Bell States:
57
58# All pure and maximally entangled states are called Bell states. For a 2 qubit system, a set of four Bell states
59# build a basis for all states. This basis is:
60#     $$\left| 00 \right\rangle , \left| 01 \right\rangle , \left| 10 \right\rangle , \left| 11 \right\rangle$$
61
62# These can be expressed in vector form as:
63# $$\left| e1 \right\rangle = 64# \begin{bmatrix} 65# 1 \\ 66# 0 \\ 67# 0 \\ 68# 0 \\ 69# \end{bmatrix} 70#$$
71
72# $$\left| e2 \right\rangle = 73# \begin{bmatrix} 74# 0 \\ 75# 1 \\ 76# 0 \\ 77# 0 \\ 78# \end{bmatrix}$$
79
80# $$\left| e3 \right\rangle = 81# \begin{bmatrix} 82# 0 \\ 83# 0 \\ 84# 1 \\ 85# 0 \\ 86# \end{bmatrix}$$
87
88# $$\left| e4 \right\rangle = 89# \begin{bmatrix} 90# 0 \\ 91# 0 \\ 92# 0 \\ 93# 1 \\ 94# \end{bmatrix}$$
95
96# The bell state we will use is
97# $$\left| \psi \right\rangle = {i\over\sqrt{2}}(\left| 01 \right\rangle + \left| 10 \right\rangle )$$
98
99# The corresponding density matrix is
100#
101# $$\left| \psi \right\rangle \left\langle \psi \right| = {1 \over 2} 102# \begin{bmatrix} 103# 0 & 0 & 0 & 0\\ 104# 0 & 1 & -1 & 0\\ 105# 0 & -1 & 1 & 0\\ 106# 0 & 0 & 0 & 0 107# \end{bmatrix} 108# 109#$$
110
111# ### The Geometry of Qubits
112
113# We will talk about the simplest case - a qubit with two degrees of freedom, spin up and down (1/2 & -1/2).
114
115# The qubit is given by:
116#     $$\left| \psi \right\rangle = \alpha \left| 0 \right\rangle + \beta \left| 1 \right\rangle$$
117
118# Using the normalization condition, we can reparameretize the qubit state vector to
119#  $$\left| \psi \right\rangle = \cos \theta /2 \left| 0 \right\rangle + \sin \theta / 2 e^{i \phi} \left| 1 \right\rangle$$
120
121# 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.
122
123# In[22]:
124
125
126psi = [0,1]
127plot_bloch_multivector(psi)
128
129
130# ## Exercise
131# ### Decoherence of a Two Qubit System
132
133# 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.
134
135# $$\frac{d\rho}{dt} = -i[H, \rho] - \mathcal{D}(\rho)$$
136
137# 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
138
139# $\mathcal{D}(\rho) = \frac{1}{2} \sum_{n=k}^{d} \lambda_k (P_k \rho + \rho P_k -2P_k \rho P_k)$
140
141# Here the $P_k$ project onto the different k's of the eigenbasis $\left| e_k \right\rangle$ as below:
142
143# $$P_k = \left| e_k \right\rangle \left\langle e_k\right|$$
144
145# Taking all this into account, and that the projeciton operators must sum to 1, $\mathcal{D}(\rho)$ becomes
146
147# $$\mathcal{D}(\rho) = \lambda (\rho - \sum_{n=k}^{d} P_k \rho P_k)$$
148
149# 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:
150
151# $$\frac{d\rho}{dt} = -\mathcal{D}(\rho)$$
152
153# 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
154
155# $$\begin{pmatrix} 156# \rho_{11} & 0 & 0 & 0 \\ 157# 0 & \rho_{22} & 0 & 0 \\ 158# 0 & 0 & \rho_{33} & 0 \\ 159# 0 & 0 & 0 & \rho_{44} 160# \end{pmatrix}$$
161
162# In this case, differential equations decouple and we are left with the solutions:
163
164# $$\frac{d\rho_{ij}}{dt} = -\lambda, 165# \frac{d\rho_{ii}}{dt} = 0$$
166
167# In the general rotated basis, we find a different story. Consider the following rotation operator:
168
169# $$\left| \alpha \right \rangle = 170# \begin{pmatrix} 171# \cos\frac{\theta}{2} & \sin\frac{\theta}{2} \\ 172# -\sin\frac{\theta}{2} & \cos\frac{\theta}{2} 173# \end{pmatrix}$$
174
175# Which operates on either of the following:
176
177# $$\begin{pmatrix} 178# \left| e_1 \right\rangle \\ 179# \left| e_3 \right\rangle 180# \end{pmatrix} 181# \begin{pmatrix} 182# \left| e_2 \right\rangle \\ 183# \left| e_4 \right\rangle 184# \end{pmatrix}$$
185
186# 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):
187
188# $$\begin{pmatrix} 189# (\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 \\ 190# 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})) \\ 191# \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 \\ 192# 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}) 193# \end{pmatrix}$$
194
195# 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.
196
197# In[25]:
198
199
200import math
201import numpy as np
202from numpy import linalg as LA
203#import scipy.integrate as integrate
204import matplotlib.pyplot as plt
205
206rho = [[1,1,1,1], [1,1,1,1], [1,1,2,1], [1,1,1,1]] # trial rho
207l = 1 # lambda
208x = (math.pi)/2 # angle
209
210change_rho = np.zeros((4,4))
211
212def densityMatrix(v):
213    return np.outer(v, v)
214
215def DecoherenceA(change_rho, rho, l):
216    for i in range(4):
217        for j in range(4):
218            if i == j:
219                change_rho[i][j] = 0
220            else:
221                change_rho[i][j] = -l*rho[i][j]
222
223    return(change_rho)
224
225def oneRotated(change_rho, rho, l):
226
227    for i in range(4):
228        for j in range(4):
229            if i == j:
230                if i < 2:
231                    change_rho[i][j] = ((math.cos(x/2)**4)+(math.sin(x/2)**4))*(rho[i][i]+rho[i+2][i+2])
232                else:
233                    change_rho[i][j] = ((math.cos(x/2)**4)+(math.sin(x/2)**4))*(rho[i][i]+rho[i-2][i-2])
234            elif j == i+2:
235                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]))
236            elif i == j+2:
237                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]))
238            else:
239                change_rho[i][j] = 0
240
241    new_rho = -l*np.subtract(rho, change_rho)
242    return(new_rho)
243
244def PlotRho(f, rho, change_rho, i, j):
245    drho = np.zeros((4,4))
246    rho_list = []
247    time = []
248
249    steps = np.arange(0, 50, 1)
250
251    for n in steps:
252        if n == 0:
253            rho_list.append(rho[i][j])
254            time.append(0)
255
256        drho = f(change_rho, rho, 1)
257        rho = drho*0.1 + rho
258        t = n*0.1
259
260        rho_list.append(rho[i][j])
261        time.append(t)
262
263    plt.title('Density matrix')
264    plt.xlabel('time / a.u.')
265    plt.ylabel('Density matrix terms')
266    plt.plot(time, rho_list, 'b')
267
268def purity(rho):
269    return np.trace(np.matmul(rho,rho))
270
271def concurrence(rho):
272    y_mat = np.array([[0, -1j], [1j, 0]], dtype=np.complex)
273    y = np.tensordot(y_mat, y_mat)
274    rho_tilda = y*rho*y
275    R2 = np.sqrt(rho)*rho_tilda*np.sqrt(rho)
276    R = np.sqrt(R2)
277
278    evals = LA.eigvals(R)
279    evals_asc = np.sort(evals)
280    evals_desc = evals_asc[::-1]
281    #print(evals_desc)
282
283    c = (evals_desc[0] - evals_desc[1] - evals_desc[2] - evals_desc[3])/2
284
285    if c > 0:
286        return c
287    else:
288        return 0
289
290def PlotPurity(rho, f, purity):
291    drho = np.zeros((4,4))
292    rho2_list = []
293    time = []
294
295    steps = np.arange(0, 50, 1)
296
297    for n in steps:
298        if n == 0:
299            rho2_list.append(purity(rho))
300            time.append(0)
301
302        drho = f(change_rho, rho, 1)
303        rho = drho*0.1 + rho
304        t = n*0.1
305
306        rho2_list.append(purity(rho))
307        time.append(t)
308
309    plt.title('Density matrix')
310    plt.xlabel('time / a.u.')
311    plt.ylabel('Density matrix terms')
312    plt.plot(time, rho2_list, 'g')
313
314def PlotConcurrence(rho, f, concurrence):
315    drho = np.zeros((4,4))
316    rho2_list = []
317    time = []
318
319    steps = np.arange(0, 50, 1)
320
321    for n in steps:
322        if n == 0:
323            rho2_list.append(concurrence(rho))
324            time.append(0)
325
326        drho = f(change_rho, rho, 1)
327        rho = drho*0.1 + rho
328        #print(rho)
329        t = n*0.1
330
331        rho2_list.append(concurrence(rho))
332        time.append(t)
333
334    plt.title('Density matrix')
335    plt.xlabel('time / a.u.')
336    plt.ylabel('Density matrix terms')
337    plt.plot(time, rho2_list, 'b')
338    #print(rho2_list)
339
340singlet =np.array([0,1/math.sqrt(2),-1/math.sqrt(2),0], dtype=np.complex)
341singlet_densmat = densityMatrix(singlet)
342PlotConcurrence(singlet_densmat, DecoherenceA, concurrence)
343PlotConcurrence(singlet_densmat, oneRotated, concurrence)
344
345
346# ## Questions
347
348# 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?
349# 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:
350
351# In[ ]:
352
353
354backend = Aer.get_backend('statevector_simulator')
355final_state = execute(qc,backend).result().get_statevector()
356print(final_state)
357plot_bloch_multivector(final_state)
358
359
360# Simulate the rotations that you applied numerically on your quantum circuit.
361
362# 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.
363
364# 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?
365
366# $$P = Tr[\rho^{2}], C = max[0, \lambda_1 - \lambda_2 - \lambda_3 - \lambda_4]$$
367
368# 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.
369
370# In[ ]:
371
372
374
375
376# In[55]:
377
378
379qc = QuantumCircuit(1)
380qc.h(0)
381qc.rz(pi/4, 0)
382# See the circuit:
383qc.draw()
384
385
386# In[56]:
387
388
389# Let's see the result
390backend = Aer.get_backend('statevector_simulator')
391final_state = execute(qc,backend).result().get_statevector()
392print(final_state)
393plot_bloch_multivector(final_state)
394
395