Qubit-CPW Stark Shift

The system described here is a transmon qubit coupled to a voltage-biased CPW. We are going to model the qubit as a multi-level artificial atom coupled to a single harmonic mode of a superconductor resonator by a Jaynes-Cummings Hamiltonian [1]:

$ \newcommand{\ket}[1]{\left|{#1}\right\rangle} \newcommand{\bra}[1]{\left\langle{#1}\right|} $

$\displaystyle H_{JC} = \hbar\omega_r (a^\dagger a + 1/2) + \hbar\Sigma\omega_j \ket{j}\bra{j} + \hbar\Sigma g_{j,j+1}(a^\dagger\ket{j}\bra{j+1}+a\ket{j+1}\bra{j}) $

where:

$\omega_r$ is the bare resonator frequency
$\omega_j$ is the frequency of the $|j>\rightarrow |j+1>$ transition
$g_{j,j+1}$ is the coupling strength

This Hamiltonian can be diagonalized in the dispersive limit, $\Delta_{j,j+1} \equiv \omega_{j,j+1} - \omega_r \gg g_{j,j+1}$, resuting in:

$\displaystyle \tilde{H}^{2}_{JC} = \hbar\tilde{\omega}_r (a^\dagger a + 1/2) + \dfrac{\hbar\tilde{\omega}_{ge}}{2}\sigma_z + \hbar\chi (a^\dagger a + 1/2) \sigma_z$

where:

$\chi \simeq \chi_{ge}-\dfrac{\chi_{ef}}{2} $
$\chi_{j,j+1} \equiv \dfrac{g^2_{j,j+1}}{\Delta_{j,j+1}}$
$\tilde{\omega}_{r} \simeq \omega_r - \dfrac{\chi_{ef}}{2}$

1.Jaynes, E. T. & Cummings, F. W. Comparison of quantum and semiclassical radiation theories with application to the beam maser. Proceedings of the IEEE 51, 89–109 (1963).



In [1]:

    
#import functions
%pylab inline
from MyUnits import *
from MyFunctions import *
from qutip import *
from MyQubit import *
import mpld3









    



Populating the interactive namespace from numpy and matplotlib



In [2]:

    
# %load_ext autoreload



In [3]:

    
# %aimport MyQubit
# %autoreload 1

Functions

Hamiltonian in the rotating frame of both drives.

$ \displaystyle H_{tot} \approx \hbar \tilde{\Delta}_d\left(a^\dagger a + 1/2\right) + \dfrac{\hbar}{2}\chi(\left(a^\dagger a + 1/2\right)\sigma_z + H_{Kerr}^4+H_{Drive}^2$

$H_{Kerr}^{(4)} \simeq \hbar \zeta \left(a^\dagger a + 1/2\right)^2\sigma_z + \hbar\zeta'\left(a^\dagger a + 1/2\right)^2$

$H_{drive}^{(2)} = \dfrac{\hbar}{2}\Omega_d (a^\dagger + a) + \dfrac{\hbar}{2}\Omega_s (\sigma^+ + \sigma^-) $

Function add to file MyQubit.py

Import data



In [4]:

    
cd /Users/rouxinol/Dropbox/Lab LaHaye/Hugo/2014/12042014









    



/Users/rouxinol/Dropbox/Lab LaHaye/Hugo/2014/12042014



In [5]:

    
!Head CavityAttn\ MWFreq\ swp\ Flux2pt8\ b

P,F,A,P= ImportXYZZ("/Users/rouxinol/Dropbox/Lab LaHaye/Hugo/2014/12042014/CavityAttn MWFreq swp Flux2pt8 b")

fig, axes = subplots(1,1, figsize=(16,6))
for i in range(6):
    plot(F,A[i,])

A0 = dBmtoV(A[0,])
A1 = dBmtoV(A[1,])
A2 = dBmtoV(A[2,])
A3 = dBmtoV(A[3,])
A4 = dBmtoV(A[4,])
A5 = dBmtoV(A[5,])

Estimated qubit parameters from fabrication

From the measurement of the resistence of one junction the maximum Josephson Energy for the qubit was estimated.



In [6]:

    
Ej = 11.55 # Maximum Josephson Energy

Ec = 0.22 # Capacitive Energy

w_ge_max = sqrt(8*Ec*Ej)-Ec
print('Qubit max frequency: ',w_ge_max, 'GHz')

f = 5e9 # GHz
T = 50e-3 # K
n_th_a = 1/(exp(h*f/(kB*T)-1))

print('<n> thermal: ',n_th_a)









    



Qubit max frequency:  4.28865833702 GHz
<n> thermal:  0.02238770514

System Parameters

From the measurements take in December 2014 the mains system parametes were determined

The resonator bare frequency is $\omega_r$ = 5.049 GHz, with a internal quality factor $Q_I \simeq 200k$, and a loaded quality facot $Q_L \simeq 20k$.

From the resistence of the Josephson resistence a maximum Josephson energy $E_{J,max}/h \approx 11.92$ GHz and a Coulomb charging energy of $E_c/h \approx $ 0.22 GHz (See comments above)

Using these parametes $\chi_{ef}$, $\chi_{ge}$ and $\omega_{ge}$ can be estimated using:

$\chi_{ef} \simeq 2\left(\omega_r - \tilde{\omega}_r\right) $

and

$\chi_{ge} \simeq \chi + \dfrac{\chi_{ef}}{2}$

and finally:

$\omega_{ge} \simeq \tilde{\omega}_{ge} - \chi_{ge}$

$\omega{ef}$ was estimated using two photon spectroscopy of que qubit. We observe $\dfrac{\tilde{\omega}_{gf}}{2} \simeq $ 4.245 GHz one gets $\omega_{gf} \simeq$ 8.49 GHz; using:

$\tilde{\omega}_{gf} \simeq \tilde{\omega}_{ge}+\tilde{\omega}_{ef}$

$\tilde{\omega}_{ef} \simeq \tilde{\omega}_{gf}-\tilde{\omega}_{ge} \approx \omega_{ef}$

Recovering data from measurements realized December 2014



In [7]:

    
# Measured data
w_r = (5.049) *2 * pi # bare resonator frequency

Q_l = 20000 # resonator loaded quality factor

w_r_tilde = (5.065) * 2 * pi # max cavity frequency change

chi = -(0.0054) * 2 * pi # Change in the spectroscopy peaks -> 2*n*Chi = 4.35 - 4.3

w_ge_tilde = (4.365) *2 * pi # Dressed qubit frequency |g> -> |e>


w_gf_tilde = 2 * (4.245) *2 * pi # Dressed qubit frequency |g> -> |f> 


# Derived data

chi_ef = 2*(w_r - w_r_tilde)
print('Chi_ef/2pi = ',chi_ef/2/pi, 'GHz')

chi_ge = chi + chi_ef/2
print('Chi_ge/2pi = ', chi_ge/2/pi, 'GHz\n')

w_ge = w_ge_tilde - chi_ge
print('qubit bare frequency: w_ge/2pi =', w_ge/2/pi, 'GHz')

w_ef = w_gf_tilde - w_ge_tilde # qubit e -> f transition 
print('qubit |e> -> |f> frequency: w_ef/2pi =', w_ef/2/pi, 'GHz\n')

Delta_ge = w_ge - w_r # g -> e detunning
print('Detunning |g> -> |e>: Delta_ge/2pi =', Delta_ge/2/pi,'GHz')

Delta_ef = w_ef - w_r # g -> e detunning
print('Detunning |e> -> |f>: Delta_ef/2pi =', Delta_ef/2/pi,'GHz\n')

g_ge = sqrt(chi_ge * Delta_ge)
print('Coupling strength |g> -> |e>: g_ge/2pi =', (g_ge/2/pi),'GHz')

g_ef = sqrt(chi_ef * Delta_ef)
print('Coupling strength |e> -> |f>: g_ef/2pi =', (g_ef/2/pi),'GHz\n')









    



Chi_ef/2pi =  -0.03200000000000027 GHz
Chi_ge/2pi =  -0.021400000000000134 GHz

qubit bare frequency: w_ge/2pi = 4.3864 GHz
qubit |e> -> |f> frequency: w_ef/2pi = 4.125 GHz

Detunning |g> -> |e>: Delta_ge/2pi = -0.6626 GHz
Detunning |e> -> |f>: Delta_ef/2pi = -0.9240000000000008 GHz

Coupling strength |g> -> |e>: g_ge/2pi = 0.119078293572 GHz
Coupling strength |e> -> |f>: g_ef/2pi = 0.171953482082 GHz

Pertubative expansion parameter for dispersive regime



In [8]:

    
# lambda expansion parameter
lambda_ge = g_ge/Delta_ge # expansion parameter lambda << 1
print('Expansion parameter |g> -> |e> << 1: ', lambda_ge, ' << 1')

lambda_ef = g_ef/Delta_ef # expansion parameter lambda << 1
print('Expansion parameter |e> -> |f> << 1: ', lambda_ef, ' << 1 \n')









    



Expansion parameter |g> -> |e> << 1:  -0.179713693891  << 1
Expansion parameter |e> -> |f> << 1:  -0.18609684208  << 1

Kerr-type nonlinearities

The general drive Hamiltonian $\displaystyle H_{drive}$ in the dispersive regime result in a Kerr-type nonlinearities Hamiltonian $H_{Kerr}^{(4)}$ and a Hamiltonian to the cavity and qubit tones $H_{drive}^{(2)}$, where:

$H_{Kerr}^{(4)} \simeq \hbar \zeta \left(a^\dagger a + 1/2\right)^2\sigma_z + \hbar\zeta'\left(a^\dagger a + 1/2\right)^2$

with

$\zeta \approx \chi_{ef}\lambda_{ef}^2 - 2\chi_{ge}\lambda_{ge}^2 + \dfrac{7\chi_{ef}\lambda_{ge}^2}{4} -\dfrac{5\chi_{ge}\lambda_{ef}^2}{4}$

and

$ \zeta' \approx \left(\chi_{ge}-\chi_{ef}\right)\left(\lambda_{ge}^2+\lambda_{ef}^2\right)$

and

$H_{drive}^{(2)} = \dfrac{\hbar}{2}\Omega_d (a^\dagger + a) + \dfrac{\hbar}{2}\Omega_s (\sigma^+ + \sigma^-) $

where

$\Omega_s/\Omega_d$ are the amplitudes of the qubit/resonator.



In [9]:

    
zeta_1 = chi_ef * lambda_ef**2 - 2 * chi_ge * lambda_ge**2 + (7 * chi_ef*lambda_ge**2 - 5 * chi_ge * lambda_ef**2)/4
print('cross-Kerr coeficient: Zeta_1/2pi =', zeta_1/2/pi, 'GHz')

zeta_2 = (chi_ge - chi_ef ) * (lambda_ge**2 + lambda_ef**2)
print('self-Kerr coeficient: Zeta_2/2pi =', zeta_2/2/pi, 'GHz\n')









    



cross-Kerr coeficient: Zeta_1/2pi = -0.000608138737206 GHz
self-Kerr coeficient: Zeta_2/2pi = 0.000709447891881 GHz

Qubit-Cavity characteristic times

Relaxation and Dephasing of the qubit and the resonator dissipation are accounted in the master equation through:

Resonator photon dissipation: $\kappa_n = \dfrac{\omega_r}{Q_l}$
Qubit relaxation: $\Gamma_n$
Qubit dephasing: $\gamma_{\varphi}$

The effects of the finite temperature are incorporated to the master equation through:

Resonator thermal excitations: $\kappa_p$
Qubit thermal excitations: $\Gamma_p$

where $\dfrac{\kappa_p}{\kappa_n}\simeq\exp\left({\dfrac{\hbar\omega_r}{k_BT}}\right)$ and $\dfrac{\Gamma_p}{\Gamma_n}\simeq\exp\left({\dfrac{\hbar\omega_{ge}}{k_BT}}\right)$

Qubit relaxation time was measured to be: $T_1 = 4.3 \mu$s

From Rabi measurements the Rabi decay time was measure to be: $T_{Rabi} = 5.8 \mu$s, and Rabi frequency of: $f_{Rabi} = 2.42$MHz with $\bar{n}_s \simeq 1500$ photons

Using

$T_1 = \dfrac{1}{\Gamma_+ + \Gamma_-}$
$\dfrac{1}{T_2^*}=\dfrac{1}{2T_1}+\dfrac{1}{T_{\varphi }}$



In [10]:

    
Temp = 20e-3 # temperature in Kelvin


T_1 = 4300 # qubit relaxation time
print('T_1=',T_1*1e-9, 's\n')

Delta_w_ge = (0.0002) *2 * pi # Spectroscopy HWHM lowest <n> Stark shift peak
T_2_star = 1/Delta_w_ge 
print('T_2_star =', T_2_star*1e-9, 's\n')

T_rabi = 5800 # Rabi decay time
print('T_rabi=',T_rabi*1e-9, 's')

w_rabi = (0.00248) * 2 * pi # Rabi frequency time
print('w_rabi/2pi =',w_rabi/2/pi, 'GHz\n')

T_dep = 2*T_2_star*T_1/(2*T_1 - T_2_star)
print('T_dephasing =', T_dep*1e-9 , 's\n')


kappa_n = w_r/Q_l # resonator dissipation
print('Resonator dissipation k-/2pi =',kappa_n/2/pi, 'GHz\n')

gamma_rel = 1/T_1
print('Qubit Relaxarion: Gamma_rel =',gamma_rel, 'GHz\n')



gamma_dep = 1/T_dep
print('Qubit dephasing: Gamma_def =',gamma_dep, 'GHz\n')









    



T_1= 4.3e-06 s

T_2_star = 7.957747154594768e-07 s

T_rabi= 5.8e-06 s
w_rabi/2pi = 0.00248 GHz

T_dephasing = 8.769176059676528e-07 s

Resonator dissipation k-/2pi = 0.00025245000000000004 GHz

Qubit Relaxarion: Gamma_rel = 0.00023255813953488373 GHz

Qubit dephasing: Gamma_def = 0.0011403579916684755 GHz



In [11]:

    
# basic parameters

w_d = 5.066 * 2 * pi # resonator drive frequency

w_s = 4.36 * 2 * pi # qubit drive frequency

# derivated parameters

# Delta tilde parameters

Delta_d_tilde = w_r_tilde - w_d
print('Detuning of the dressed resonator/2pi: ', Delta_d_tilde/2/pi)

Delta_s_tilde = w_ge_tilde - w_s
print('Detuning of the dressed qubit/2pi: ', Delta_s_tilde/2/pi, '\n')









    



Detuning of the dressed resonator/2pi:  -0.0009999999999993723
Detuning of the dressed qubit/2pi:  0.005000000000000254

Power vs. cavity occupancy

The cavity mean occupancy is:

$$ \bar{n}=\dfrac{P_{out}Q}{2\pi f hf} $$

Using Blais et al. (PRA 062320 (2004))

$$ \left<n\right> = \dfrac{2P}{\hbar\omega_r\kappa} $$

Also from Blais et al.

$$ \left< V_{out} \right> = \sqrt{R\hbar\omega_r\kappa}\dfrac{\left(a^{\dagger}+a\right)}{2} $$



In [12]:

    
p_out = linspace(-111,-151,9)

f = w_r/2/pi*1e9
Q_c = 20000 # External coupling # from arxiv 1308.2245v1 Suri et alli (??? add by hand)
n_mean = dBmtoW(p_out)*Q_l/(pi*f*h*f) # standard equation 
n_mean_2 = 2*Q_l*dBmtoW(p_out)/(hbar*w_r_tilde*1e9*kappa_n*1e9)/Q_c # from arxiv 1308.2245v1 Suri et alli
n_mean_3 = Q_c*kappa_n*1e9*dBmtoW(p_out)/(2*Q_l*hbar*w_r_tilde*1e9)/((kappa_n*1e9/2)**2+(Delta_d_tilde*1e9-chi*1e9)**2) # from arxiv 1308.2245v1 Suri et alli
n_mean_4 = 2*dBmtoW(p_out)/(hbar*w_r*1e9*kappa_n*1e9)



In [13]:

    
print('Power   =>     <n>    => <n2>   => <n3> => <n4> ')
for i in range(len(p_out)):
    print(p_out[i], ' => ',n_mean[i],' => ', n_mean_2[i] ,' => ', n_mean_3[i],' => ', n_mean_4[i])









    



Power   =>     <n>    => <n2>   => <n3> => <n4> 
-111.0  =>  2993.73524412  =>  2984.27823249  =>  2.45395984213  =>  2993.73524412
-116.0  =>  946.702208294  =>  943.711638633  =>  0.776010238772  =>  946.702208294
-121.0  =>  299.373524412  =>  298.427823249  =>  0.245395984213  =>  299.373524412
-126.0  =>  94.6702208294  =>  94.3711638633  =>  0.0776010238772  =>  94.6702208294
-131.0  =>  29.9373524412  =>  29.8427823249  =>  0.0245395984213  =>  29.9373524412
-136.0  =>  9.46702208294  =>  9.43711638633  =>  0.00776010238772  =>  9.46702208294
-141.0  =>  2.99373524412  =>  2.98427823249  =>  0.00245395984213  =>  2.99373524412
-146.0  =>  0.946702208294  =>  0.943711638633  =>  0.000776010238772  =>  0.946702208294
-151.0  =>  0.299373524412  =>  0.298427823249  =>  0.000245395984213  =>  0.299373524412

Drive amplitude qubit

$\Omega_d \simeq \left(\dfrac{Q_C\kappa_- P_{RF}}{2Q_L\hbar\tilde{\omega}_r}\right)^{1/2}$



In [14]:

    
P_qubit = -140 # Input power qubit (0-20-20-76 - cavity(?))

P_cavity = -(106 + 32) # Input power cavity

Gamma_d = sqrt(Q_c*kappa_n*1e9*dBmtoW(P_cavity)/(2*Q_l*hbar*w_r_tilde*1e9))/1e9  # Drive amplitude qubit
print('Drive amplitude cavity =',Gamma_d/2/pi, 'GHz')

Gamma_s = sqrt(Q_c*kappa_n*1e9*dBmtoW(P_qubit)/(2*Q_l*hbar*w_ge_tilde*1e9))/1e9 # Drive amplitude cavity ????
print('Drive amplitude qubit =',Gamma_s/2/pi,'GHz')









    



Drive amplitude cavity = 0.000308010154756 GHz
Drive amplitude qubit = 0.0002635497665 GHz

Simulation



In [15]:

    
tlist = None

Plot <\Sigma_z>, <\n> as a function Pump frequency



In [101]:

    
N = 20 # number of states cavity

f_d_i = 5.058 # cavity initial frequency 
f_d_f = 5.062# cavity final frequency
steps_d = 200
w_d = 2 * pi * linspace(f_d_i,f_d_f,steps_d)

f_s_i = 4.2 # qubit initial freq
f_s_f = 4.4  # qubit final freq
steps_s = 300
w_s = 2 * pi * linspace(f_s_i,f_s_f,steps_s)

amp_sz = []
amp_n = []
amp_x = []
for wr in w_d:
    
    for wq in w_s:
        
        Delta_d_tilde = w_r_tilde - wr
        Delta_s_tilde = w_ge_tilde - wq
        rho,temp, temp2, temp3 = calc_spectrum(N,tlist,Delta_d_tilde,Delta_s_tilde,chi,Gamma_d,Gamma_s, zeta_1,zeta_2,kappa_n,gamma_rel,gamma_dep,n_th_a)
        amp_sz.append((temp.tr()))
        amp_n.append((temp2.tr()))
        amp_x.append((temp3.tr()))



In [102]:

    
Tr_sz = reshape(amp_sz,(-1,len(w_s)))
Tr_n = reshape(amp_n,(-1,len(w_s)))



In [102]:



In [103]:

    
fig, ax = subplots(1,2, figsize=(16,10))
im = ax[0].pcolor(w_s/2/pi, w_d/2/pi ,(abs(Tr_sz)))#,vmin=0, vmax=1)
ax[0].set_xlim(f_s_i,f_s_f)
ax[0].set_ylim(f_d_i,f_d_f)
ax[0].set_ylabel(r'Cavity tone frequency (GHz)')
ax[0].set_xlabel(r'Qubit tone frequency (GHz)')

fig.colorbar(im, ax=ax[0])

im2 = ax[1].pcolor(w_s/2/pi, w_d/2/pi ,(abs(Tr_n)))#,vmin=0, vmax=1)
fig.colorbar(im2, ax=ax[1])
ax[1].set_xlim(f_s_i,f_s_f)
ax[1].set_ylim(f_d_i,f_d_f)
ax[1].set_ylabel(r'Cavity tone frequency (GHz)')
ax[1].set_xlabel(r'Qubit tone frequency (GHz)')









    Out[103]:





<matplotlib.text.Text at 0x11ef8cd30>



In [ ]:

    
# mpld3.display(fig)

Plot <\Sigma_z>, <\n> as a function Spectroscopy power (qubit tone)



In [97]:

    
N = 20 # number of states cavity


P_cavity = -(106 + 32) # Input power cavity

Gamma_d = sqrt(Q_c*kappa_n*1e9*dBmtoW(P_cavity)/(2*Q_l*hbar*w_r_tilde*1e9))/1e9  # Drive amplitude qubit
print('Drive amplitude cavity =',Gamma_d/2/pi, 'GHz')

# w_s = 2 * pi * linspace(f_s_i,f_s_f,steps_s)

wr = (5.0596)  * 2* pi
Delta_d_tilde = w_r_tilde - wr 

# qubit freq scan
f_s_i = 4.2 # qubit initial freq
f_s_f = 4.4  # qubit final freq
steps_s = 200
w_s = 2 * pi * linspace(f_s_i,f_s_f,steps_s)

# Power scan
P_i = -110 # qubit tone start
P_f = -160 # qubit tone final
steps_p = 50 # 
P_qubit = linspace(P_i,P_f, steps_p)




amp_sz = []
amp_n = []
amp_x = []
amp_a = []
for Pq in P_qubit:
    
    for wq in w_s:
        
        
        Delta_s_tilde = w_ge_tilde - wq
        Gamma_s = sqrt(Q_c*kappa_n*1e9*dBmtoW(Pq)/(2*Q_l*hbar*w_ge_tilde*1e9))/1e9

        rho,temp, temp2, temp3, temp4 = calc_spectrum(N,tlist,Delta_d_tilde,Delta_s_tilde,chi,Gamma_d,Gamma_s, zeta_1,zeta_2,kappa_n,gamma_rel,gamma_dep,n_th_a)
        amp_sz.append((temp.tr()))
        amp_n.append((temp2.tr()))
        amp_x.append((temp3.tr()))
        amp_a.append((temp4.tr()))









    



Drive amplitude cavity = 0.000308010154756 GHz



In [98]:

    
Tr_sz = reshape(amp_sz,(-1,len(w_s+1)))
Tr_n = reshape(amp_n,(-1,len(w_s+1)))
Tr_x = reshape(amp_x,(-1,len(w_s+1)))
Tr_a = reshape(amp_a,(-1,len(w_s+1)))



In [99]:

    
Tr_sz.shape, P_qubit.shape,w_s.shape









    Out[99]:





((50, 200), (50,), (200,))



In [100]:

    
fig, ax = subplots(2,2, figsize=(16,10))

im = ax[0,0].pcolor( w_s/2/pi,P_qubit,abs(Tr_sz))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax[0,0])
ax[0,0].set_xlim(4.27,4.39)
ax[0,0].set_ylim(P_i,P_f)
ax[0,0].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[0,0].set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=10)
ax[0,0].set_title(r'$Tr[\rho\sigma_z]$',fontsize=20)


im2 = ax[0,1].pcolor(w_s/2/pi,P_qubit,abs(Tr_n))#,vmin=0, vmax=1)
fig.colorbar(im2, ax=ax[0,1])
ax[0,1].set_ylim(P_i,P_f)
ax[0,1].set_xlim(4.27,4.39)
ax[0,1].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[0,1].set_xlabel(r'Qubit Tone Frequency (GHz)', fontsize=10)
ax[0,1].set_title(r'$Tr[\rho {a^{\dagger}a}]$',fontsize=20)

im3 = ax[1,0].pcolor(w_s/2/pi,P_qubit,-abs(Tr_x))#,vmin=0, vmax=1)
fig.colorbar(im2, ax=ax[1,0])
ax[1,0].set_ylim(P_i,P_f)
ax[1,0].set_xlim(4.27,4.39)
ax[1,0].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[1,0].set_xlabel(r'Qubit Tone Frequency (GHz)', fontsize=10)
ax[1,0].set_title(r'$Tr[\rho (a^{\dagger} + a)]$',fontsize=20)

im3 = ax[1,1].pcolor(w_s/2/pi,P_qubit,abs(Tr_a))#,vmin=0, vmax=1)
fig.colorbar(im2, ax=ax[1,1])
ax[1,1].set_ylim(P_i,P_f)
ax[1,1].set_xlim(4.27,4.39)
ax[1,1].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[1,1].set_xlabel(r'Qubit Tone Frequency (GHz)', fontsize=10)
ax[1,1].set_title(r'$Tr[\rho a]$',fontsize=20)









    Out[100]:





<matplotlib.text.Text at 0x10f5656d8>



In [101]:

    
# mpld3.display(fig)
fig, ax = subplots(1,1, figsize=(16,10))

im = ax.pcolor( w_s/2/pi,P_qubit,abs(Tr_sz))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax)
ax.set_xlim(4.27,4.39)
ax.set_ylim(P_i,P_f)
ax.set_ylabel(r'Qubit Tone Power(dBm)',fontsize=20)
ax.set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=20)
ax.set_title(r'$Tr[\rho\sigma_z]$',fontsize=20)









    Out[101]:





<matplotlib.text.Text at 0x12a2ae128>

<\Sigma_z>, <\n> as a function cavity power (cavity tone)



In [61]:

    
N = 20 # number of states cavity

P_qubit = -140 # Input power qubit (0-20-20-76 - cavity(?))

# P_cavity = -(106 + 30) # Input power cavity

# Gamma_d = sqrt(Q_c*kappa_n*1e9*dBmtoW(P_cavity)/(2*Q_l*hbar*w_r_tilde*1e9))/1e9  # Drive amplitude qubit
# print('Drive amplitude cavity =',Gamma_d/2/pi, 'GHz')

Gamma_s = sqrt(Q_c*kappa_n*1e9*dBmtoW(P_qubit)/(2*Q_l*hbar*w_ge_tilde*1e9))/1e9 # Drive amplitude cavity ????
print('Drive amplitude qubit =',Gamma_s/2/pi,'GHz')



# # qubit tone 
# f_s_i = 4.3 # qubit initial freq
# f_s_f = 4.  # qubit final freq
# steps_s = 50 # steps qubit tone

# w_s = 2 * pi * linspace(f_s_i,f_s_f,steps_s)

wr = (5.0596)  * 2* pi
Delta_d_tilde = w_r_tilde - wr 

# qubit freq scan
f_s_i = 4.2 # qubit initial freq
f_s_f = 4.4  # qubit final freq
steps_s = 500
w_s = 2 * pi * linspace(f_s_i,f_s_f,steps_s)

# Power scan
P_i = -160 # qubit tone start
P_f = -130 # qubit tone final
steps_p = 50 # 
P_cavity = linspace(P_i,P_f, steps_p)



rho_matrix = []
amp_sz = []
amp_n = []
amp_x = []
amp_a = []
for Pc in P_cavity:
    
    for wq in w_s:
        
        
        Delta_s_tilde = w_ge_tilde - wq
        Gamma_d = sqrt(Q_c*kappa_n*1e9*dBmtoW(Pc)/(2*Q_l*hbar*w_r_tilde*1e9))/1e9  # Drive amplitude qubit

        rho,temp, temp2, temp3, temp4 = calc_spectrum(N,tlist,Delta_d_tilde,Delta_s_tilde,chi,Gamma_d,Gamma_s, zeta_1,zeta_2,kappa_n,gamma_rel,gamma_dep,n_th_a)
#         rho_matrix.append(rho)
        amp_sz.append((temp.tr()))
        amp_n.append((temp2.tr()))
        amp_x.append((temp3.tr()))
        amp_a.append((temp4.tr()))









    



Drive amplitude qubit = 0.0002635497665 GHz



In [62]:

    
Tr_sz = reshape(amp_sz,(-1,len(w_s+1)))
Tr_n = reshape(amp_n,(-1,len(w_s+1)))
Tr_x = reshape(amp_x,(-1,len(w_s+1)))
Tr_a = reshape(amp_a,(-1,len(w_s+1)))



In [63]:

    
Tr_sz.shape,Tr_n.shape,Tr_x.shape, Tr_a.shape,w_s.shape,P_cavity.shape









    Out[63]:





((50, 500), (50, 500), (50, 500), (50, 500), (500,), (50,))



In [64]:

    
fig, ax = subplots(2,2, figsize=(16,10))

im = ax[0,0].pcolor(w_s/2/pi,P_cavity,abs(Tr_sz))#,vmin=0, vmax=1)
fig.colorbar(im, ax=ax[0,0])
ax[0,0].set_xlim(4.27,4.39)
ax[0,0].set_ylim(P_i,P_f)
ax[0,0].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[0,0].set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=10)
ax[0,0].set_title(r'$Tr[\rho\sigma_z]$',fontsize=20)


im2 = ax[0,1].pcolor(w_s/2/pi,P_cavity,abs(Tr_n))#,vmin=0, vmax=1)
fig.colorbar(im2, ax=ax[0,1])
ax[0,1].set_ylim(P_i,P_f)
ax[0,1].set_xlim(4.27,4.39)
ax[0,1].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[0,1].set_xlabel(r'Qubit Tone Frequency (GHz)', fontsize=10)
ax[0,1].set_title(r'$Tr[\rho {a^{\dagger}a}]$',fontsize=20)

im3 = ax[1,0].pcolor(w_s/2/pi,P_cavity,abs(Tr_x))#,vmin=0, vmax=1)
fig.colorbar(im2, ax=ax[1,0])
ax[1,0].set_ylim(P_i,P_f)
ax[1,0].set_xlim(4.27,4.39)
ax[1,0].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[1,0].set_xlabel(r'Qubit Tone Frequency (GHz)', fontsize=10)
ax[1,0].set_title(r'$Tr[\rho (a^{\dagger} + a)]$',fontsize=20)

im3 = ax[1,1].pcolor(w_s/2/pi,P_cavity,abs(Tr_a))#,vmin=0, vmax=1)
fig.colorbar(im2, ax=ax[1,1])
ax[1,1].set_ylim(P_i,P_f)
ax[1,1].set_xlim(4.27,4.39)
ax[1,1].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=10)
ax[1,1].set_xlabel(r'Qubit Tone Frequency (GHz)', fontsize=10)
ax[1,1].set_title(r'$Tr[\rho a]$',fontsize=20)


# fig, ax = subplots(1,3, figsize=(16,5))

# im = ax[0].pcolor( w_s/2/pi,P_cavity,abs(Tr_sz))#,vmin=0, vmax=1)
# fig.colorbar(im, ax=ax[0])
# ax[0].set_xlim(4.27,4.39)
# ax[0].set_ylim(P_i,P_f)
# ax[0].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=20)
# ax[0].set_xlabel(r'Qubit Tone Frequency (GHz)',fontsize=20)
# ax[0].set_title(r'$Tr[\rho\sigma_z]$',fontsize=20)

# im2 = ax[1].pcolor(w_s/2/pi,P_cavity,abs(Tr_n))#,vmin=0, vmax=1)
# fig.colorbar(im2, ax=ax[1])
# ax[1].set_ylim(P_i,P_f)
# ax[1].set_xlim(4.27,4.39)
# ax[1].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=20)
# ax[1].set_xlabel(r'Qubit Tone Frequency (GHz)', fontsize=20)
# ax[1].set_title(r'$Tr[\rho {a^{\dagger}a}]$',fontsize=20)

# im3 = ax[2].pcolor(w_s/2/pi,P_cavity,abs(Tr_x))#,vmin=0, vmax=1)
# fig.colorbar(im2, ax=ax[1])
# ax[2].set_ylim(P_i,P_f)
# ax[2].set_xlim(4.27,4.39)
# ax[2].set_ylabel(r'Qubit Tone Power(dBm)',fontsize=20)
# ax[2].set_xlabel(r'Qubit Tone Frequency (GHz)', fontsize=20)
# ax[2].set_title(r'$Tr[\rho {a^{\dagger} + a}]$',fontsize=20)









    Out[64]:





<matplotlib.text.Text at 0x121e0d0b8>



In [95]:

    
# Test
t = 20
print(P_cavity[t])
File = A0

fig, axes = subplots(1,1, figsize=(16,6),sharex=True)

# fig.subplots_adjust(left=0, wspace=0.25

axes.plot(1e9*w_s/2/pi+4e6,-(Tr_x[t]),'k-',linewidth=3)
axes.set_xlabel(r'f (GHz)', fontsize=18)
axes.set_ylabel(r'$Tr[\rho \sigma_z]$ or $Tr[\rho a^{\dagger}a]$ ', fontsize=18)
axes.set_title(r'$Tr[\rho \sigma_z]$',fontsize=20)
#####
axes0 = axes.twinx()
axes0.plot(F,(File)*1e6,'ro', label='data' )

# axes.legend(loc=2)
#left plot
# axes.set_ylim(0.3,1.5) # simulation 
# axes0.set_ylim(0,2) # data

axes0.set_xlim(4.25*1e9,4.40*1e9);
plot(w_s/2/pi,abs(Tr_sz[t]))









    



-147.755102041






    Out[95]:





[<matplotlib.lines.Line2D at 0x1253b6c50>]



In [96]:

    
###A0
t = 30
print(P_cavity[t])
File = A0

fig, axes = subplots(1,1, figsize=(16,6),sharex=True)

# fig.subplots_adjust(left=0, wspace=0.25

axes.plot(1e9*w_s/2/pi+4e6,abs(Tr_sz[t]),'k-',linewidth=3)
axes.set_xlabel(r'f (GHz)', fontsize=18)
axes.set_ylabel(r'$Tr[\rho \sigma_z]$ or $Tr[\rho a^{\dagger}a]$ ', fontsize=18)
axes.set_title(r'$Tr[\rho \sigma_z]$',fontsize=20)
#####
axes0 = axes.twinx()
axes0.plot(F,(File)*1e6,'ro', label='data' )

# axes.legend(loc=2)
#left plot
axes.set_ylim(0.3,1.5) # simulation 
axes0.set_ylim(0,2) # data

axes0.set_xlim(4.25*1e9,4.40*1e9);
plot(w_s/2/pi,abs(Tr_sz[t]))









    



-141.632653061






    Out[96]:





[<matplotlib.lines.Line2D at 0x126385cf8>]



In [46]:

    
###A1
t = 31
print(P_cavity[t])
File = A1

fig, axes = subplots(1,1, figsize=(16,6),sharex=True)

# fig.subplots_adjust(left=0, wspace=0.25

axes.plot(1e9*w_s/2/pi+4e6,Tr_sz[t],'k-',linewidth=3)
axes.set_xlabel(r'f (GHz)', fontsize=18)
axes.set_ylabel(r'$Tr[\rho \sigma_z]$ or $Tr[\rho a^{\dagger}a]$ ', fontsize=18)
axes.set_title(r'$Tr[\rho \sigma_z]$',fontsize=20)
#####
axes0 = axes.twinx()
axes0.plot(F,(File)*1e6,'ro', label='data' )

# axes.legend(loc=2)
#left plot
axes.set_ylim(0.25,1.35) # simulation 
axes0.set_ylim(0,3) # data

axes0.set_xlim(4.25*1e9,4.40*1e9);
plot(w_s/2/pi,abs(Tr_sz[t]))









    



-141.020408163






    



/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/numpy/core/numeric.py:462: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)






    Out[46]:





[<matplotlib.lines.Line2D at 0x11d74afd0>]



In [47]:

    
###A2
t = 43
print(P_cavity[t])
File = A2

fig, axes = subplots(1,1, figsize=(16,6),sharex=True)

# fig.subplots_adjust(left=0, wspace=0.25

axes.plot(1e9*w_s/2/pi+15e6,Tr_sz[t],'k-',linewidth=3)
axes.set_xlabel(r'f (GHz)', fontsize=18)
axes.set_ylabel(r'$Tr[\rho \sigma_z]$ or $Tr[\rho a^{\dagger}a]$ ', fontsize=18)
axes.set_title(r'$Tr[\rho \sigma_z]$',fontsize=20)
#####
axes0 = axes.twinx()
axes0.plot(F,(File)*1e6,'ro', label='data' )

# axes.legend(loc=2)
#left plot
axes.set_ylim(0.65,1.02) # simulation 
axes0.set_ylim(0,5) # data

axes0.set_xlim(4.25*1e9,4.40*1e9);
plot(w_s/2/pi,abs(Tr_sz[t]))









    



-133.673469388






    Out[47]:





[<matplotlib.lines.Line2D at 0x11dd1b5f8>]



In [48]:

    
###A3
t = 44
print(P_cavity[t])
File = A3

fig, axes = subplots(1,1, figsize=(16,6),sharex=True)

# fig.subplots_adjust(left=0, wspace=0.25

axes.plot(1e9*w_s/2/pi+15e6,Tr_sz[t],'k-',linewidth=3)
axes.set_xlabel(r'f (GHz)', fontsize=18)
axes.set_ylabel(r'$Tr[\rho \sigma_z]$ or $Tr[\rho a^{\dagger}a]$ ', fontsize=18)
axes.set_title(r'$Tr[\rho \sigma_z]$',fontsize=20)
#####
axes0 = axes.twinx()
axes0.plot(F,(File)*1e6,'ro', label='data' )

# axes.legend(loc=2)
#left plot
axes.set_ylim(0.7,1.0) # simulation 
axes0.set_ylim(0,7) # data

axes0.set_xlim(4.25*1e9,4.40*1e9);
plot(w_s/2/pi,abs(Tr_sz[t]))









    



-133.06122449






    Out[48]:





[<matplotlib.lines.Line2D at 0x10f4f6d68>]



In [49]:

    
###A4
t = 45
print(P_cavity[t])
File = A4

fig, axes = subplots(1,1, figsize=(16,6),sharex=True)

# fig.subplots_adjust(left=0, wspace=0.25

axes.plot(1e9*w_s/2/pi+13.7e6,Tr_sz[t],'k-',linewidth=3)
axes.set_xlabel(r'f (GHz)', fontsize=18)
axes.set_ylabel(r'$Tr[\rho \sigma_z]$ or $Tr[\rho a^{\dagger}a]$ ', fontsize=18)
axes.set_title(r'$Tr[\rho \sigma_z]$',fontsize=20)
#####
axes0 = axes.twinx()
axes0.plot(F,(File)*1e6,'ro', label='data' )

# axes.legend(loc=2)
#left plot
axes.set_ylim(0.7,1.03) # simulation 
axes0.set_ylim(0,9) # data

axes0.set_xlim(4.25*1e9,4.40*1e9);
plot(w_s/2/pi,abs(Tr_sz[t]))









    



-132.448979592






    Out[49]:





[<matplotlib.lines.Line2D at 0x11d636940>]



In [50]:

    
###A5
t = 46
print(P_cavity[t])
File = A5

fig, axes = subplots(1,1, figsize=(16,6),sharex=True)

# fig.subplots_adjust(left=0, wspace=0.25

axes.plot(1e9*w_s/2/pi+12e6,Tr_sz[t],'k-',linewidth=3)
axes.set_xlabel(r'f (GHz)', fontsize=18)
axes.set_ylabel(r'$Tr[\rho \sigma_z]$ or $Tr[\rho a^{\dagger}a]$ ', fontsize=18)
axes.set_title(r'$Tr[\rho \sigma_z]$',fontsize=20)
#####
axes0 = axes.twinx()
axes0.plot(F,(File)*1e6,'ro', label='data' )

# axes.legend(loc=2)
#left plot
axes.set_ylim(0.74,1.01) # simulation
axes0.set_ylim(0,10) # data

axes0.set_xlim(4.25*1e9,4.40*1e9);
plot(w_s/2/pi,abs(Tr_sz[t]))









    



-131.836734694






    Out[50]:





[<matplotlib.lines.Line2D at 0x11e2501d0>]

Fix cavity tone spectroscopy



In [263]:

    
P_qubit = -140 # Input power qubit (0-20-20-76 - cavity(?))

P_cavity = -(106 + 30) # Input power cavity

Gamma_d = sqrt(Q_c*kappa_n*1e9*dBmtoW(P_cavity)/(2*Q_l*hbar*w_r_tilde*1e9))/1e9  # Drive amplitude qubit
print('Drive amplitude cavity =',Gamma_d/2/pi, 'GHz')

Gamma_s = sqrt(Q_c*kappa_n*1e9*dBmtoW(P_qubit)/(2*Q_l*hbar*w_ge_tilde*1e9))/1e9 # Drive amplitude cavity ????
print('Drive amplitude qubit =',Gamma_s/2/pi,'GHz')


N = 20 # number of states cavity
wr = (5.0596)  *2*pi
Delta_d_tilde = w_r_tilde - wr
w_s = 2 * pi * linspace(f_s_i,f_s_f,500)
amp1 = []
amp2 = []
for wq in w_s:

    Delta_s_tilde = w_ge_tilde - wq
    rho, temp, temp2 = calc_spectrum(N,tlist,Delta_d_tilde,Delta_s_tilde,chi,Gamma_d,Gamma_s, zeta_1,zeta_2,kappa_n,gamma_rel,gamma_dep,n_th_a)
    amp1.append((temp.tr()))
    amp2.append((temp2.tr()))
    
fig, axes = plt.subplots(1, 2, figsize=(16,6))

axes[0].plot(w_s/2/pi,amp1)
axes[0].set_ylabel('Cavity Pump frequency (GHz)')
axes[0].set_xlabel('Qubit Pump frequency (GHz)')


axes[1].plot(w_s/2/pi,amp2)
axes[1].set_ylabel('Cavity Pump frequency (GHz)')
axes[1].set_xlabel('Qubit Pump frequency (GHz)')
# axes[1].title()

#___________________________________________________#######

rho.ptrace(0).diag()

fig, axes = plt.subplots(1, 2, figsize=(16,6))

xvec = np.linspace(-5,5,200)

rho_cavity = rho.ptrace(0)
W = wigner(rho_cavity, xvec, xvec)
wlim = abs(W).max()
axes[1].contourf(xvec, xvec, W, 100, norm=mpl.colors.Normalize(-wlim,wlim), cmap=plt.get_cmap('RdBu'))
axes[1].set_xlabel(r'Im $\alpha$', fontsize=18)
axes[1].set_ylabel(r'Re $\alpha$', fontsize=18)

axes[0].bar(arange(0, N), real(rho_cavity.diag()), color="blue", alpha=0.6)
axes[0].set_xlabel(r'$n$', fontsize=18)
axes[0].set_ylabel(r'Occupation probability', fontsize=18)
axes[0].set_ylim(0, 1)
axes[0].set_xlim(0, N);









    



Drive amplitude cavity = 0.000387761810913 GHz
Drive amplitude qubit = 0.0002635497665 GHz






    



---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-263-63c3e7fd8036> in <module>()
     19 
     20     Delta_s_tilde = w_ge_tilde - wq
---> 21     rho, temp, temp2 = calc_spectrum(N,tlist,Delta_d_tilde,Delta_s_tilde,chi,Gamma_d,Gamma_s, zeta_1,zeta_2,kappa_n,gamma_rel,gamma_dep,n_th_a)
     22     amp1.append((temp.tr()))
     23     amp2.append((temp2.tr()))

/Users/rouxinol/Python/Modules/MyQubit.py in calc_spectrum(N, tlist, Delta_d_tilde, Delta_s_tilde, chi, Gamma_d, Gamma_s, zeta_1, zeta_2, kappa_n, gamma_rel, gamma_dep, n_th_a)
     76 #     spec = []
     77 
---> 78     rho    = steadystate(H,c_op_list)
     79     rho_sz = rho * sz
     80     rho_n  = rho * (a.dag()*a + I/2)

/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/qutip/steadystate.py in steadystate(A, c_op_list, **kwargs)
    213     if ss_args['method'] == 'direct':
    214         if ss_args['sparse']:
--> 215             return _steadystate_direct_sparse(A, ss_args)
    216         else:
    217             return _steadystate_direct_dense(A, ss_args)

/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/qutip/steadystate.py in _steadystate_direct_sparse(L, ss_args)
    347         lu = splu(L, permc_spec=ss_args['permc_spec'],
    348                   diag_pivot_thresh=ss_args['diag_pivot_thresh'],
--> 349                   options=dict(ILU_MILU=ss_args['ILU_MILU']))
    350         v = lu.solve(b)
    351         _direct_end = time.time()

/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/scipy/sparse/linalg/dsolve/linsolve.py in splu(A, permc_spec, diag_pivot_thresh, drop_tol, relax, panel_size, options)
    240         _options.update(options)
    241     return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
--> 242                           ilu=False, options=_options)
    243 
    244 

KeyboardInterrupt:



In [628]:

    
# 
fig, axes = subplots(1,2, figsize=(16,6),sharex=True)
# a = 4200
# b = 4400
File = A5

fig.subplots_adjust(left=0, wspace=0.25)


axes0 = axes[0].twinx()
axes0.plot(F,(File)*1e6,'r*', label='data' )
axes0.set_ylabel(r'$\mu V$', fontsize=18)
axes0.legend(loc=0)

axes[0].plot(1e9*w_s/2/pi+4e6, (amp1),'k-',linewidth=2,label='Simulation')
axes[0].set_xlabel(r'f (GHz)', fontsize=18)
axes[0].set_ylabel(r'$Tr[\rho \sigma_z]$ or $Tr[\rho a^{\dagger}a]$ ', fontsize=18)
axes[0].set_title(r'$Tr[\rho \sigma_z]$',fontsize=20)
axes[0].legend(loc=2)

# Graphic right


axes[1].plot(1e9*w_s/2/pi+4e6, (amp2),'k-',linewidth=2, label='Simulation')
axes[1].set_xlabel(r'f (GHz)', fontsize=18)
axes[1].set_ylabel(r'$Tr[\rho \sigma_z]$ or $Tr[\rho a^{\dagger}a]$ ', fontsize=18)
axes[1].set_title(r'$Tr[\rho {a^{\dagger}a}]$',fontsize=20)
axes[1].legend(loc=0)

axes1 = axes[1].twinx()
axes1.plot(F,(File)*1e6,'r*', label='data')
axes1.set_ylabel(r'$\mu V$', fontsize=18)
axes1.legend(loc=2)


#left plot
axes[0].set_ylim(0.5,1.08) # simulation 
axes0.set_ylim(2,9) # data

# right plot
axes[1].set_ylim(3.9,5.35) #simulation
axes1.set_ylim(2,9) # Data

axes1.set_xlim(4.25*1e9,4.40*1e9);

Fix cavity tone with 10kHz spreed



In [36]:

    
P_qubit = -140 # Input power qubit (0-20-20-76 - cavity(?))

P_cavity = -(106 + 26) # Input power cavity

Gamma_d = sqrt(Q_c*kappa_n*1e9*dBmtoW(P_cavity)/(2*Q_l*hbar*w_r_tilde*1e9))/1e9  # Drive amplitude qubit
print('Drive amplitude cavity =',Gamma_d/2/pi, 'GHz')

Gamma_s = sqrt(Q_c*kappa_n*1e9*dBmtoW(P_qubit)/(2*Q_l*hbar*w_ge_tilde*1e9))/1e9 # Drive amplitude cavity ????
print('Drive amplitude qubit =',Gamma_s/2/pi,'GHz')


N = 20 # number of states cavity

f_center = 5.0597
Delta_f  = 0.0002
f_d_i = f_center - Delta_f # cavity initial frequency 
f_d_f = f_center + Delta_f# cavity final frequency
steps_d = 6

w_d = 2 * pi * linspace(f_d_i,f_d_f,steps_d)

f_s_i = 4.2 # qubit initial freq
f_s_f = 4.4  # qubit final freq
steps_s = 500
w_s = 2 * pi * linspace(f_s_i,f_s_f,steps_s)

amp_1a = []
amp_1b = []
for wr in w_d:
    
    for wq in w_s:
        
        Delta_d_tilde = w_r_tilde - wr
        Delta_s_tilde = w_ge_tilde - wq
        rho,temp, temp2 = calc_spectrum(N,tlist,Delta_d_tilde,Delta_s_tilde,chi,Gamma_d,Gamma_s, zeta_1,zeta_2,kappa_n,gamma_rel,gamma_dep,n_th_a)
        amp_1a.append((temp.tr()))
        amp_1b.append((temp2.tr()))

Tr_sz = reshape(amp_1a,(-1,len(w_s+1)))
Tr_n = reshape(amp_1b,(-1,len(w_s+1)))









    



Drive amplitude cavity = 0.000614561054412 GHz
Drive amplitude qubit = 0.0002635497665 GHz



In [37]:

    
fig, ax = subplots(1,2, figsize=(16,10))
im = ax[0].pcolor(w_s/2/pi, w_d/2/pi ,(abs(Tr_sz)))#,vmin=0, vmax=1)
ax[0].set_xlim(f_s_i,f_s_f)
ax[0].set_ylim(f_d_i,f_d_f)
ax[0].set_ylabel(r'Cavity tone frequency (GHz)')
ax[0].set_xlabel(r'Qubit tone frequency (GHz)')

fig.colorbar(im, ax=ax[0])









    Out[37]:





<matplotlib.colorbar.Colorbar at 0x117ba17f0>



In [71]:

    
mean_sz = []
mean_sz = Tr_sz.mean(axis=0,)
Tr_sz.shape









    Out[71]:





(6, 500)



In [39]:

    
fig, axes = plt.subplots(1, 1, figsize=(16,6))
for i in range(steps_d):
    axes.plot(w_s/2/pi,Tr_sz[i])
axes.plot(w_s/2/pi,mean_sz,'ko')









    Out[39]:





[<matplotlib.lines.Line2D at 0x1177d8e10>]



In [40]:

    
# 
fig, axes = subplots(1,1, figsize=(16,6),sharex=True)
# a = 4200
# b = 4400
File = A3

fig.subplots_adjust(left=0, wspace=0.25)


axes0 = axes.twinx()
axes0.plot(F,(File)*1e6,'r*', label='data' )
# axes0.set_ylabel(r'$\mu V$', fontsize=18)
# axes0.legend(loc=0)

axes.plot(1e9*w_s/2/pi+4e6, (mean_sz),'k--',linewidth=2,label='Simulation')
for i in range(steps_d):
    axes.plot(1e9*w_s/2/pi+3e6,Tr_sz[i])
    
axes.set_xlabel(r'f (GHz)', fontsize=18)
axes.set_ylabel(r'$Tr[\rho \sigma_z]$ or $Tr[\rho a^{\dagger}a]$ ', fontsize=18)
axes.set_title(r'$Tr[\rho \sigma_z]$',fontsize=20)
# axes.legend(loc=2)

# Graphic right


#left plot
axes.set_ylim(0,1.15) # simulation 
axes0.set_ylim(0,9) # data


axes0.set_xlim(4.25*1e9,4.40*1e9);



In [112]:

    
mean_sz = (Tr_sz[4]+ Tr_sz[5] )/2#+ Tr_sz[2]+Tr_sz[3]+Tr_sz[4] +Tr_sz[5])/6



In [113]:

    
# 
fig, axes = subplots(1,1, figsize=(16,6),sharex=True)
# a = 4200
# b = 4400
File = A2

fig.subplots_adjust(left=0, wspace=0.25)


axes0 = axes.twinx()
axes0.plot(F,(File)*1e6,'r*', label='data' )
# axes0.set_ylabel(r'$\mu V$', fontsize=18)
# axes0.legend(loc=0)

axes.plot(1e9*w_s/2/pi+4e6, (mean_sz),'k--',linewidth=2,label='Simulation')
# for i in range(steps_d):
#     axes.plot(1e9*w_s/2/pi+3e6,Tr_sz[i])
    
axes.set_xlabel(r'f (GHz)', fontsize=18)
axes.set_ylabel(r'$Tr[\rho \sigma_z]$ or $Tr[\rho a^{\dagger}a]$ ', fontsize=18)
axes.set_title(r'$Tr[\rho \sigma_z]$',fontsize=20)
# axes.legend(loc=2)

# Graphic right


#left plot
axes.set_ylim(0.65,1.1) # simulation 
axes0.set_ylim(0,5) # data


axes0.set_xlim(4.25*1e9,4.40*1e9);

New setup hamiltonian



In [102]:

    
N = 20 # number of states cavity


P_cavity = -(106 + 32) # Input power cavity

Gamma_d = sqrt(Q_c*kappa_n*1e9*dBmtoW(P_cavity)/(2*Q_l*hbar*w_r_tilde*1e9))/1e9  # Drive amplitude qubit
print('Drive amplitude cavity =',Gamma_d/2/pi, 'GHz')

# w_s = 2 * pi * linspace(f_s_i,f_s_f,steps_s)

wr = (5.0596)  * 2* pi
Delta_d_tilde = w_r_tilde - wr 

# qubit freq scan
f_s_i = 4.2 # qubit initial freq
f_s_f = 4.4  # qubit final freq
steps_s = 200
w_s = 2 * pi * linspace(f_s_i,f_s_f,steps_s)

# Power scan
P_i = -110 # qubit tone start
P_f = -160 # qubit tone final
steps_p = 50 # 
P_qubit = linspace(P_i,P_f, steps_p)




amp_sz = []
amp_n = []
amp_x = []
amp_a = []
for Pq in P_qubit:
    
    for wq in w_s:
        
        
        Delta_s_tilde = w_ge_tilde - wq
        Gamma_s = sqrt(Q_c*kappa_n*1e9*dBmtoW(Pq)/(2*Q_l*hbar*w_ge_tilde*1e9))/1e9

        rho,temp, temp2, temp3, temp4 = calc_spectrum_2(N,Delta,w_q,Delta_d_tilde,Delta_s_tilde,chi,Gamma_d,Gamma_s, zeta_1,zeta_2,kappa_n,gamma_rel,gamma_dep, n_th_a)
        amp_sz.append((temp.tr()))
        amp_n.append((temp2.tr()))
        amp_x.append((temp3.tr()))
        amp_a.append((temp4.tr()))









    



Drive amplitude cavity = 0.000308010154756 GHz






    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-102-aa2870cf6f00> in <module>()
     39         Gamma_s = sqrt(Q_c*kappa_n*1e9*dBmtoW(Pq)/(2*Q_l*hbar*w_ge_tilde*1e9))/1e9
     40 
---> 41         rho,temp, temp2, temp3, temp4 = calc_spectrum_2(N,Delta,w_q,Delta_d_tilde,Delta_s_tilde,chi,Gamma_d,Gamma_s, zeta_1,zeta_2,kappa_n,gamma_rel,gamma_dep, n_th_a)
     42         amp_sz.append((temp.tr()))
     43         amp_n.append((temp2.tr()))

NameError: name 'calc_spectrum_2' is not defined



In [ ]: