CoCalc Public FilesQuantum Dynamics Project (1).pyOpen with one click!
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
7
from sympy import *
8
from numpy import sqrt
9
from math import pi
10
from qiskit import *
11
from 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} \begin{equation}
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
# \end{equation}
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
126
psi = [0,1]
127
plot_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
200
import math
201
import numpy as np
202
from numpy import linalg as LA
203
#import scipy.integrate as integrate
204
import matplotlib.pyplot as plt
205
206
rho = [[1,1,1,1], [1,1,1,1], [1,1,2,1], [1,1,1,1]] # trial rho
207
l = 1 # lambda
208
x = (math.pi)/2 # angle
209
210
change_rho = np.zeros((4,4))
211
212
def densityMatrix(v):
213
return np.outer(v, v)
214
215
def 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
225
def 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
244
def 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
268
def purity(rho):
269
return np.trace(np.matmul(rho,rho))
270
271
def 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
290
def 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
314
def 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
340
singlet =np.array([0,1/math.sqrt(2),-1/math.sqrt(2),0], dtype=np.complex)
341
singlet_densmat = densityMatrix(singlet)
342
PlotConcurrence(singlet_densmat, DecoherenceA, concurrence)
343
PlotConcurrence(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
354
backend = Aer.get_backend('statevector_simulator')
355
final_state = execute(qc,backend).result().get_statevector()
356
print(final_state)
357
plot_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
373
1. Answer:
374
375
376
# In[55]:
377
378
379
qc = QuantumCircuit(1)
380
qc.h(0)
381
qc.rz(pi/4, 0)
382
# See the circuit:
383
qc.draw()
384
385
386
# In[56]:
387
388
389
# Let's see the result
390
backend = Aer.get_backend('statevector_simulator')
391
final_state = execute(qc,backend).result().get_statevector()
392
print(final_state)
393
plot_bloch_multivector(final_state)
394
395