Reproduce: Stabalizing Feedback Controls for quantum systems¶

Robert Johansson ([email protected]). Following Mirrahimi and Handel (arxiv quanth-ph/0510066 2005)

In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from qutip import *


Section 5: Two-qubit systems¶

In Sec. Mirrahimi and Handel, the authors give a feedback control law that stabalizes two interesting entangled qubit states. Here we will implement this feedback law and demonstrate it numerically.

In Eq. 5.1 in Mirrahimi and Handel, the quantum filtering equation for a two-qubit system is given by

$$d\rho_t(t) = - i u_1(t) [\sigma_y^{(1)}, \rho_t(t)] dt - i u_2(t) [\sigma_y^{(2)}, \rho_t(t)] dt -\frac{1}{2} [F_z, [F_z, \rho_t(t)]]dt + \sqrt{\eta} \left\{F_z\rho_t(t) + \rho_t(t)F_z - 2{\rm Tr}[F_z \rho_t(t)]\rho_t(t)\right\} dW_t$$

or $\mathcal{D}$ and $\mathcal{H}$ we can write:

$$d\rho_t(t) = - i u_1(t) [\sigma_y^{(1)}, \rho_t(t)] dt - i u_2(t) [\sigma_y^{(2)}, \rho_t(t)] dt + \mathcal{D}[F_z] \rho_t(t) dt + \sqrt{\eta} \mathcal{H}[F_z] \rho_t(t) dW_t$$

Implementation in QuTiP¶

In QuTiP we write the stochastic master equation on the form:

$$d\rho(t) = -i[H, \rho(t)] dt + \mathcal{D}[B] \rho(t) dt + D_{1}[A]\rho(t) dt + D_{2}[A]\rho(t) dW$$

With $H = 0$, $B = F_z$, and $A = F_z$, we can identify:

$$D_{1}[A]\rho_t(t) = -i u_1(t) [\sigma_y^{(1)}, \rho_t(t)] -i u_2(t) [\sigma_y^{(2)}, \rho_t(t)] + \mathcal{D}[A] \rho_t(t)$$$$D_{2}[A]\rho_t(t) = \sqrt{\eta} \mathcal{H}[A] \rho_t(t)$$

Setup operators and states¶

In [2]:
sx1, sy1, sz1 = [tensor(op, identity(2)) for op in (sigmax(), sigmay(), sigmaz())]
sx2, sy2, sz2 = [tensor(identity(2), op) for op in (sigmax(), sigmay(), sigmaz())]

In [3]:
Fx, Fy, Fz = sx1 + sx2, sy1 + sy2, sz1 + sz2

In [4]:
psi_ud = tensor(basis(2, 1), basis(2, 0))
psi_du = tensor(basis(2, 0), basis(2, 1))

psi_s = (psi_ud + psi_du).unit()
psi_a = (psi_ud - psi_du).unit()


The purpose of Sec. 5 in Mirrahimi and Handel is to give a feedback control law that stabalize the following two entangled qubit states:

In [5]:
rho_s = psi_s * psi_s.dag()
rho_a = psi_a * psi_a.dag()

In [6]:
# problem parameters
eta = 0.4
gamma = 0.15

In [7]:
# simulation paramters
T = 25
nsteps = 1000
nsubsteps = 100
times = np.linspace(0, T, nsteps)
ntraj = 1

In [8]:
# a random initial state
rho0 = tensor(ket2dm(rand_ket(2, density=1)), ket2dm(rand_ket(2, density=1)))

In [9]:
# expectation value operators
e_ops = [sx1, sy1, sz1, sx2, sy2, sz2]

In [10]:
# Hamiltonian, deterministic and stochastic collapse operators
H = 0 * Fz
c_ops = []
sc_ops = [Fz]


To implement the control laws for the two target states in QuTiP we need to defined the d1 and d2 functions for the stochastic master equation:

$D_1$¶

$$D_{1}[F_z]\rho_t(t) = -i u_1(t) [\sigma_y^{(1)}, \rho_t(t)] -i u(t) [\sigma_y^{(2)}, \rho_t(t)] + \mathcal{D}[F_z] \rho_t(t) = u_1(t) L[\sigma_y^{(1)}] \rho_t(t) + u_2(t) L[\sigma_y^{(2)}] \rho_t(t) + \mathcal{D}[F_z] \rho_t(t)$$
In [11]:
from qutip.expect import expect_rho_vec

In [12]:
# precompute required matrices and vectors used in u_a_t and u_s_t
L_sy1_data = liouvillian(sy1, []).data
L_sy2_data = liouvillian(sy2, []).data

rho_s_vec = mat2vec(rho_s.full()).ravel()
rho_s_mat = vec2mat(rho_s_vec)

rho_a_vec = mat2vec(rho_a.full()).ravel()
rho_a_mat = vec2mat(rho_a_vec)

# define functions for the control laws. Theorem 5 in Mirrahimi and Handel

_u_a_prev_region = 0
def u_a_t(A, rho_vec):
"""Feedback control from Mirrahimi and Handel:
Compute and return the u_1 and u_2 functions for
the case when the target state is rho_a.
"""
global _u1_prev_region

tr_rr = np.diag(np.dot(rho_a_mat, vec2mat(rho_vec))).sum().real

if tr_rr < gamma / 2:
_u1_a_prev_region = 0
u_1 = 1
u_2 = 0

elif tr_rr >= gamma or _u_a_prev_region == 1:
_u1_a_prev_region = 1

x =  - L_sy1_data * rho_vec
y = np.diag(np.dot(vec2mat(x), rho_a_mat)).sum().real
u_1 = 1 - y

x = - L_sy2_data * rho_vec
y = np.diag(np.dot(vec2mat(x), rho_a_mat)).sum().real
u_2 = 1 - y

else:
u_1 = 1
u_2 = 0

return u_1, u_2

_u_s_prev_region = 0
def u_s_t(A, rho_vec):
"""Feedback control from Mirrahimi and Handel.
Compute and return the u_1 and u_2 functions for
the case when the target state is rho_s.
"""
global _u1_prev_region

tr_rr = np.diag(np.dot(rho_s_mat, vec2mat(rho_vec))).sum().real

if tr_rr < gamma / 2:
_u1_s_prev_region = 0
u_1 = 1
u_2 = 0

elif tr_rr >= gamma or _u_s_prev_region == 1:
_u1_s_prev_region = 1

x =  - L_sy1_data * rho_vec
y = np.diag(np.dot(vec2mat(x), rho_s_mat)).sum().real
u_1 = 1 - y

x = - L_sy2_data * rho_vec
y = np.diag(np.dot(vec2mat(x), rho_s_mat)).sum().real
u_2 = - 1 - y

else:
u_1 = 1
u_2 = 0

return u_1, u_2

def d1_s(t, rho_vec, A, args):
# rho_s as target
u_1, u_2 = u_s_t(A, rho_vec)

return (u_1 * L_sy1_data * rho_vec +
u_2 * L_sy2_data * rho_vec +
A[7] * rho_vec)

def d1_a(t, rho_vec, A, args):
# rho_a as target
u_1, u_2 = u_a_t(A, rho_vec)

return (u_1 * L_sy1_data * rho_vec +
u_2 * L_sy2_data * rho_vec +
A[7] * rho_vec)


$D_2$¶

$$D_{2}[F_z]\rho_t(t) = \sqrt{\eta} \mathcal{H}[F_z] \rho_t(t) = \sqrt{\eta}(F_z\rho_t(t) + \rho_t(t) F_z - \mathrm{Tr}[F_z\rho_t(t) + \rho_t(t) F_z] \rho_t(t)) \rightarrow \sqrt{\eta} \left((A_L + A_R^\dagger)\rho_t - \mathrm{Tr}[(A_L + A_R^\dagger)\rho_t] \rho_t\right)$$
In [13]:
def d2(t, rho_vec, A, args):
n_sum = A[0] + A[3]
e1 = expect_rho_vec(n_sum, rho_vec, False)
return [np.sqrt(eta) * (n_sum * rho_vec - e1 * rho_vec)]


Simulation¶

In [14]:
# solve using the d2 implementation that stabilized rho_s
result_s = smesolve(H, rho0, times, [], sc_ops, e_ops,
ntraj=ntraj, nsubsteps=nsubsteps, d1=d1_s, d2=d2,
distribution='normal',
options=Options(store_states=True, average_states=True))

Total run time: 105.36s

In [15]:
# solve using the d2 implementation that stabilized rho_a
result_a = smesolve(H, rho0, times, [], sc_ops, e_ops,
ntraj=ntraj, nsubsteps=nsubsteps, d1=d1_a, d2=d2,
distribution='normal',
options=Options(store_states=True, average_states=True))

Total run time: 108.19s


Functions for visualizing the results¶

In [16]:
def visualize_dynamics(result, rho_f):

fig, axes = plt.subplots(1, 3, figsize=(16,5))

## Expectation values of the collective degrees of freedom
axes[0].plot(times, result.expect[0] + result.expect[3], 'r', lw=2, label=r'$\langle F_x\rangle$')
axes[0].plot(times, result.expect[1] + result.expect[4], 'b', lw=2, label=r'$\langle F_y\rangle$')
axes[0].plot(times, result.expect[2] + result.expect[5], 'k', lw=2, label=r'$\langle F_z\rangle$')

axes[0].plot(times, np.ones_like(times) * expect(Fx, rho_f), 'r--', lw=2)
axes[0].plot(times, np.ones_like(times) * expect(Fy, rho_f), 'b--', lw=2)
axes[0].plot(times, np.ones_like(times) * expect(Fz, rho_f), 'k--', lw=2)

axes[0].set_ylim(-2.1, 2.1)
axes[0].set_xlim(0, times.max())
axes[0].set_xlabel('time', fontsize=12)
axes[0].legend(loc=4);

## Expectation values of qubit 1
axes[1].plot(times, result.expect[0], 'r', lw=2, label=r'$\langle\sigma_x^{(1)}\rangle$')
axes[1].plot(times, result.expect[1], 'b', lw=2, label=r'$\langle\sigma_y^{(1)}\rangle$')
axes[1].plot(times, result.expect[2], 'k', lw=2, label=r'$\langle\sigma_z^{(1)}\rangle$')

axes[1].plot(times, np.ones_like(times) * expect(sx1, rho_f), 'r--', lw=2)
axes[1].plot(times, np.ones_like(times) * expect(sy1, rho_f), 'b--', lw=2)
axes[1].plot(times, np.ones_like(times) * expect(sz1, rho_f), 'k--', lw=2)

axes[1].set_ylim(-1.1, 1.1)
axes[1].set_xlim(0, times.max())
axes[1].set_xlabel('time', fontsize=12)
axes[1].legend(loc=4);

## Expectation values of qubit 2
axes[2].plot(times, result.expect[3], 'r', lw=2, label=r'$\langle\sigma_x^{(2)}\rangle$')
axes[2].plot(times, result.expect[4], 'b', lw=2, label=r'$\langle\sigma_y^{(2)}\rangle$')
axes[2].plot(times, result.expect[5], 'k', lw=2, label=r'$\langle\sigma_z^{(2)}\rangle$')

axes[2].plot(times, np.ones_like(times) * expect(sx2, rho_f), 'r--', lw=2)
axes[2].plot(times, np.ones_like(times) * expect(sy2, rho_f), 'b--', lw=2)
axes[2].plot(times, np.ones_like(times) * expect(sz2, rho_f), 'k--', lw=2)

axes[2].set_ylim(-1.1, 1.1)
axes[2].set_xlim(0, times.max())
axes[2].set_xlabel('time', fontsize=12)
axes[2].legend(loc=4)

return fig, axes

In [17]:
def visualize_dynamics_blochsphere(result, rho_f):
fig, axes = plt.subplots(1, 2, figsize=(16, 8), subplot_kw={'projection': '3d'})

b = Bloch(fig=fig, axes=axes[0])
b.make_sphere()

b = Bloch(fig=fig, axes=axes[1])
b.make_sphere()

return fig, axes

In [18]:
def visualize_fidelity(result, rho_target):

fig, ax = plt.subplots(1, 1, figsize=(8,5))

f = [fidelity(rho_target, rho) for rho in result.states]
c = [concurrence(rho) for rho in result.states]
ax.plot(result.times, f, 'r', lw=2, label=r'fidelity')
ax.plot(result.times, c, 'b', lw=2, label=r'concurrence')
ax.set_ylim(0, 1)
ax.set_xlim(0, times.max())
ax.set_xlabel('time', fontsize=12)
ax.legend(loc=4)

return fig, ax


Stabilizing $\rho_s$¶

In [19]:
visualize_dynamics(result_s, rho_s);

In [20]:
rho = result_s.states[-1].tidyup(atol=1e-3)
rho

Out[20]:
Quantum object: dims = [[2, 2], [2, 2]], shape = [4, 4], type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.500 & 0.500 & 0.0\\0.0 & 0.500 & 0.500 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\\end{array}\right)\end{equation*}
In [21]:
#target
rho_s

Out[21]:
Quantum object: dims = [[2, 2], [2, 2]], shape = [4, 4], type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.500 & 0.500 & 0.0\\0.0 & 0.500 & 0.500 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\\end{array}\right)\end{equation*}
In [22]:
# fidelity at final time in simulation
fidelity(rho_s, rho)

Out[22]:
0.9999989773045586
In [23]:
# evolution of fidelity from the random initial state
visualize_fidelity(result_a, rho_a);

In [24]:
visualize_dynamics_blochsphere(result_s, rho_s);


Stabilizing $\rho_a$¶

In [25]:
visualize_dynamics(result_a, rho_a);

In [26]:
rho = result_a.states[-1].tidyup(atol=1e-3)
rho

Out[26]:
Quantum object: dims = [[2, 2], [2, 2]], shape = [4, 4], type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.500 & -0.500 & 0.0\\0.0 & -0.500 & 0.500 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\\end{array}\right)\end{equation*}
In [27]:
rho_a

Out[27]:
Quantum object: dims = [[2, 2], [2, 2]], shape = [4, 4], type = oper, isherm = True\begin{equation*}\left(\begin{array}{*{11}c}0.0 & 0.0 & 0.0 & 0.0\\0.0 & 0.500 & -0.500 & 0.0\\0.0 & -0.500 & 0.500 & 0.0\\0.0 & 0.0 & 0.0 & 0.0\\\end{array}\right)\end{equation*}
In [28]:
# fidelity at final time in simulation
fidelity(rho_a, rho)

Out[28]:
1.0000000234736035
In [29]:
# evolution of fidelity from the random initial state
visualize_fidelity(result_a, rho_a);

In [30]:
visualize_dynamics_blochsphere(result_a, rho_a);


Versions¶

In [31]:
from qutip.ipynbtools import version_table

version_table()

Out[31]:
SoftwareVersion
QuTiP3.1.0.dev-2cc9b57
OSposix [linux]
IPython2.3.0
SciPy0.14.0
Cython0.20.2
Python3.4.0 (default, Apr 11 2014, 13:05:11) [GCC 4.8.2]
Numpy1.8.2
matplotlib1.4.0
Fri Nov 07 14:15:19 2014 JST