Calculations of the Flux qubit wave function¶

In [3]:
%matplotlib inline

In [4]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm

In [5]:
from wavefunction import *
from wavefunction.wavefunction2d import *


Introduction¶

Here follow Orlando et al., PRB 60, 15398 (1999) and numerically calculate the wavefunctions and energy-levels of the Flux qubit device using a two-dimensional potential of the superconducting Flux qubit circuit. The two generalized coordinates are the gauge invariant phases across the Josephson junctions in the device.

Original coordinates¶

The Hamiltonian in the original phase coordinates is

$H_t = \frac{1}{2} P^T M^{-1} P + E_J(2 + \alpha - \cos\varphi_1 - \cos\varphi_2 - \alpha\cos(2\pi f + \varphi_1 - \varphi_2))$

where $P = (-i\hbar\partial_{\varphi_1}, -i\hbar\partial_{\varphi_2})^T$ is the generalized momentum vector and the mass matrix is

$M = (\Phi_0/2\pi)^2 C \begin{pmatrix} 1 + \alpha + \gamma & -\alpha \\ -\alpha & 1 + \alpha + \gamma\end{pmatrix}$

and

$M^{-1} = \frac{1}{(\Phi_0/2\pi)^2 C} \frac{1}{(1 + \gamma)(1 + \gamma + 2\alpha)} \begin{pmatrix} 1 + \alpha + \gamma & \alpha \\ \alpha & 1 + \alpha + \gamma\end{pmatrix}$

This gives

$H_t = - \frac{4E_C }{(1 + \gamma)(1 + \gamma + 2\alpha)} \left[ (1 + \gamma + \alpha)\left(\frac{\partial}{\partial\varphi_1}\right)^2 + 2\alpha\frac{\partial}{\partial\varphi_1}\frac{\partial}{\partial\varphi_2} + (1 + \gamma + \alpha)\left(\frac{\partial}{\partial\varphi_2}\right)^2 \right] • E_J(2 + \alpha - \cos\varphi_1 - \cos\varphi_2 - \alpha\cos(2\pi f + \varphi_1 - \varphi_2))$

For the numerical calculations we use units of $E_J$

$H_t = - \frac{4}{(1 + \gamma)(1 + \gamma + 2\alpha)}\frac{E_C}{E_J} \left[ (1 + \gamma + \alpha)\left(\frac{\partial}{\partial\varphi_1}\right)^2 + 2\alpha\frac{\partial}{\partial\varphi_1}\frac{\partial}{\partial\varphi_2} + (1 + \gamma + \alpha)\left(\frac{\partial}{\partial\varphi_2}\right)^2 \right] • (2 + \alpha - \cos\varphi_1 - \cos\varphi_2 - \alpha\cos(2\pi f + \varphi_1 - \varphi_2)))$
In [6]:
args = {'Ec': 1/80.0, 'Ej': 1.0, 'alpha': 0.8, 'gamma': 0.0, 'f': 0.50}


Plot the potential¶

In [7]:
x1 = np.linspace(-2*np.pi, 2*np.pi, 100)
x2 = np.linspace(-2*np.pi, 2*np.pi, 100)
X1, X2 = np.meshgrid(x1, x2)

In [8]:
def U_flux_qubit(x1, x2, args):
globals().update(args)
return 2 + alpha - np.cos(x1) - np.cos(x2) - alpha * np.cos(2 * np.pi * f + x1 - x2)

In [9]:
Z = np.real(U_flux_qubit(X1, X2, args))
fig, ax = plt.subplots()
p = ax.pcolor(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
ax.set_xlabel(r'$x_1$', fontsize=16)
ax.set_ylabel(r'$x_2$', fontsize=16)
cb = fig.colorbar(p, ax=ax)


Assembling Kinetic matrix¶

$H_t = - \frac{4}{(1 + \gamma)(1 + \gamma + 2\alpha)}\frac{E_C}{E_J} \left[ (1 + \gamma + \alpha)\left(\frac{\partial}{\partial\varphi_1}\right)^2 + 2\alpha\frac{\partial}{\partial\varphi_1}\frac{\partial}{\partial\varphi_2} + (1 + \gamma + \alpha)\left(\frac{\partial}{\partial\varphi_2}\right)^2 \right] • (2 + \alpha - \cos\varphi_1 - \cos\varphi_2 - \alpha\cos(2\pi f + \varphi_1 - \varphi_2)))$

$\displaystyle K_{n_1, n_2}^{m_1, m_2} = \delta_{n_1,m_1} \delta_{n_2,m_2} \left(k_{11}\left(\frac{2\pi m_1}{T_{x_1}}\right)^2 + k_{12}\frac{(2\pi)^2 m_1m_2}{T_{x_1}T_{x_2}} + k_{22}\left(\frac{2\pi m_2}{T_{x_2}}\right)^2\right)$

In [10]:
# pick a truncation of the fourier series: n1 = [-L1, ..., L1], n2 = [-L2, ..., L2]
L1 = 5
L2 = 5

# pick periods for the coordinates
Tx1 = Tx2 = 2 * np.pi

#
k11 = k22 = 4 * (Ec / Ej) * (1 + alpha + gamma) / ((1 + gamma) * (1 + 2 * alpha + gamma))
k12 = 4 * (Ec / Ej) * 2 * alpha / ((1 + gamma) * (1 + 2 * alpha + gamma))

In [11]:
K = assemble_K(L1, L2, k11, k12, k22, Tx1, Tx2)


Assembling the potential contribution¶

The flux qubit potential we consider here is

$\displaystyle U(\varphi_1, \varphi_2)/E_J = 2 + \alpha - \cos(\varphi_1) - \cos(\varphi_2) - \alpha \cos(2\pi f + \varphi_1 - \varphi_2)$

To obtain the Fourier series expansion of $U(\varphi_1, \varphi_2)$, we first write the $\cos$ expressions as exponential functions

$\displaystyle U(\varphi_1, \varphi_2)/E_J = (2 + \alpha) • \frac{1}{2}(e^{i\varphi_1} + e^{-i\varphi_1}) • \frac{1}{2}(e^{i\varphi_2} + e^{-\varphi_2}) • \alpha\frac{1}{2}(e^{i2\pi f}e^{i\varphi_1}e^{-i\varphi_2} + e^{-i2\pi f}e^{-i\varphi_1}e^{i\varphi_2})$

$\displaystyle U(\varphi_1, \varphi_2)/E_J = \sum_{n_1}\sum_{n_2} u_{n_1n_2} e^{in_1\varphi_1}e^{in_2\varphi_2}$

we can identity

$\displaystyle u{n_1n_2} = (2 + \alpha) \delta{0,n1}\delta{0,n_2} • \frac{1}{2}(\delta{1, n_1} + \delta{-1, n1})\delta{0,n_2} • \frac{1}{2}(\delta{1, n_2} + \delta{-1, n2})\delta{0,n_1} • \frac{\alpha}{2} (e^{i2\pi f} \delta{1,n_1}\delta{-1,n2} + e^{-i2\pi f}\delta{-1,n1}\delta{+1,n_2})$
In [12]:
def assemble_u_flux_qubit(L1, L2, args, sparse=False):

globals().update(args)

L1n = 2 * L1 + 1
L2n = 2 * L2 + 1

u = np.zeros((L1n, L2n), dtype=np.complex)

for n1 in np.arange(-L1, L1+1):
N1 = n1 + L1
for n2 in np.arange(-L2, L2+1):
N2 = n2 + L2
u[N1, N2] = (2 + alpha) * delta(0, n1) * delta(0, n2) + \
- 0.5 * (delta(1, n1) + delta(-1, n1)) * delta(0, n2) + \
- 0.5 * (delta(1, n2) + delta(-1, n2)) * delta(0, n1)+ \
- 0.5 * alpha * np.exp(+2j * np.pi * f) * delta(+1, n1) * delta(-1, n2) + \
- 0.5 * alpha * np.exp(-2j * np.pi * f) * delta(-1, n1) * delta(+1, n2)
return u

In [13]:
u = assemble_u_flux_qubit(L1, L2, args)


Evaluate and plot the Fourier series representation of the flux qubit potential¶

$\displaystyle U(x_p, x_m) = \sum_{n_1}\sum_{n_2} u_{n_1n_2} e^{in_1x_p}e^{in_2x_m}$

In [14]:
Z = np.real(evalute_fourier_series(X1, X2, L1, L2, u))
fig, ax = plt.subplots()
p = ax.pcolor(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
cb = fig.colorbar(p, ax=ax)


Potential contribution to the eigenvalue problem¶

Given this Fourier series we can calcuate the potential contribution to the eigenvalue problem as

$\displaystyle V_{n_1, n_2}^{m_1, m_2} = u_{n_1-m_1,n_2-m_2}$

In [15]:
V = assemble_V(L1, L2, u)


and specifically for the flux qubit potential

$\displaystyle V{n_1, n_2}^{m_1, m_2} = (2 + \alpha) \delta{0,n1-m_1}\delta{0,n_2-m_2} • \frac{1}{2}(\delta{1, n_1-m_1} + \delta{-1, n1-m_1})\delta{0,n_2-m_2} • \frac{1}{2}(\delta{1, n_2-m_2} + \delta{-1, n2-m_2})\delta{0,n_1-m_1} • \frac{\alpha}{2} (e^{i2\pi f} \delta{1,n_1-m_1}\delta{-1,n2-m_2} + e^{-i2\pi f}\delta{-1,n1-m_1}\delta{+1,n_2-m_2})$
In [16]:
def assemble_V_flux_qubit(L1, L2, args, sparse=False):

globals().update(args)

L1n = 2 * L1 + 1
L2n = 2 * L2 + 1

V = np.zeros((L1n*L1n, L2n*L2n), dtype=np.complex)

for n1 in np.arange(-L1, L1+1):
for n2 in np.arange(-L1, L1+1):
N = index_m2v(L1, n1, n2)
for m1 in np.arange(-L2, L2+1):
for m2 in np.arange(-L2, L2+1):
M = index_m2v(L2, m1, m2)

V[N,M] = (2 + alpha) * delta(m1, n1) * delta(m2, n2) + \
- 0.5 * (delta(m1 + 1, n1) + delta(m1 - 1, n1)) * delta(m2, n2) + \
- 0.5 * (delta(m2 + 1, n2) + delta(m2 - 1, n2)) * delta(m1, n1) + \
- 0.5 * alpha * np.exp(+2j * np.pi * f) * delta(m1 + 1, n1) * delta(m2 - 1, n2) + \
- 0.5 * alpha * np.exp(-2j * np.pi * f) * delta(m1 - 1, n1) * delta(m2 + 1, n2)

return V

In [17]:
V_fq = assemble_V_flux_qubit(L1, L2, args)

In [18]:
abs(V - V_fq).max()

Out[18]:
0.0

Solving the eigenvalue problem¶

We now want to solve the eigenvalue problem

$H \Psi = E \Psi$

where $H$ is the matrix representation of the Hamiltonian assembled from $K$ and $V$ above.

In [19]:
H = K + V

In [20]:
vals, vecs = solve_eigenproblem(H)


Plot energy levels and potential¶

In [21]:
fig, ax = plt.subplots()

U_x = np.real(U_flux_qubit(x1, -x1, args))

for val in vals:
ax.plot(x1, np.real(val) * np.ones_like(x1), 'k')

ax.plot(x1, U_x, label="potential", lw='2')

ax.axis('tight')
ax.set_ylim(min(vals[0], U_x.min()), U_x.max())
ax.set_ylabel(r'$U(x_1, x_2=x_1)$', fontsize=16)
ax.set_xlabel(r'$x_1$', fontsize=16);

In [22]:
Nstates = 4

fig, axes = plt.subplots(Nstates, 3, figsize=(16, 5 * Nstates))

for n in range(Nstates):

psi = convert_v2m(L1, L2, vecs[n])

Z = np.real(evalute_fourier_series(X1, X2, L1, L2, u))
PSI = evalute_fourier_series(X1, X2, L1, L2, psi)

p = axes[n, 0].pcolor(X1/(2*np.pi), X2/(2*np.pi), np.real(PSI),
cmap=cm.RdBu, vmin=-abs(PSI.real).max(), vmax=abs(PSI.real).max())
c = axes[n, 0].contour(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
cb = fig.colorbar(p, ax=axes[n, 0])

p = axes[n, 1].pcolor(X1/(2*np.pi), X2/(2*np.pi), np.imag(PSI),
cmap=cm.RdBu, vmin=-abs(PSI.imag).max(), vmax=abs(PSI.imag).max())
c = axes[n, 1].contour(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
cb = fig.colorbar(p, ax=axes[n, 1])

p = axes[n, 2].pcolor(X1/(2*np.pi), X2/(2*np.pi), np.abs(PSI), cmap=cm.Blues, vmin=0, vmax=abs(PSI).max())
c = axes[n, 2].contour(X1/(2*np.pi), X2/(2*np.pi), Z, cmap=cm.RdBu, vmin=Z.min(), vmax=Z.max())
cb = fig.colorbar(p, ax=axes[n, 2])

axes[n, 0].set_ylabel("%d state" % n);

axes[0, 0].set_title(r"$\mathrm{re}\;\Psi(\varphi_p, \varphi_m)$", fontsize=16);
axes[0, 1].set_title(r"$\mathrm{im}\;\Psi(\varphi_p, \varphi_m)$", fontsize=16);
axes[0, 2].set_title(r"$|\Psi(\varphi_p, \varphi_m)|^2$", fontsize=16);