Wavefunctions for a current-biased Josephson-junction phase qubit

Robert Johansson ([email protected])

In [1]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
In [2]:
from scipy import *
from scipy import optimize
In [3]:
from wavefunction import *
from wavefunction.wavefunction1d import *

Problem Parameters

In [4]:
h = 6.626e-34
h_ = h/(2*pi)
e = 1.602e-19
Phi0 = h / (2 * e)
cf = 1

ceta = 0.05
J = 0.9725
Ic = 13.3e-6
Ib = J * Ic
Cj = 4.3e-12

mm = Cj * (1 + ceta) * (Phi0/(2*pi))**2 * cf
Ej = Phi0/(2*pi) * Ic / cf

args = {'Ej': Ej, 'Ic': Ic, 'Ib': Ib}

k = -h_ ** 2 / (2 * mm)

x_min = -pi
x_max =  pi
N = 750 

Potential

In [5]:
def U_current_biased(x, args):
    """
    Potential for a current-biased phase qubit potential
    (the washboard potential)
    """
    Ej = args['Ej']
    Ic = args['Ic']
    Ib = args['Ib']

    u = - Ej * (cos(x) + Ib / Ic * x)

    return u
In [6]:
x = linspace(x_min, x_max, N)
In [7]:
U = U_current_biased(x, args)
In [11]:
# analytical formula for minima points of potential
x_opt_min = arcsin(Ib/Ic)
print("x_opt_min = %f " % x_opt_min)

x_opt_max = -math.asin(Ib/Ic) + pi
print("x_opt_max = %f" % x_opt_max)

U_min = U_current_biased(x_opt_min, args)
U_max = U_current_biased(x_opt_max, args)
dU = U_max - U_min
dx = x_opt_max - x_opt_min

print("dU = %f" % dU)
x_opt_min = 1.335735 
x_opt_max = 1.805858
dU = 0.000000

Plot potential

In [12]:
fig, ax = subplots()

ax.plot(x, U/Ej)
ax.plot(x_opt_min, U_current_biased(x_opt_min, args) / Ej, 'o')

ax.set_ylim(-1.55, -1.5)
ax.set_xlim(0.8,2.2)
ax.set_xlabel(r'$x$', fontsize=18)
ax.set_ylabel(r'$U(x)/E_J$', fontsize=18);

Calculate the eigenfunctions

In [13]:
u = assemble_u_potential(N, U_current_biased, x, args)

K = assemble_K(N, k, x_min, x_max)

V = assemble_V(N, u, x_min, x_max)

H = K + V

evals, evecs = solve_eigenproblem(H)

Identify bound states

In [14]:
nbound = 0
boundidx = []

for i in range(0,N):
    if evals[i] > U_min - 0.5 * dU and evals[i] < U_max + 0.5 * dU:
        if inner(evecs[i] * (x_opt_min-dx < x) * (x < x_opt_max), evecs[i]) > 0.95:
            nbound = nbound + 1
            boundidx.append(i)

print "Found bound states: ", boundidx
Found bound states:  [91, 92, 94, 96, 97, 99, 101, 102, 104]

Plot bounded eigenstates

In [15]:
NN = len(boundidx)

fig, axes = subplots(NN, 1, figsize=(10, NN * 1), sharex=True, sharey=True)

for n, m in enumerate(boundidx):
    axes[n].plot(x, real(evecs[m]))

axes[n].set_xlim(0.5, 2.5)
axes[n].set_xlabel(r'$x$', fontsize=18)

fig.tight_layout();

Plot bounded eigenstates together with the potential

In [16]:
fig, ax = subplots(figsize=(12,8))

ax.plot(x, U/Ej, 'k')
ax.plot(x_opt_min, U_current_biased(x_opt_min, args) / Ej, '.')

for n, m in enumerate(boundidx):
    Y0 = evals[m]/Ej * ones(shape(x))
    Y = 0.001 * evecs[m] + Y0

    ax.plot(x, Y0.real, 'k--')
    ax.plot(x, Y.real)
    
ax.set_ylim(-1.535, -1.52)
ax.set_xlim(0.8,2.2)
ax.set_xlabel(r'$x$', fontsize=18)
ax.set_ylabel(r'$U(x)$', fontsize=18);

Evaluate transition matrix elements

In [17]:
inner_prod = zeros((len(boundidx), len(boundidx))).astype(np.complex)
expect_pos = zeros((len(boundidx), len(boundidx))).astype(np.complex)
expect_kin = zeros((len(boundidx), len(boundidx))).astype(np.complex)

for i, l in enumerate(boundidx):
    for j, k in enumerate(boundidx):
    
        psi_l = wavefunction_normalize(x, evecs[l])
        psi_k = wavefunction_normalize(x, evecs[k])
    
        inner_prod[i,j] = inner_product(x, psi_l, psi_k)                 # <psi_l|psi_k>
        expect_pos[i,j] = inner_product(x, psi_l, x * psi_k)             # <psi_l|x|psi_k>
        expect_kin[i,j] = inner_product(x, psi_l, derivative(x, psi_k))  # <psi_l|p|psi_k>
In [18]:
print("Bound energy-levels:")

(evals[boundidx] - evals[boundidx[0]]).real / Ej
Bound energy-levels:
Out[18]:
array([ 0.        ,  0.00107872,  0.00213317,  0.003161  ,  0.00415908,
        0.005123  ,  0.00604603,  0.00691584,  0.00770281])
In [19]:
print("Harmonicity:")

((evals[boundidx] - evals[boundidx[0]]) / (evals[boundidx[1]]-evals[boundidx[0]])).real
Harmonicity:
Out[19]:
array([ 0.        ,  1.        ,  1.97749688,  2.9303132 ,  3.85555056,
        4.74913082,  5.60479901,  6.41112924,  7.14067069])
In [20]:
print "inner_prod:\n"
print_matrix(inner_prod)
inner_prod:

 1.000   -0.000   -0.000   -0.000   -0.000    0.000    0.000    0.000   -0.000  
-0.000    1.000    0.000   -0.000   -0.000   -0.000   -0.000    0.000    0.000  
-0.000    0.000    1.000    0.000   -0.000   -0.000    0.000    0.000    0.000  
-0.000   -0.000    0.000    1.000    0.000    0.000   -0.000   -0.000    0.000  
-0.000   -0.000   -0.000    0.000    1.000   -0.000    0.000    0.000   -0.000  
 0.000   -0.000   -0.000    0.000   -0.000    1.000   -0.000    0.000    0.000  
 0.000   -0.000    0.000   -0.000    0.000   -0.000    1.000   -0.000   -0.000  
 0.000    0.000    0.000   -0.000    0.000    0.000   -0.000    1.000   -0.000  
-0.000    0.000    0.000    0.000   -0.000    0.000   -0.000   -0.000    1.000  
In [21]:
print "expect_pos:\n"
print_matrix(expect_pos)
expect_pos:

 1.341    0.049   -0.002    0.000   -0.000    0.000   -0.000    0.000   -0.000  
 0.049    1.351    0.070   -0.005    0.000   -0.000    0.000   -0.000    0.000  
-0.002    0.070    1.363    0.086   -0.007    0.001   -0.000    0.000   -0.000  
 0.000   -0.005    0.086    1.376    0.100   -0.010    0.001   -0.000    0.000  
-0.000    0.000   -0.007    0.100    1.390    0.113   -0.013    0.002   -0.000  
 0.000   -0.000    0.001   -0.010    0.113    1.407    0.125   -0.017    0.002  
-0.000    0.000   -0.000    0.001   -0.013    0.125    1.428    0.137   -0.023  
 0.000   -0.000    0.000   -0.000    0.002   -0.017    0.137    1.456    0.146  
-0.000    0.000   -0.000    0.000   -0.000    0.002   -0.023    0.146    1.515  
In [22]:
print "expect_kin:\n"
print_matrix(expect_kin)
expect_kin:

-0.000    10.132   -1.024    0.103   -0.019    0.004   -0.001    0.000   -0.000  
-10.132   -0.000    14.099   -1.830    0.215   -0.044    0.012   -0.003    0.001  
 1.024   -14.099    0.000    16.956   -2.683    0.361   -0.081    0.024   -0.007  
-0.103    1.830   -16.956    0.000    19.171   -3.617    0.551   -0.133    0.044  
 0.019   -0.215    2.683   -19.171    0.000    20.899   -4.670    0.810   -0.207  
-0.004    0.044   -0.361    3.617   -20.899    0.000    22.167   -5.909    1.194  
 0.001   -0.012    0.081   -0.551    4.670   -22.167   -0.000    22.850   -7.450  
-0.000    0.003   -0.024    0.133   -0.810    5.909   -22.850   -0.000    22.068  
 0.000   -0.001    0.007   -0.044    0.207   -1.194    7.450   -22.068    0.000  

Software

In [23]:
%reload_ext version_information
%version_information numpy, scipy, matplotlib, wavefunction
Out[23]:
SoftwareVersion
Python2.7.4 (default, Apr 19 2013, 18:28:01) [GCC 4.7.3]
IPython2.0.0-dev
OSposix [linux2]
numpy1.8.0.dev-895866d
scipy0.13.0.dev-7c6d92e
matplotlib1.3.x
wavefunction1.0.0
Mon Sep 02 15:16:18 2013 JST
In [ ]: