đŁ Cavity dephasing from transmon thermal excitations and methods to reduce it
import numpy as np
import qutip as qt
import matplotlib.pyplot as plt
import datetime
import csv
from scipy.optimize import curve_fit
from tools import proj
from multiprocessing import Pool
from itertools import product
from matplotlib import ticker, cm
plt.style.use('ggplot')
%matplotlib inline
đ What is the Lindblad master equation?
Normally, we evaluate the time evolution of a quantum system with the Schrodinger equation. Unfortunately, although it is possible to describe everything with the schrodinger equation, it is not always feasible to do so. Take, for example, interaction with the environment outside our experiment, we could model this interaction with the schrodinger equation but to do so we need to be able to model the quantum behaviour of the environment, and that is hard to impossible. A better approach is to use a master equation (the Lindblad master equation, to be precise), it allows usu to model the quantum dynamics of an open system, i.e. non-unitary operation on our quantum system. If we had, for example, a qubit that has some chance to collapse on the ground state every so often (non-coherently, that is) we could model such thing with the master equation.
This is the Lindblad master equation in itâs full form:
\[\dot{\rho} = -\frac{i}{\hbar}[H,\rho] + \sum_{n,m=1}^{N^2-1} h_{nm}(A_n \rho A^\dagger_m - \frac{1}{2}\{A^\dagger_mA_n, \rho\})\]But is better solved (for us) with itâs diagonal form
\[\dot{\rho} = -\frac{i}{\hbar}[H,\rho] + \sum_{i=1}^{N^2-1} \gamma_i(L_i \rho L^\dagger_i - \frac{1}{2}\{L^\dagger_i L_i, \rho\})\]Where $\gamma_i$ is the rate at which this process occurs and $L_i$ are the Lindblad operators (or jump operators),
\[L_i = \sum_{j=1}^{N^2-1}u_{ji}A_j\]$u$ is the unitary transformation that diagnelizes $h$, i.e. $u^\dagger h u$ is a diagonal matrix.
â But we donât care about any of this crap â
Instead, weâll approximate the expression using the Born, Markov, and Secular approximations, and get
\[\boxed{L_A[\rho]=2A O A^\dagger - A^\dagger A O - O A^\dagger A}\] \[\boxed{\dot{\rho} = - \frac{i}{\hbar}[H, \rho] + \sum_{A\in {A}} \frac{\gamma_A}{2}L_A[\rho]}\]Where each ${A}$ is a set of non-unitary operators, which happen with rate $\gamma$.
đš Defining parameters
Most parameters were obtained from Philip Reinholdâs lovely phd dissertation. Also check out Exploring the Quantum by Serge Haroche, mainly chapter 4, to see an explenation of all the parameters, this is a great book in general.
wc = 0 # cavity freq (MHz), it's 0 becuase it doesn't affect the dephasing time, and would make our life easier later
wa = 0 # atom freq (MHz)
chi_e = 2*np.pi*0.05 # Dispersive coupling strength between the cavity and |e> level
chi_f = 2*np.pi*0.1 # Dispersive coupling strength between the cavity and |f> level
kappa = 1e-6 # cavity dissipation rate
gamma = 0.03 # atom dissipation rate
Nc = 2 # number of cavity levels (|0>, |1>, ..., |Nc-1>)
Nq = 3 # number of atom levels (|g>, |e>, ..., |Na-1>)
n_th_a = 0.02 # avg number of atom bath excitation
n_th_c = 0.01 # avg number of cavity bath excitation
mu = [0.999, 0.98, 0.96] # Measurement error rate [|g> error, |e> error |f> error]
# Non linear components (not really important smh :/ )
alpha_T = -2*np.pi*200 # Transmon anharmonicity (E_f - E_e = E_e - E_g - alpha_T) ~ MHz
K = -2*np.pi*4e-3 # Cavity anharmonicity ~ KHz
chi2 = -2*np.pi*20e-3 #
tlist = np.linspace(0, 100, 500) # MUST: max(tlist) >> 1/chi
đš Defining Hamiltonian, initial state and operators
To do the simulations we need to define an initial condition (should be normalized and stuffâŚ)
\[\lvert \psi_{initial} \rangle = (\lvert 0 \rangle + \lvert 1 \rangle) \otimes \lvert g \rangle\]and an Hamiltonian
\[\mathcal{H} = \overbrace{\underbrace{\omega_c \hat{a}^\dagger \hat{a}}_{\text{cavity}} + \underbrace{\omega_a \hat{\sigma}_-^\dagger \hat{\sigma}_-}_{\text{atom}} + \underbrace{\hat{a}^\dagger \hat{a} ( \chi_e \lvert e \rangle \langle e \rvert + \chi_f \lvert f \rangle \langle f \rvert )}_{\text{atom-cavity disspersive coupling}}}^{\text{linear}} + \overbrace{\underbrace{\chi_{nl} {\hat{a}^\dagger}^2 \hat{a}^2 {\hat{\sigma}^\dagger_-}^2 \hat{\sigma}_-^2}_{\text{non-linear coupling}} + \underbrace{\frac{\alpha_T}{2}{\hat{\sigma}^\dagger_-}^2 {\hat{\sigma}_-}^2 + \frac{\mathcal{K}}{2} {\hat{a}^\dagger}^2 {\hat{a}}^2}_{\text{atom/cavity anhamonicity}}}^{\text{non-linear (non-important...)}}\]# intial state
cav0 = (qt.basis(Nc, 1) + qt.basis(Nc, 0))/np.sqrt(2) # |0> + |1>
qub0 = qt.basis(Nq, 0) # |g>
psi0 = qt.ket2dm(qt.tensor(cav0, qub0)) # |0> + |1> x |g>
psi_e = qt.ket2dm(qt.tensor(cav0, qt.basis(Nq, 1))) # |0> + |1> x |e>
psi_f = qt.ket2dm(qt.tensor(cav0, qt.basis(Nq, 2))) # |0> + |1> x |f>
# operators
a = qt.tensor(qt.destroy(Nc), qt.qeye(Nq))
sm = qt.tensor(qt.qeye(Nc), qt.destroy(Nq))
e = qt.tensor(qt.qeye(Nc), qt.basis(Nq, 1))
f = qt.tensor(qt.qeye(Nc), qt.basis(Nq, 2))
ad = a.dag()
smd = sm.dag()
ed = e.dag()
fd = f.dag()
# Dispersive hamiltonian
H = wc*ad*a + wa*smd*sm + ad*a*( chi_e*e*ed + chi_f * f*fd )\
+ chi2*ad*ad*a*a*smd*smd*sm*sm + (alpha_T/2)*smd*smd*sm*sm + (K/2)*ad*ad*a*a # Non-linear
đš Defining collapse operators
Now this is the interesting part. Each operators in the list c_ops will be evaluated as if itâs operation is occouring at some specified rate. So the third one, for example, the $\hat{\sigma}_z$ is responsible for qubit non-coherently relaxing.
In this example, these are the collapse operators:
- Cavity relaxation: $\sqrt{\frac{K}{2}(1 + n_{th})} \hat{a}\ \rightarrow$ The cavity loses a photon due to thermal noise
- Cavity excitation: $\sqrt{\frac{K}{2}n_{th}} \hat{a}^\dagger \quad\quad\rightarrow$ The cavity gains a photon due to thermal noise
- Qubit relaxation: $\sqrt{\frac{\gamma}{2}(1 + n_{th})} \hat{\sigma}_-\ \rightarrow$ The atom loses a photon due to thermal noise
- Qubit excitation: $\sqrt{\frac{\gamma}{2} n_{th}} \hat{\sigma}_+ \quad\ \ \quad\rightarrow$ The atom gains a photon due to thermal noise
You can find this result in Harocheâs book, chapter 4.5
Note: Here, because of QuTip architeture, the occourence rate $\gamma$ is embedded into the operators themself (hence the factor of $\sqrt{\frac{\gamma}{2}}$).
c_ops = [
np.sqrt(0.5*kappa*(1+n_th_c)) * a, # cavity relaxation (Thermal)
np.sqrt(0.5*kappa*n_th_c) * a.dag(), # cavity excitation (Thermal)
np.sqrt(0.5*gamma*(1+n_th_a)) * sm, # qubit relaxation (Thermal)
np.sqrt(0.5*gamma*n_th_a) * sm.dag(), # qubit excitation (Thermal)
]
đš Soooo⌠what are you doing here?
Glad you asked :) We want to decrease the dephasing as much as possible, since the thermal excitations cause the problem, one may cool down the system as much as possible, but obviouly we would reach our limit pretty quickly and need to find new ways.
A different approach would be to measure the qubit repeatedly, and if we measure it to be in the excited state, put it back into the ground state. We can do this by projecting the qubit-cavity state into this state:
\[\rho_{new} \leftarrow (I \otimes \lvert g \rangle\langle g \rvert)\ \rho_{old}\ (I\otimes\lvert g \rangle\langle g \rvert) + (I \otimes \lvert e \rangle\langle e \rvert)\ \rho_{old}\ (I \otimes \lvert e \rangle\langle e \rvert)\]It looked compicated but it comes down to putting the qubit into the ground state, and adding the cavity state of the ground part (of the qubit, i.e. the cavity state if the qubit were grounded) together with the cavity state of the excited part, to make the cavity state unchanged
\[\text{qubit} \otimes\text{cavity} \xrightarrow{\textit{measure}} \text{ground} \otimes (\text{cavity}_{ground} + \text{cavity}_{excited})\]đš Simulate an experiment
With the sim_meas function we simulate an experiment, where psi0 is the initial state, H is the Hamiltonian and dt is the time between measurements (in $\mu s$).
def sim_meas(psi0, H, dt):
"""Simulate the system and do a measurment each dt time step, each measurment project the qubit state to the ground state"""
states = []
for t in np.arange(0, len(tlist), dt):
psi_meas = proj(states[-1], Nc, Nq, mu) if states else psi0
time_range = tlist[t:min(t+dt, len(tlist))]
output = qt.mesolve(H, psi_meas, time_range, c_ops)
states.extend(output.states)
return states
Letâs test it out with some experiment
states = sim_meas(psi0, H, len(tlist))
n = [abs(np.trace(state*(psi0+psi_e+psi_f))) for state in states] # Prob of |0> + |1> (don't be confused by the psi0+psi_e+psi_f)
fig, ax = plt.subplots(figsize=(10,5))
ax.plot(tlist, n) # P(|0>+|1> X |g> + |0>+|1> X |e>)
ax.legend([r"$P(| g \rangle + | e \rangle)$"])
ax.set_xlabel(r'Time [$\mu s$]')
ax.set_ylabel(r'Occupation probability of the $| g \rangle + | f \rangle$ state')
ax.set_title('Dephasing of a superconducting cavity coupled to a transmon');
Fitting the curve
We can use a linear fit to find the life-time of the cavity. Since we know that
\[e^{-\frac{1}{\tau}x} \approx 1 - \frac{1}{\tau}x + \dots\]We can fit the curve with a line $y = mx + b$ and from that get the life time as
\[\tau = -\frac{1}{m}\]Letâs put it into code
delta_t = tlist[-1] - tlist[0]
m = -(n[0]-n[-1])/delta_t # Slope of the line from the first to the last point
plt.figure(figsize=(10,5))
plt.plot(tlist, n)
plt.plot(tlist, m*tlist+n[0], '--')
plt.legend([r"$P(| g \rangle + | e \rangle)$", "fit"])
plt.xlabel(r'Time [$\mu s$]')
plt.ylabel(r'Occupation probability of the $| g \rangle + | f \rangle$ state')
plt.title('Dephasing of a superconducting cavity coupled to a transmon - with linear fit');
plt.text(delta_t/2, 0.8*m*delta_t/2+n[0], f'$y={m*1e4:.2f}\\cdot 10^4t+{n[0]:.2f}$', fontsize=12, color='RoyalBlue')
print(f'The life-time of the cavity is \33[35m\33[1m{-.5/m:.3f}\33[0m')
The life-time of the cavity is 3396.817
Whatâs going on?
note: This is for an old version of the simulation, back when I set the atom and cavity frequencies (wc & wa) to be a number different than zero, so the cavity oscillated around the $\hat{Z}$ axis (I think itâs called Ramsey oscilations or something?). Itâs not the case anymore but the explenation of the dephasing is still valid so I left it here :)
In this example, the transmon starts at the ground state, while the cavity starts at the superpositon $\lvert 0 \rangle + \lvert 1 \rangle$, since we give the cavity energy, the relative phase of the cavity changes and the cavity state ârotatesâ from $\lvert 0 \rangle + \lvert 1 \rangle$ to $\lvert 0 \rangle - \lvert 1 \rangle$ and so forth. While we can easily see the oscilations, it is not so clear why would the cavity decay, after all, we set the cavity dissipation to be 0, it shouldnât have losses, so whatâs going on? Finally, we get to the real problem with dephasing. The reason the cavity oscilation decay is that we loose our knowladge of the frequency of the cavity, and if we donât know the frequency, we canât know after a while weather the cavity is in the state $\lvert 0 \rangle + \lvert 1 \rangle$ or $\lvert 0 \rangle - \lvert 1 \rangle$, so eventally it becomes a constant 50/50% of either state (i.e. the cavity approaches the state $\frac{1}{2}(\lvert 0 \rangle+\lvert 1 \rangle)(\langle 0 \rvert+\langle 1 \rvert) + \frac{1}{2}(\lvert 0 \rangle-\lvert 1 \rangle)(\langle 0 \rvert-\langle 1 \rvert)$).
But why do we loose our knowledge of the cavity frequency? Well, even though the cavity and the transmon are off resonance (so there arenât rabi oscilations between them) they are still coupled. In the disspersive regime, we can think of it as if the state of the transmon affects the frequency of the cavity, so if for example the cavity frequency is $\omega$ when the qubit is at the ground level, it would be $\omega + \chi$ when the transom is excited. And the transmon getâs excited (probabiliticly) due to thermal fluctuation, so the frequency of the cavity changes and we get the dephasing.
The effect of the time between measurments and $\chi_e$ on the dephasing rate
We can plot the dephasing rate (and the life-time) of the cavity vs the time between measurments ($dt$) and $\chi_e$ and we see an interesting behivour. We can see that the life-time approachas infinity when $dt \ll 1/\chi_e$. Or, more precisly, we would see exactly that if we werenât considering the third atom level $\lvert f \rangle$ (since $\chi_f \ne 0$), but more about that later.
def calc_lifetime(tup):
"""Calculate the lifetime of the cavity given the coupling (chi) and the time between measurements (dt). Given in the tuple (dt, chi)"""
dt, chi_e = tup
H = wc*ad*a + wa*smd*sm + ad*a*( chi_e*e*ed + chi_f * f*fd )\
+ chi2*ad*ad*a*a*smd*smd*sm*sm + (alpha_T/2)*smd*smd*sm*sm + (K/2)*ad*ad*a*a # Non-linear
states = sim_meas(psi0, H, dt)
qub = psi0 + psi_e + psi_f # |g> + |e> + |f>
# Initial and final states of the cavity (projected onto the)
init = abs(np.trace(states[0]*qub))
final = abs(np.trace(states[-1]*qub))
m = (final-init)/tlist[-1] # m = dy/dx
return -0.5/m
Calculating the life-time of different $\chi_e$ and $dt$ values
t_tot = 10 # The largest value of dt to simulate (in micro-seconds)
t_max = np.argmin(abs(tlist-t_tot))
dt_pnts = min(20*2, t_max) # Number of values of dt to simulate (can't be more than the amount of points there are in tlist)
chi_pnts = 10 # Number of values of chi_e to simulate
dt_list = np.linspace(2, t_max, dt_pnts, dtype=int)
chi_e_list = np.linspace(0, 0.15, chi_pnts)
param_range = list(product(dt_list,chi_e_list))
with Pool(20) as p:
results = p.map(calc_lifetime, param_range)
tau = np.reshape(results, (dt_pnts, chi_pnts)) # The life time of the cavity
â Windows warning â
Iâm pretty sure that the multi-processing functionality doesnât work on Windows (in the lab I ran it from a Linux subsystem on Windows). If thatâs the case, simply replace ~~~python with Pool(20) as p: results = p.map(calc_lifetime, param_range) ~~~ with ~~~python results = list(map(calc_lifetime, param_range)) ~~~ be warned that it would take much longer though.
If you have a mac, I have no idea what would happen.
Plotting the results
fig, ax=plt.subplots(figsize=(23, 11))
x = chi_e_list
y = np.linspace(tlist[0], tlist[t_max], len(tau))
min_order = np.floor(np.log10(np.min(tau)))
max_order = np.ceil(np.log10(np.max(tau)))
levels = [(1+4*(i%2))*10**(min_order+i//2) for i in range(int(max_order-min_order)*2+1)]
cp1 = ax.contourf(x*1e3, y, tau, levels=levels, locator=ticker.LogLocator(), cmap=cm.Pastel1)
chis = np.linspace(1e3/(t_tot*2*np.pi), chi_e_list[-1]*1e3,100)
ax.plot(chis, 1e3/(2*np.pi*chis), '--', label=r'$\frac{1}{\chi_e}$')
# ax.plot(chis, [1]*len(chis), '--', label=r'$\frac{1}{\chi_e}$')
ax.set_title('Life-time of a cavity coupled to a transmon dispersively', fontsize=25)
ax.set_xlabel(r'Coupling $\chi_e$ [KHz]', fontsize=20)
ax.set_ylabel(r'Time between measurements $dt$ [$\mu s$]', fontsize=20)
ax.legend(fontsize=20)
ax.tick_params(labelsize=15)
cbar = fig.colorbar(cp1);
cbar.ax.tick_params(labelsize=18)
cbar.set_label('Life time of the cavity [$\\mu s$]', fontsize=20)
đš Case when not considering the third (âforbiddenâ) atom level ($\chi_f=0$)
I mentioned before that we need to do the repeated measurments only if we consider the third level. That is becuase if there were only two levels we could set $\chi_e$ to be in the order of $1 KHz$ so the dephasing time would be in the order of second more or less (as youâll see later). Letâs simulate that now to see Iâm not lying (I will not go into details since itâs basically the same thing but now $\chi_f=0$).
def calc_lifetime(tup):
"""Calculate the lifetime of the cavity given the coupling (chi) and the time between measurements (dt). Given in the tuple (dt, chi)"""
dt, chi_e = tup
H = wc*ad*a + wa*smd*sm + chi_e*ad*a*e*ed\
+ chi2*ad*ad*a*a*smd*smd*sm*sm + (alpha_T/2)*smd*smd*sm*sm + (K/2)*ad*ad*a*a # Non-linear
states = sim_meas(psi0, H, dt)
qub = psi0 + psi_e + psi_f # |g> + |e> + |f>
# Initial and final states of the cavity (projected onto the)
init = abs(np.trace(states[0]*qub))
final = abs(np.trace(states[-1]*qub))
m = (final-init)/tlist[-1] # m = dy/dx
return -0.5/m
# ===========================================================================================
t_tot = 10 # The largest value of dt to simulate (in micro-seconds)
t_max = np.argmin(abs(tlist-t_tot))
dt_pnts = min(20*2, t_max) # Number of values of dt to simulate (can't be more than the amount of points there are in tlist)
chi_pnts = 10 # Number of values of chi_e to simulate
dt_list = np.linspace(2, t_max, dt_pnts, dtype=int)
chi_e_list = np.linspace(0, 0.15, chi_pnts)
param_range = list(product(dt_list,chi_e_list))
with Pool(20) as p:
results = p.map(calc_lifetime, param_range)
tau = np.reshape(results, (dt_pnts, chi_pnts)) # The life time of the cavity
# ===========================================================================================
fig, ax=plt.subplots(figsize=(23, 11))
x = chi_e_list
y = np.linspace(tlist[0], tlist[t_max], len(tau))
min_order = np.floor(np.log10(np.min(tau)))
max_order = np.ceil(np.log10(np.max(tau)))
levels = [(1+4*(i%2))*10**(min_order+i//2) for i in range(int(max_order-min_order)*2+1)]
cp1 = ax.contourf(x*1e3, y, tau, levels=levels, locator=ticker.LogLocator(), cmap=cm.Pastel1)
chis = np.linspace(1e3/(t_tot*2*np.pi), chi_e_list[-1]*1e3,100)
ax.plot(chis, 1e3/(2*np.pi*chis), '--', label=r'$\frac{1}{\chi_e}$')
# ax.plot(chis, [1]*len(chis), '--', label=r'$\frac{1}{\chi_e}$')
ax.set_title('Life-time of a cavity coupled to a transmon dispersively', fontsize=25)
ax.set_xlabel(r'Coupling $\chi_e$ [KHz]', fontsize=20)
ax.set_ylabel(r'Time between measurements $dt$ [$\mu s$]', fontsize=20)
ax.legend(fontsize=20)
ax.tick_params(labelsize=15)
cbar = fig.colorbar(cp1);
cbar.ax.tick_params(labelsize=18)
cbar.set_label('Life time of the cavity [$\\mu s$]', fontsize=20)
And just as expected, for low enough $\chi_e$, changing the time between meaurments stops effecting much.