Ondrej Cernotik
23 February 2015
Using a cavity mode dispersively coupled to two qubits and measuring the output field, it is possible to entangle the two qubits. The stochastic master equation for the qubits, in the presence of thermal noise for the cavity field, takes the form $$ \mathrm{d}\rho = \frac{2\chi^2}{\kappa}(2\bar{n}+1)\mathcal{D}[\sigma_z^1+\sigma_z^2]\rho\mathrm{d}t+\sqrt{\frac{2\chi^2}{\kappa(2\bar{n}+1)}}\mathcal{H}[\sigma_z^1+\sigma_z^2]\rho\mathrm{d}W, $$ where we assume that both qubits have the same coupling to the cavity. Moreover, the measurement current evolves according to the equation $$ \mathrm{d}I = \sqrt{\frac{8\chi^2}{\kappa(2\bar{n}+1)}}(\langle\sigma_z^1\rangle+\langle\sigma_z^2\rangle)\mathrm{d}t+\mathrm{d}W. $$ See the file Dispersive qubit readout.ipynb and Ref. [1] for details.
Following Ref. [2], the entanglement is generated by letting the system evolve from the state $|\psi_0\rangle = \frac{1}{2}(|0\rangle+|1\rangle)\otimes(|0\rangle+|1\rangle)$ for a time $T_m$ and looking at the integrated current $J(T_m) = \int_0^{T_m}\mathrm{d}t\,I(t)$; if the magnitude of the integrated current is smaller than a threshold value $\nu$, $|J|\le\nu$, the resulting state is the Bell state $|\Psi_+\rangle=\frac{1}{\sqrt{2}}(|01\rangle+|10\rangle)$; otherwise, it is one of the separable states $|00\rangle$, $|11\rangle$. Below, we investigate how the presence of thermal noise affects the resulting entanglement, as quantified by the logarithmic negativity [3], and the success probability of the postselection procedure, extending the results of Ref. [2].
[1] O. Cernotik, D. Vasilyev, and K. Hammerer, arXiv:1503.07484.
[2] C. L. Hutchison, J. Gambetta, A. Blais, and F. K. Wilhelm, Canadian Journal of Physics 87, 225 (2009).
[3] G. Vidal and R. F. Werner, Physical Review A 65, 032314 (2002).
In [1]:
%matplotlib inline
In [2]:
from qutip import qeye, sigmaz, tensor, ket2dm, basis, smesolve, Options, partial_transpose
from matplotlib.pyplot import subplots, plot
from numpy import linspace, cumsum, log2, sqrt, zeros, real
In [3]:
def logneg(rho):
#***** function logneg(rho) *****
#Calculates the logarithmic negativity of a given state.
#Input parameters:
# rho: quantum state.
#Returns:
# EN: logarithmic negativity.
rhoPT = partial_transpose(rho, [0,1])
EN = log2((rhoPT*rhoPT.dag()).sqrtm().tr())
return EN
In [4]:
def entanglement(chi, kappa, nbar, Ntraj, Nsteps, Nsub, tmax):
#***** function entanglement() *****
#Used to calculate the logarithmic negativity, corresponding
#success probability and histograms of integrated currents
#for given system parameters.
#Input parameters:
# chi: qubit-cavity coupling,
# kappa: cavity measurement and decay rate,
# nbar: number of thermal excitations,
# Ntraj: number of generated quantum trajectories,
# Nsteps: number of time steps in the time evolution,
# Nsub: number of substeps,
# tmax: measurement time.
#Returns:
# E: logarithmic negativity after postselection.
# E[i] gives the logarithmic negativity of the state
# obtained by using the threshold value i.
# Psucc: success probability of the postselection,
# Psucc[i] gives probability with threshold i.
# M: data for plotting histogram of the integrated current.
#***** System operators *****
tlist = linspace(0, tmax, Nsteps)
dt = tmax/float(Nsteps)
Id2 = qeye(2)
sz1 = tensor(sigmaz(), Id2)
sz2 = tensor(Id2, sigmaz())
H = 0*tensor(Id2, Id2)
c_op = [sqrt(8*chi*chi*nbar*(nbar+1)/(kappa*(2*nbar+1)))*(sz1+sz2)]
sc_op = [sqrt(2*chi*chi/(kappa*(2*nbar+1)))*(sz1+sz2)]
e_op = [sz1+sz2]
ket_plus = (basis(2,0)+basis(2,1)).unit()
rho_plus = ket2dm(ket_plus)
rho0 = tensor(rho_plus, rho_plus)
#***** Quantum trajectories *****
opts = Options(store_states = True)
E = zeros(100)
Psucc = zeros(100)
M = zeros(Ntraj)
t = zeros(20)
for i in xrange(0, 20):
t[i] = (i+1)*5
for i in xrange(0, Ntraj):
print i
sol = smesolve(H, rho0, tlist, c_op, sc_op, e_op, nsubsteps=Nsub, method='homodyne', solver='fast-milstein',
store_measurement=True, options=opts)
meas = dt*real(cumsum(sol.measurement[0][:,0]))
M[i] = meas[Nsteps-1]
En = logneg(sol.states[0][Nsteps-1])
for j in xrange(0,100):
if (abs(M[i])<=float(j)):
E[j] += En
Psucc[j] += 1
for i in xrange(1,100):
E[i] /= Psucc[i]
Psucc /= float(Ntraj)
E[0] = E[1]
return E, Psucc, M
In [5]:
chi = 0.1
kappa = 1.
nbar = 1.
Ntraj = 200
Nsteps = 1000
Nsub = 250
tmax = 50
E, P, M = entanglement(chi, kappa, nbar, Ntraj, Nsteps, Nsub, tmax)
In [7]:
th = linspace(1, 100, 100)
fig, axes = subplots(nrows=1, ncols=2, figsize=(12,5))
axes[0].plot(th, E, label='Logarithmic negativity')
axes[0].plot(th, P, label='Success probability')
axes[0].set_xlabel('Threshold')
axes[0].set_ylabel('Log. negativity/Succ. probability')
axes[0].legend(loc=0)
axes[1].hist(M, bins=51)
axes[1].set_xlim((-max(abs(M)), max(abs(M))))
axes[1].set_xlabel('Integrated current')
axes[1].set_ylabel('Counts')
Out[7]: