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,n_1}\delta_{0,n_2} - \frac{1}{2}(\delta_{1, n_1} + \delta_{-1, n_1})\delta_{0,n_2} - \frac{1}{2}(\delta_{1, n_2} + \delta_{-1, n_2})\delta_{0,n_1} - \frac{\alpha}{2} (e^{i2\pi f} \delta_{1,n_1}\delta_{-1,n_2} + e^{-i2\pi f}\delta_{-1,n_1}\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,n_1-m_1}\delta_{0,n_2-m_2} - \frac{1}{2}(\delta_{1, n_1-m_1} + \delta_{-1, n_1-m_1})\delta_{0,n_2-m_2} - \frac{1}{2}(\delta_{1, n_2-m_2} + \delta_{-1, n_2-m_2})\delta_{0,n_1-m_1} - \frac{\alpha}{2} (e^{i2\pi f} \delta_{1,n_1-m_1}\delta_{-1,n_2-m_2} + e^{-i2\pi f}\delta_{-1,n_1-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);