Author: J. R. Johansson (robert@riken.jp), https://jrjohansson.github.io/
This lecture series was developed by J.R. Johannson. The original lecture notebooks are available here.
This is a slightly modified version of the lectures, to work with the current release of QuTiP. You can find these lectures as a part of the qutip-tutorials repository. This lecture and other tutorial notebooks are indexed at the QuTiP Tutorial webpage.
import matplotlib.pyplot as plt
import numpy as np
from qutip import Qobj, about, plot_energy_levels, ket2dm, mesolve
%matplotlib inline
The Hamiltonian for a Josephson charge qubit is
$\displaystyle H = \sum_n 4 E_C (n_g - n)^2 \left|n\right\rangle\left\langle n\right| - \frac{1}{2}E_J\sum_n\left(\left|n+1\right\rangle\left\langle n\right| + \left|n\right\rangle\left\langle n+1\right| \right)$
where $E_C$ is the charge energy, $E_J$ is the Josephson energy, and $\left| n\right\rangle$ is the charge state with $n$ Cooper-pairs on the island that makes up the charge qubit.
Below we will repeatedly need to obtain the charge qubit Hamiltonian for different parameters, and to plot the eigenenergies, so here we define two functions to do these tasks.
def hamiltonian(Ec, Ej, N, ng):
"""
Return the charge qubit hamiltonian as a Qobj instance.
"""
m = np.diag(4 * Ec * (np.arange(-N, N + 1) - ng) ** 2) + 0.5 * Ej * (
np.diag(-np.ones(2 * N), 1) + np.diag(-np.ones(2 * N), -1)
)
return Qobj(m)
def plot_energies(ng_vec, energies, ymax=(20, 3)):
"""
Plot energy levels as a function of bias parameter ng_vec.
"""
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
for n in range(len(energies[0, :])):
axes[0].plot(ng_vec, energies[:, n])
axes[0].set_ylim(-2, ymax[0])
axes[0].set_xlabel(r"$n_g$", fontsize=18)
axes[0].set_ylabel(r"$E_n$", fontsize=18)
for n in range(len(energies[0, :])):
axes[1].plot(
ng_vec,
(energies[:, n] - energies[:, 0]) /
(energies[:, 1] - energies[:, 0]),
)
axes[1].set_ylim(-0.1, ymax[1])
axes[1].set_xlabel(r"$n_g$", fontsize=18)
axes[1].set_ylabel(r"$(E_n-E_0)/(E_1-E_0)$", fontsize=18)
return fig, axes
def visualize_dynamics(result, ylabel):
"""
Plot the evolution of the expectation values stored in result.
"""
fig, ax = plt.subplots(figsize=(12, 5))
ax.plot(result.times, result.expect[0])
ax.set_ylabel(ylabel, fontsize=16)
ax.set_xlabel(r"$t$", fontsize=16);
N = 10
Ec = 1.0
Ej = 1.0
ng_vec = np.linspace(-4, 4, 200)
energies = np.array([hamiltonian(Ec, Ej, N, ng).eigenenergies()
for ng in ng_vec])
plot_energies(ng_vec, energies);
ng_vec = np.linspace(-1, 1, 200)
energies = np.array([hamiltonian(Ec, Ej, N, ng).eigenenergies()
for ng in ng_vec])
plot_energies(ng_vec, energies, ymax=(7.5, 3.0));
ng_vec = np.linspace(-4, 4, 200)
Ec = 1.0
Ej = 5.0
energies = np.array([hamiltonian(Ec, Ej, N, ng).eigenenergies()
for ng in ng_vec])
plot_energies(ng_vec, energies, ymax=(50, 3));
Ec = 1.0
Ej = 10.0
energies = np.array([hamiltonian(Ec, Ej, N, ng).eigenenergies()
for ng in ng_vec])
plot_energies(ng_vec, energies, ymax=(50, 3));
Ec = 1.0
Ej = 50.0
energies = np.array([hamiltonian(Ec, Ej, N, ng).eigenenergies()
for ng in ng_vec])
plot_energies(ng_vec, energies, ymax=(50, 3));
Note that the energy-level splitting is essentially independent of the gate bias $n_g$, at least for the lowest few states. This device insensitive to charge noise. But at the same time the two lowest energy states are no longer well separated from higher states (it has become more like an harmonic oscillator). But some anharmonicity still remains, and it can still be used as a qubit if the leakage of occupation probability of the higher states can be kept under control.
Let's go back to the charge regime, and look at the lowest few energy levels again:
N = 10
Ec = 1.0
Ej = 1.0
ng_vec = np.linspace(-1, 1, 200)
energies = np.array([hamiltonian(Ec, Ej, N, ng).eigenenergies()
for ng in ng_vec])
plot_energies(ng_vec, energies, ymax=(10, 3));
We can see that around $n_g = 0.5$ we have two lowest energy levels that are well separated for the higher energy levels:
ng_vec = np.linspace(0.25, 0.75, 200)
energies = np.array([hamiltonian(Ec, Ej, N, ng).eigenenergies()
for ng in ng_vec])
plot_energies(ng_vec, energies, ymax=(10, 1.1));
Let's tune the system to $n_g = 0.5$ and look at the Hamiltonian and its eigenstates in detail
H = hamiltonian(Ec, Ej, N, 0.5)
H
evals, ekets = H.eigenstates()
The eigenenergies are sorted:
evals
array([ 0.47065435, 1.46676684, 9.01371984, 9.01760693, 25.00520901, 25.00520943, 49.00260427, 49.00260427, 81.00156252, 81.00156252, 121.00104167, 121.00104167, 169.00074405, 169.00074405, 225.00055804, 225.00055804, 289.00043403, 289.00043411, 361.00034728, 361.00347214, 441.00312494])
Only two states have a significant weight in the two lowest eigenstates:
ekets[0].full() > 0.1
array([[False], [False], [False], [False], [False], [False], [False], [False], [False], [False], [ True], [ True], [False], [False], [False], [False], [False], [False], [False], [False], [False]])
abs(ekets[1].full()) > 0.1
array([[False], [False], [False], [False], [False], [False], [False], [False], [False], [False], [ True], [ True], [False], [False], [False], [False], [False], [False], [False], [False], [False]])
We can use these two isolated eigenstates to define a qubit basis:
psi_g = ekets[0] # basis(2, 0)
psi_e = ekets[1] # basis(2, 1)
# psi_g = basis(2, 0)
# psi_e = basis(2, 1)
and corresponding Pauli matrices:
sx = psi_g * psi_e.dag() + psi_e * psi_g.dag()
sz = psi_g * psi_g.dag() - psi_e * psi_e.dag()
and an effective qubit Hamiltonian
evals[1] - evals[0]
0.9961124875822172
H0 = 0.5 * (evals[1] - evals[0]) * sz
A = 0.25 # some driving amplitude
Hd = 0.5 * A * sx # obtained by driving ng(t),
# but now H0 is in the eigenbasis so the drive becomes a sigma_x
Doing this we have a bunch of extra energy levels in the system that aren't involved in the dynamics, but so far they are still in the Hamiltonian.
qubit_evals = H0.eigenenergies()
qubit_evals - qubit_evals[0]
array([0. , 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.49805624, 0.99611249])
fig = plt.figure(figsize=(4, 2))
plot_energy_levels([H0, Hd], fig=fig);
Imagine that we also can drive a $\sigma_x$ type of interaction (e.g., external field):
Heff = [H0, [Hd, "sin(wd*t)"]]
args = {"wd": (evals[1] - evals[0])}
Let's look at the Rabi oscillation dynamics of the qubit when initially placed in the ground state:
psi0 = psi_g
tlist = np.linspace(0.0, 100.0, 500)
result = mesolve(Heff, psi0, tlist, [], [ket2dm(psi_e)], args=args)
/usr/share/miniconda3/envs/test-environment/lib/python3.10/site-packages/qutip/core/coefficient.py:386: UserWarning: Both `cython` and `filelock` are required for compilation of string coefficents. Falling back on `eval`. warnings.warn(
visualize_dynamics(result, r"$\rho_{ee}$");
We can see that only the two selected states are included in the dynamics, and very little leakage to other levels occur.
Instead of keeping all the inactive quantum states in the calculation we can eliminate them, so that we obtain a true two-level system.
np.where(abs(ekets[0].full().flatten()) > 0.1)[0]
array([10, 11])
np.where(abs(ekets[1].full().flatten()) > 0.1)[0]
array([10, 11])
keep_states = np.where(abs(ekets[1].full().flatten()) > 0.1)[0]
H0 = Qobj(H0.full()[keep_states, :][:, keep_states])
H0
Hd = Qobj(Hd.full()[keep_states, :][:, keep_states])
Hd
And if we look at the energy level diagram now we see that we only have two states in the system, as desired.
fig = plt.figure(figsize=(4, 2))
plot_energy_levels([H0, Hd], fig=fig);
Heff = [H0, [Hd, "sin(wd*t)"]]
args = {"wd": (evals[1] - evals[0])}
psi0 = Qobj(psi0.full()[keep_states, :])
psi_e = Qobj(psi_e.full()[keep_states, :])
tlist = np.linspace(0.0, 100.0, 500)
result = mesolve(Heff, psi0, tlist, [], [ket2dm(psi_e)], args=args)
visualize_dynamics(result, r"$\rho_{ee}$");
about()
QuTiP: Quantum Toolbox in Python ================================ Copyright (c) QuTiP team 2011 and later. Current admin team: Alexander Pitchford, Nathan Shammah, Shahnawaz Ahmed, Neill Lambert, Eric Giguère, Boxi Li, Jake Lishman, Simon Cross and Asier Galicia. Board members: Daniel Burgarth, Robert Johansson, Anton F. Kockum, Franco Nori and Will Zeng. Original developers: R. J. Johansson & P. D. Nation. Previous lead developers: Chris Granade & A. Grimsmo. Currently developed through wide collaboration. See https://github.com/qutip for details. QuTiP Version: 5.1.0.dev0+50fc9b7 Numpy Version: 1.22.4 Scipy Version: 1.13.0 Cython Version: 3.0.10 Matplotlib Version: 3.5.2 Python Version: 3.10.4 Number of CPUs: 4 BLAS Info: Generic INTEL MKL Ext: False Platform Info: Linux (x86_64) Installation path: /usr/share/miniconda3/envs/test-environment/lib/python3.10/site-packages/qutip ================================================================================ Please cite QuTiP in your publication. ================================================================================ For your convenience a bibtex reference can be easily generated using `qutip.cite()`