Ondrej Cernotik
24 February 2015
Here we consider a qubit dispersively coupled to a transducer consisting of two coupled oscillators. The first oscillator interacts with the qubit and, at the same time, couples to a thermal bath, while the other oscillator is used to read out the state of the qubit. The system is described by the equation $$ \mathrm{d}\rho = -i[\chi\sigma_zp_1+\omega(a^\dagger a+b^\dagger b)+gq_1q_2,\rho]\mathrm{d}t +\gamma(\bar{n}+1)\mathcal{D}[a]\rho\mathrm{d}t+\gamma\bar{n}\mathcal{D}[a^\dagger]\rho\mathrm{d}t +\kappa\mathcal{D}[b]\rho\mathrm{d}t+\sqrt{\kappa}\mathcal{H}[ib]\rho\mathrm{d}W. $$ Here, the first (i.e., the thermal) oscillator is described by the annihilation operator $a$ and the canonical operators $q_1$, $p_1$, and the second (i.e., the readout) oscillator is described using the operators $b$, $q_2$, $p_2$. The qubit couples to the transducer at the rate $\chi$, the oscillators interact with strength $g$ and both have frequency $\omega$. The thermal oscillator has decay rate $\gamma$ and couples to a bath with $\bar{n}$ thermal excitations, the readout oscillator couples to vacuum bath at rate $\kappa$.
Since this system is too complex to be treated analytically, we follow the formal treatment of Ref. [1] to obtain an effective equation of motion for the qubit. To compare the effective equation with the exact model, we use the trace distance as in the other examples; in the simulations, the unconditional steady state is taken as the initial state of the transducer.
References
[1] O. Cernotik, D. Vasilyev, and K. Hammerer, arXiv:1503.07484.
In [1]:
%matplotlib inline
In [2]:
from qutip import qeye, tensor, destroy, sigmax, sigmay, sigmaz, steadystate, ket2dm, basis, smesolve
from numpy import linspace, array, dot, sqrt, diag
from scipy.linalg import inv, solve_lyapunov, solve_continuous_are
from matplotlib.pyplot import subplots, plot
In [3]:
def trace_dist(chi, omega, g, kappa, gamma, nbar, tmax, Nsteps, Nsub, NFock1, NFock2, Ntraj, rhoq0):
#***** function trace_dist *****
#Input parameters:
# chi: coupling between the qubit and the thermal oscillator,
# omega: frequency of the oscillators,
# g: coupling between the two oscillators,
# kappa: measurement (and decay) rate,
# gamma: decay rate of the thermal oscillator,
# nbar: thermal occupation,
# tmax: maximum time,
# Nsteps: number of time steps in the evolution,
# Nsub: number of substeps,
# NFock1: Fock space cutoff for the thermal oscillator,
# NFock2: Fock space cutoff for the measurement oscillator,
# Ntraj: number of generated trajectories,
# rhoq0: initial qubit state.
#Returns:
# t: time list,
# D: trace distance between the exact and approximate solutions.
#***** System operators *****
tlist = linspace(0, tmax, Nsteps)
Id2 = qeye(2)
IdN1 = qeye(NFock1)
IdN2 = qeye(NFock2)
a = tensor(Id2, destroy(NFock1), IdN2)
xa = (a+a.dag())/sqrt(2.)
pa = 1j*(a.dag()-a)/sqrt(2.)
b = tensor(Id2, IdN1, destroy(NFock2))
xb = (b+b.dag())/sqrt(2.)
pb = 1j*(b.dag()-b)/sqrt(2.)
sx = tensor(sigmax(), IdN1, IdN2)
sy = tensor(sigmay(), IdN1, IdN2)
sz = tensor(sigmaz(), IdN1, IdN2)
#***** Transducer operators for determining the
# transducer steady state *****
a2 = tensor(destroy(NFock1), IdN2)
xa2 = (a2+a2.dag())/sqrt(2.)
b2 = tensor(IdN1, destroy(NFock2))
xb2 = (b2+b2.dag())/sqrt(2.)
#***** Full model *****
H_full = chi*sz*pa+omega*(a.dag()*a+b.dag()*b)+g*xa*xb
c_op_full = [sqrt(gamma*(nbar+1))*a, sqrt(gamma*nbar)*a.dag()]
sc_op_full = [1j*sqrt(kappa)*b]
e_op_full = [sx, sy, sz]
rhoT = steadystate(omega*(a2.dag()*a2+b2.dag()*b2)+g*xa2*xb2,
[sqrt(gamma*(nbar+1))*a2, sqrt(gamma*nbar)*a2.dag(), sqrt(kappa)*b2]) # determining the steady state
rho0 = tensor(rhoq0, rhoT)
#***** Gaussian adiabatic elimination *****
sigma = array([[0,1,0,0], [-1,0,0,0], [0,0,0,1], [0,0,-1,0]])
A = array([[-gamma/2., omega, 0., 0.], [-omega, -gamma/2., -g/2., 0.],
[0., 0., -kappa/2., omega], [-g/2., 0., -omega, -kappa/2.]])
N = diag([gamma*(nbar+0.5), gamma*(nbar+0.5), kappa/2., kappa/2.])
c = sqrt(kappa/2.)*array([[0., 0., 0.,-1.]]).transpose()
m = sqrt(kappa/2.)*array([[0., 0., 1., 0.]]).transpose() # definitions following Ref. [1]
covU = solve_lyapunov(A, -2*N)
covC = solve_continuous_are((A+2*dot(dot(sigma, m), c.transpose())).transpose(), c,
2*(N-dot(dot(sigma, m), dot(m.transpose(), sigma.transpose()))), array([[2]]))
# unconditional and conditional covariance matrices
Q = A-2*dot(dot(covC, c)-dot(sigma, m), c.transpose())
invA = inv(A)
invQ = inv(Q)
s = array([[0., chi*sigmaz(), 0., 0.]]).transpose() # vector of the system operators in the interaction
Lambda = dot(covC-1j*sigma, dot(invQ.transpose(), c))-dot(invA, dot(covC, c)-dot(sigma, m))
M = dot(invA, covU+1j*sigma)-dot(covU-1j*sigma.transpose(), invA.transpose())
Sigma = -0.5*(dot(invA, covU-1j*sigma)+dot(covU+1j*sigma.transpose(), invA.transpose()))
H_gae = 0.25*1j*sum(sum(dot(s.transpose(), dot(M, s))))
c_op_gae = [chi*sqrt(Sigma[1,1]-dot(Lambda, Lambda.conjugate().transpose())[1,1])*sigmaz()]
sc_op_gae = [1j*sum(sum(dot(Lambda.transpose(), s)))]
e_op_gae = [sigmax(), sigmay(), sigmaz()] # Hamiltonian, collapse and measurement operators, evolution operators
#***** Quantum trajectories *****
D = 0
for i in xrange(0, Ntraj):
print i
sol_full = smesolve(H_full, rho0, tlist, c_op_full, sc_op_full, e_op_full, nsubsteps=Nsub, method='homodyne',
solver='fast-milstein', store_measurement=True)
sol_gae = smesolve(H_gae, rhoq0, tlist, c_op_gae, sc_op_gae, e_op_gae, nsubsteps=Nsub, method='homodyne',
solver='fast-milstein', noise=sol_full.noise)
s1x = sol_full.expect[0]
s1y = sol_full.expect[1]
s1z = sol_full.expect[2]
s2x = sol_gae.expect[0]
s2y = sol_gae.expect[1]
s2z = sol_gae.expect[2]
D = D+sqrt((s1x-s2x)**2+(s1y-s2y)**2+(s1z-s2z)**2)
D = D/Ntraj
return tlist, D
In [4]:
chi = 0.2
omega = 5.
g = 0.5
kappa = 1.
gamma = 0.1
nbar = 0.5
tmax = 10.
Nsteps = 500
Nsub = 250
NFock1 = 10
NFock2 = 8
Ntraj = 10
t, D = trace_dist(chi, omega, g, kappa, gamma, nbar, tmax, Nsteps, Nsub, NFock1, NFock2, Ntraj, ket2dm((basis(2,0)+basis(2,1)).unit()))
In [5]:
fig, axes = subplots(figsize=(12,8))
axes.plot(t, D)
axes.set_xlabel('Time')
axes.set_ylabel('Trace distance')
Out[5]:
In [6]:
D_time = D.sum()/Nsteps
print D_time