Note
This page was generated from tut//4Analysis//4.34TransmonqubitCPBhamiltonianchargebasis.ipynb.
Modeling transmon qubit Cooperpair box Hamiltonian in the charge basis¶
(Zlatko Minev, Christopher Warren, Nick Lanzillo 2021)
This module models the transmon qubit in the cooperpair charge basis, assuming wrapped junction phase variable. The Hamiltonian is given by:
where \(E_{C}\) is the charging energy, \(E_{J}\) is the Josephson energy, \(\hat n\) is the number of Cooper pairs transferred between charge islands, \(\hat{\phi}\) is the gaugeinvariant phase difference between charge islands, and \(n_{g}\) is effective offset charge of the device. Expressions for the charging energy, Josephson energy and offset charge can be written as:
where \(C_{\Sigma} = C_{J}+C_{B}+C_{g}\) (the sum of the Josephson capacitance, shunting capacitance and gate capacitance), \(L_{J}\) is the inductance of the Josephson junction, and \(\phi\) is the magnetic flux.
The variables are
Observe that \(\hat \phi\) and \(\hat n\) are both dimensiuonless, and they obey the commutation relationship:
The Hamiltonian can be written in the charge (\(\hat n\)) basis as:
Where \(\hat{n} = \sum_{n=\inf}^{\inf} n\rangle\langle n\)
Hcpb class¶
Hamiltonianmodel Cooper pair box (Hcpb) class.
Used to model analytically the CPB Hamiltonian quickly and efficiently. Solves in charge basis tridiagonal eigenvalue problem for arbitrary Ej, Ec, ng values.
As long as nlevels remains fixed the number of charge states considered does not change and it does not recreate the arrays, just recomputes the properties
Returns all properties of interest for the CPB.
This model is closer to the analytic solution than the Duffing oscillator model. Can work backwards from target qubit parameters to get the Ej, Ec or use input Ej, Ec to find the spectrum of the Cooper Pair Box.
@author: Christopher Warren (Chalmers University of Technology), updated by Zlatko K. Minev (IBM Quantum)
@date: 2020, 2021
Let’s start by importing the key modules for this demo:
[2]:
%load_ext autoreload
%autoreload 2
%config IPCompleter.greedy = True
%matplotlib inline
%config Completer.use_jedi = False
%config InlineBackend.figure_format = 'svg'
from qiskit_metal.analyses.hamiltonian.transmon_charge_basis import Hcpb
from qiskit_metal.analyses.hamiltonian.transmon_CPB_analytic import Hcpb_analytic
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
The autoreload extension is already loaded. To reload it, use:
%reload_ext autoreload
Let’s model a transmon¶
Energy levels¶
We can get a feel for how to use the Hcpb class by plotting the first few eigenvalues as a function of offset charge, similar to the plots reported in Phys. Rev. A 76, 042319 (2007.) Let’s start by defining the range of offset charge from 2.0 to +2.0 and also by defining a normalization for the eigenvalues, which will be the transition energy between the first two states evaluated at the degenercy point where ng=0.5. For this exercise, we’ll take Josephson Energy to be equal to the charging energy:
[3]:
x = np.linspace(2.0,2.0,101) # this represents the charging energy (ng)
H_norm = Hcpb(nlevels=2, Ej=1000.0, Ec=1000.0, ng=0.5) # Hamiltonian definition
norm = H_norm.fij(0,1) # normalization constant
Next, we’ll sweep over the offset charge and calculate the first three eigenvalues for a given value of ng. We’ll need to define a new Hamiltonian for this.
[4]:
E0, E1, E2 = [], [], []
# For a given value of offset charge (ng, represented by x) we will calculate the CPB Hamiltonian using the previously assigned values of E_J and E_C. Then we calculate the eigenvalue for a given value of m.
for i in x:
H = Hcpb(nlevels=3, Ej=1000.0, Ec=1000.0, ng=i)
E0.append(H.evalue_k(0)/norm)
E1.append(H.evalue_k(1)/norm)
E2.append(H.evalue_k(2)/norm)
# define the minimum of E0 and set this to E=0
floor = min(E0)
plt.plot(x, E0  floor, 'k', label="m=0")
plt.plot(x, E1  floor, 'r', label="m=1")
plt.plot(x, E2  floor, 'b', label="m=2")
plt.xlabel("ng")
plt.ylabel("Em/E01")
plt.legend(title="Energy Level:", loc='upper right')
[4]:
<matplotlib.legend.Legend at 0x1c1fdccfbc8>
Comparing with the Analytic Expressions for Energy¶
We can compare calculated eigenvalues with the analytic solutions by using the “Hcpb_analytic” class, which calculates the transmon eigenvalues analytically using Mathieu characteristic values instead of a matrixbased approach. Let’s compare the calculated values of the lowest energy at zero offset charge in both cases:
[5]:
# this is using the Hcpb approach as above, solving the charge basis tridiagonal eigenvalue problem:
H_CPB = Hcpb(nlevels=15, Ej=13971.3, Ec=295.2, ng=0.0)
# this using the Hcpb_analytic class, which solves using the exact (analytic) solutions in terms of Mathieu characteristic values:
H_CPB_analytic = Hcpb_analytic(Ej=13971.3, Ec=295.2, ng=0.0)
# print and compare energies
print("E0 (HCPB):", H_CPB.evalue_k(0))
print("E0 (HCPB analytic):", H_CPB_analytic.evalue_k(0))
print("Error:", 100*(H_CPB_analytic.evalue_k(0)  H_CPB.evalue_k(0)) / H_CPB_analytic.evalue_k(0))
E0 (HCPB): 11175.114908534233
E0 (HCPB analytic): 11175.114908534231
Error: 1.6277142726798518e14
As we can see above, the calculated eigenvalues using the Hcbp class match the analytic values extremely well.
Wavefunctions¶
Let’s define a new Hamiltonian, this time with \(E_{J} >> E_{C}\) and an offset charge of \(ng=0.001\). We can calculate the transition energy between the lowest two states as well as the anharmonicity with the following:
[6]:
H = Hcpb(nlevels=3, Ej=13971.3, Ec=295.2, ng=0.001)
print(f"""
Transmon frequencies
ω01/2π = {H.fij(0,1): 6.0f} MHz
α/2π = {H.anharm(): 6.0f} MHz
""")
Transmon frequencies
ω01/2π = 5604 MHz
α/2π = 11 MHz
Note that both the transition energy and the anharmonicity are read out in units of Megahertz (MHz.)
We can plot the eigenstates (wavefunctions) of the transmon qubit using the commands below:
[7]:
import matplotlib.pyplot as plt
for k in range (3):
ψ, θ = H.psi_k(k,101)
plt.plot(θ, ψ.real+ψ.imag, label=f"{k}>") # it's in either quadrature, but not both
plt.xlabel("Junction phase θ (wrapped in the interval [π, π])")
plt.ylabel("Re(ψ(θ))")
plt.legend(title="Level")
[7]:
<matplotlib.legend.Legend at 0x1c1fdd969c8>
Verifying Orthonormality of the Wavefunctions¶
Here, we can verify the orthonormality of the wavefunctions. Let’s first take the first two eigenstates and verify that their inner product is zero, thereby confirming orthogonality. Note that since the wavefunctions can be complex, we need to take the complex conjugate of \(\Psi1\).
[8]:
Psi0, theta0 = H.psi_k(0)
Psi1, theta1 = H.psi_k(1)
print(np.dot(Psi0,Psi1.conj()))
(2.2354805733859444e093.4875983780408155e17j)
Indeed, we see that the dot product is essentially zero (within numerical precision.) Next, let’s take the inner product of the first eigenstate with itself, checking that we get an output of unity:
[9]:
print(np.dot(Psi0, Psi0.conj()))
print(np.dot(Psi1, Psi1.conj()))
(1+0j)
(1+0j)
Indeed we see that the dot products are essentially equal to unity, confirming that the states are appropriately normalized.
Additional Analysis: Charge Dispersion, Energy Level Differences, Anharmonicity and Dephasing Time (T2)¶
Charge Dispersion¶
The peaktopeak value of the charge dispersion for the mth energy level is given by the expression: \(\epsilon_{m} = E_{m}(n_{g}=0.5)  E_{m}(n_{g}=0.0)\). We can plot \(\epsilon_{m}/E_{01}\) as a function of \(E_{J}/E_{C}\) for the first few energy levels and reproduce the figure published in Phys. Rev. A 76, 042319 (2007) (Figure 4(a)).
We can start by defining a value of charging energy and creating empty lists for \(\epsilon_{0}\) through \(\epsilon_{4}\):
[10]:
E_c=100.0 # charging energy
epsilon0, epsilon1, epsilon2, epsilon3 = [], [], [], [] # charge dispersion for m=0 through m=4
x = np.linspace(1,140,101) # this this ratio of Ej/Ec which will go on the xaxis.
Next, we simply evaluate the expression given above for \(\epsilon_{m}\) based on \(E_{m}\) and \(E_{0}\). We use two separate Hamiltonians to do this; one evaluated at \(n_{g}=0.5\) and one evaluated at \(n_{g}=0.0\). We also normalize by the transition energy between the lowest two states evaluated at the degeneracy point (\(E_{01}\).) Finally, we populate the lists each \(\epsilon_{m}\).
[11]:
for i in x:
E_j = i*E_c
H_zero = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.0)
H_half = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)
H_norm = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)
norm = H_norm.fij(0,1) # normalization constant
epsilon0.append(abs(H_half.evalue_k(0)  H_zero.evalue_k(0))/norm)
epsilon1.append(abs(H_half.evalue_k(1)  H_zero.evalue_k(1))/norm)
epsilon2.append(abs(H_half.evalue_k(2)  H_zero.evalue_k(2))/norm)
epsilon3.append(abs(H_half.evalue_k(3)  H_zero.evalue_k(3))/norm)
We can plot these values to see the exponential decrease in the charge dispersion with increasing \(E_{J}/E_{C}\):
[12]:
plt.plot(x, epsilon0, 'k', label="m=0")
plt.plot(x, epsilon1, 'b', label="m=1")
plt.plot(x, epsilon2, 'r', label="m=2")
plt.plot(x, epsilon3, 'g', label="m=3")
plt.yscale("log")
plt.xlabel("EJ/Ec")
plt.ylabel("Epsilon/E01")
plt.legend(title="Energy Level", loc='upper right')
[12]:
<matplotlib.legend.Legend at 0x1c1fde55f48>
Energy Level Differences¶
We can also evaluate the energy level difference (\(E_{m0} = E_{m}  E_{0}\)) evaluated at the degeneracy point (\(n_{g}=0.5\)) and compare to Fig. 4(b) of Phys. Rev. A 76, 042319 (2007). To do this, we just need to create empty lists for the energy levels (\(E_{0}\) through \(E_{3}\)) as well as the energy level differences (\(E_{00}\) through \(E_{30}\).)
[13]:
E00, E10, E20, E30 = [], [], [], []
E0, E1, E2, E3 = [], [], [], []
Next, we’ll sweep over \(E_{J}/E_{C}\) from 0 to 140 (using the variable x defined above) and at each point, we’ll construct the Hamiltonian and take the difference in eigenvalues evaluated at the degeneracy point \(n_{g}=0.5\). We’ll also normalize the energy level differences by the charging energy to be consistent with Fig. 4(b) in the above reference:
[14]:
for i in x:
H = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)
E0 = H.evalue_k(0)
E1 = H.evalue_k(1)
E2 = H.evalue_k(2)
E3 = H.evalue_k(3)
E00.append((E0  E0)/E_c)
E10.append((E1  E0)/E_c)
E20.append((E2  E0)/E_c)
E30.append((E3  E0)/E_c)
We can plot these results and see how the energy level differences increase with increasing \(E_{J}/E_{C}\) ratio:
[15]:
plt.plot(x,E00,'k',label="m0")
plt.plot(x,E10,'b',label="m=1")
plt.plot(x,E20,'g',label="m=2")
plt.plot(x,E30,'r',label="m=3")
plt.xlabel("Ej/Ec")
plt.ylabel("Em0/Ec")
plt.legend(title="Energy Levels", loc='upper right')
[15]:
<matplotlib.legend.Legend at 0x1c1fdf815c8>
Anharmonicity¶
We know that for the transmon qubit, having the Josephson Energy much larger than the charging energy (\(E_{J} >> E_{C}\)) results in a decrease in anharmonicity. The latter is critical for a functional qubit in which the energy difference between the lowest two states (\(E_{01}\)) is sufficiently different than the energy difference between the second and third states, \(E_{12}\). The absolute anharmonicity is defined as \(\alpha = E_{12}  E_{01}\), while the relative anharmonicity is defined as \(\alpha_{r} = \alpha/E_{01}\).
We can easily make a plot of the anharmonicity as a function of \(E_{J}/E_{C}\) using the Hcpb class. Let’s have the ratio of \(E_{J}/E_{C}\) (which we’ll call x) vary from 0 to 80, and then we’ll create empty lists for \(\alpha\) and \(\alpha_r\).
[16]:
x = np.linspace(0,80,101) #EJ/EC
alpha = []
alpha_r = []
Next, we’ll just sweep over x and at each value we’ll calculate both the absolute and relative anharmonicity.
[17]:
for i in x:
H_anharm = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)
alpha.append(H_anharm.anharm())
alpha_r.append(H_anharm.anharm()/H_anharm.fij(0,1))
WARNING:py.warnings:C:\Users\NICKLANZILLO\Anaconda3\envs\metal2\lib\sitepackages\ipykernel_launcher.py:4: RuntimeWarning: divide by zero encountered in double_scalars
after removing the cwd from sys.path.
Note we get a warning here because the relative anharmonicity blows up as Ej/Ec goes to zero. Then we can plot the results:
[18]:
plt.figure(1)
plt.subplot(131)
plt.plot(x,alpha)
plt.xlabel("Ej/Ec")
plt.ylabel("alpha")
plt.subplot(133)
plt.plot(x,alpha_r)
plt.ylim(0.2, 1.0)
plt.xlabel("Ej/Ec")
plt.ylabel("alpha_r")
[18]:
Text(0, 0.5, 'alpha_r')
Indeed we see that the anharmonicity decays with the inverse of \(E_{J}/E_{C}\) for small values of \(E_{J}/E_{C}\) before reaching a minimum just before \(E_{J}/E_{C} \approx 20.0\), then changing sign and approaching zero as \(E_{J}/E_{C}\) approaches infinity. This matches very closely to the results found in Figure 5 of Phys. Rev. A 76, 042319 (2007).
Dephasing Time (T2)¶
We can estimate the qubit dephasing time (T2) due to charge noise by the following expression: \(T_{2} = \frac{\hbar}{A\pi \epsilon_{1}}\) where \(A\) is on the order of \(1E4\) according to Phys. Rev. A 76, 042319 (2007). Since this is essentially just the inverse of the charge dispersion for \(\epsilon_{1}\), we can easily calculate T2 as a function of \(E_{J}/E_{C}\) with the following:
[19]:
x = np.linspace(0,80,101) # ratio of Ej/Ec, varying from 0 to 80
T2 = [] # empty list for T2
Next, we’ll just calculate the T2 time as a function of \(E_{J}/E_C{}\), noting that the output of the Hcpb eigenvalue caluclation is in units of E/h ~ MHz (1E6 Hz). So our calculation simplifies to: \(T_{2} = \frac{1.0}{2*(1E4)*(1E6) \epsilon_{1}}\) where \(\epsilon_{1}\) is the direct output of the Hcpb eigenvalue calculation:
[20]:
E_c = 1000.0
for i in x:
H_half = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)
H_zero = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.0)
eps = abs(H_half.evalue_k(1)  H_zero.evalue_k(1))
T2.append(1.0/(2.0*(1E4)*(1E6)*eps) )
plt.plot(x, T2)
plt.yscale("log")
plt.xlabel("Ej/Ec")
plt.ylabel("T2 (sec)")
[20]:
Text(0, 0.5, 'T2 (sec)')
Indeed, we see that the dephasing time increases exponentially with increasing \(E_{J}/E_{C}\), which is one of the critical features of the transmon qubit. By increasing the \(E_{J}/E_{C}\) ratio, we reduce sensitivity to charge noise without sacrificing too much anharmonicity, resulting in greatly improved dephasing time. This plot matches very closely to Fig. 5(c) in Phys. Rev. A 76, 042319 (2007).
Qutip simulation¶
Wrapper around Qutip to output the diagonalized Hamiltonian truncated up to n levels of the transmon for modeling
[21]:
H.h0_to_qutip(6)
[21]:
Wrapper around Qutip to output the number operator (charge) for the Transmon Hamiltonian in the energy eigenbasis. Used for computing the coupling between other elements in the system.
[22]:
H.n_to_qutip(6)
[22]:
We can use these matrices to perform some additional analysis. As an example, we can plot the fluctuation in the number of Cooper pairs as a function of \(E_{J} / E_{C}\) for the first three energy levels. This requires calculating the variance in the number of Cooper pairs, defined as: \(\sqrt{<n^2>  <n>^2 }\) where \(n\) is calculated for a given energy level m.
So let’s first set our xaxis (representing \(E_{J}/E_{C}\)) to run from 1 to 120 and we’ll create empty lists to represent the variances for the m=0, m=1, and m=2 states. We’ll also define a value of \(E_{C}\).
[23]:
x = np.linspace(1,120, 101)
E_c=350.0
var0 = []
var1 = []
var2 = []
Next, we’ll define the Hamiltonian matrix and calculate both \(<n>\) and \(<n^2>\). Then we’ll take the square root of the difference between the two for each value of \(E_{J}/E_{C}\). Lastly, these values get appended to the empty lists that we defined above:
[24]:
for i in x:
H = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)
n_squared = H.n_to_qutip(6)*H.n_to_qutip(6)
n = H.n_to_qutip(6)
fluc0 = np.sqrt(np.real(n_squared[0,0]  (n[0,0])*(n[0,0])))
fluc1 = np.sqrt(np.real(n_squared[1,1]  (n[1,1])*(n[1,1])))
fluc2 = np.sqrt(np.real(n_squared[2,2]  (n[2,2])*(n[2,2])))
var0.append(fluc0)
var1.append(fluc1)
var2.append(fluc2)
We can plot these and check that the resulting curves match quite closely to those found in Fig. 6 of Koch et al. (2007) PRA:
[25]:
plt.plot(x,var0,'k', label="m=0")
plt.plot(x,var1,'r', label="m=1")
plt.plot(x,var2,'b', label="m=2")
plt.xlabel("Ej/Ec")
plt.ylabel("Variance")
plt.legend(title="Energy Levels", loc='upper right')
[25]:
<matplotlib.legend.Legend at 0x1c1fe01f8c8>
Physically, the above plot demonstrates that the variations in the number of Cooper pairs is on the order of 12 for gorund and first excited states of the CPB operated in the transmon regime.
We can also investigate how the offdiagonal matrix elements of the CPB Hamiltonian evolve as a function of the \(E_{J}/E_{C}\) ratio. Let’s create four empty lists (\(m_{10}\), \(m_{20}\), \(m_{30}\) and \(m_{21}\)) where \(m_{ij} = <i  n  j>\). Then we can loop over the ratio of \(E_{J}/E_{C}\) and at each value, we calculate the offdiagonal matrix element:
[28]:
m10 = []
m20 = []
m30 = []
m21 = []
# For a given value of offset charge (ng, represented by x) we will calculate the CPB Hamiltonian using the previously assigned values of E_J and E_C. Then we calculate the eigenvalue for a given value of m.
for i in x:
H = Hcpb(nlevels=15, Ej=i*E_c, Ec=E_c, ng=0.5)
m10.append(np.real(H.n_to_qutip(6)[1,0]))
m20.append(np.real(H.n_to_qutip(6)[2,0]))
m30.append(np.real(H.n_to_qutip(6)[3,0]))
m21.append(np.real(H.n_to_qutip(6)[2,1]))
Plotting these offdiagonal matrix elements is then straightforward:
[29]:
plt.plot(x,m10,'k',label="<1n0>")
plt.plot(x,m20,'r',label="<2n0>")
plt.plot(x,m30,'y',label="<3n0>")
plt.plot(x,m21,'b',label="<2n1>")
plt.xlabel("Ej/Ec")
plt.ylabel("<inj>")
plt.legend(title="Energy Levels", loc='upper right')
[29]:
<matplotlib.legend.Legend at 0x1c1fe351ec8>
We see that the above plot closely matches Fig. 7 of Koch et al. (2007) PRA. Physically, we see that the coupling between adjacent transmon states (\(<1n0>\) and \(<2n1>\)) is the only significant coupling in the \(E_{J} >> E_{C}\) regime. The coupling between nonadjacent states (like m=0 and m=2) is essentially zero at all values of \(E_{J}/E_{C}\).
Experimental¶
Let’s use the “params_from_spectrum” function to calculate the target Ej and Ec values for a desired transmon frequency and anharmonicity. We’ll use values of Ej=13971.3 MHz and Ec=295.2 MHz.
[30]:
# 13971.3, Ec=295.2
ω, α = 5431, 341
EjEc = H.params_from_spectrum(ω, α) # set self.Ej, Cj
print(EjEc)
print("transmon frequency:", H.fij(0,1), "anharmonicity:", H.anharm())
[13952.97584426 295.61920264]
transmon frequency: 5431.003285965414 anharmonicity: 340.9980226844036
We can also calculate the value of Ej given the value of Ec and the transmon frequency:
[31]:
Ej = H.params_from_freq_fixEC(ω, 295.17)
print("Ej:", Ej)
Ej: 14424.680512350775
We see that this is close to the starting vaue of 13971.3 MHz that we began with.
New section on integrating sc_qubits.¶
This section is in development.
[ ]:
import scqubits as scq
[ ]:
qubit = scq.Transmon(
EJ=13.97124102543,
EC=0.295179,
ng=0.001,
ncut=40,
truncated_dim=4 # after diagonalization, we will keep 3 levels
)
evals = qubit.eigenvals(evals_count=12)
print(f"freq = {(evals[1]  evals[0])* 1000:.0f} MHz")
print(f"alpha = {((evals[2]  evals[1])  (evals[1]  evals[0]))* 1000:.0f} MHz")
[ ]:
qubit.plot_n_wavefunction()
qubit.plot_phi_wavefunction(which=[0,1,2], mode='real')
qubit.hamiltonian()
[ ]:
import numpy as np
ng_list = np.linspace(2, 2, 220)
qubit.plot_evals_vs_paramvals('ng', ng_list, evals_count=6, subtract_ground=False);
Transmon and the Oscillator¶
This section is currently in development.
For more information, review the Introduction to Quantum Computing and Quantum Hardware lectures below

Lecture Video  Lecture Notes  Lab 

Lecture Video  Lecture Notes  Lab 

Lecture Video  Lecture Notes  Lab 

Lecture Video  Lecture Notes  Lab 

Lecture Video  Lecture Notes  Lab 

Lecture Video  Lecture Notes  Lab 