# Wave functions for a flux-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 [22]:
h = 6.626e-34
h_ = h/(2*pi)
e = 1.602e-19
Phi0 = h / (2 * e)
cf = h_

Ic = 1.7e-6
Cj = 700e-15
L = 0.72e-9

mm = Cj * (Phi0/(2*pi))**2 * cf
beta = 2*pi*L*Ic / Phi0
Ej = Phi0/(2*pi) * Ic / cf

phi = 5.0;
gamma = phi / beta;

args = {'Ej': Ej, 'beta': beta, 'gamma': gamma}

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

x_min = -0.7*pi
x_max =  4*pi
N = 750


## Potential¶

In [23]:
def U_flux_biased(x, args):
"""
Flux-biased phase qubit potential
"""
Ej = args['Ej']
beta = args['beta']
gamma = args['gamma']

u = -Ej * (cos(x) - 1 / (2 * beta) * x ** 2 + gamma * x)

return u

In [6]:
x = linspace(x_min, x_max, N)

In [7]:
U = U_flux_biased(x, args)

In [8]:
x_opt_min = optimize.fmin(U_flux_biased, 0.5, (args,))
print("\nFound dmin = %f " % x_opt_min)

Optimization terminated successfully.
Current function value: -9492758355749.710938
Iterations: 29
Function evaluations: 58

Found dmin = 1.361984

In [9]:
x_opt_max = optimize.fmin(lambda x, args: -U_flux_biased(x, args), [2.5], (args,))
print("\nFound dmax = %f" % x_opt_max)

Optimization terminated successfully.
Current function value: 9096370971095.232422
Iterations: 23
Function evaluations: 46

Found dmax = 2.347406

In [10]:
U_min = U_flux_biased(x_opt_min, args)
U_max = U_flux_biased(x_opt_max, args)

dU = U_max - U_min
dx = x_opt_max - x_opt_min

print("Barrier: dU = %f" % (dU / Ej))

Barrier: dU = 0.074707


### Plot potential¶

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

ax.plot(x, U/Ej)
ax.plot(x_opt_min, U_flux_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);


### Plot potential: zoom on local minima¶

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

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

ax.set_ylim(-2, -1.6)
ax.set_xlim(0.5, 3.0)
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_flux_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:  [199, 202, 204, 205, 208, 210, 212]


### Plot 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 eigenstates together with their eigenenergies and the potential¶

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

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

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

ax.plot(x, Y0.real, 'k--')
ax.plot(x, Y.real)

ax.set_ylim(-1.8, -1.7)
ax.set_xlim(0.5,3.0)
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>


### Harmonicity¶

In [18]:
print "Bound energy-levels:"

real((evals[boundidx] - evals[boundidx[0]]) / (evals[boundidx[1]]-evals[boundidx[0]]))

Bound energy-levels:

Out[18]:
array([ 0.        ,  1.        ,  1.97180713,  2.91207597,  3.81601571,
4.6760567 ,  5.4776277 ])
In [19]:
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    1.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    1.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    1.000   -0.000
0.000   -0.000    0.000   -0.000   -0.000   -0.000    1.000

In [20]:
print "expect_pos:\n"
print_matrix(expect_pos)

expect_pos:

1.374    0.109   -0.006    0.000   -0.000    0.000   -0.000
0.109    1.401    0.156   -0.012    0.001   -0.000    0.000
-0.006    0.156    1.430    0.193   -0.018    0.002   -0.000
0.000   -0.012    0.193    1.464    0.225   -0.025    0.003
-0.000    0.001   -0.018    0.225    1.503    0.255   -0.036
0.000   -0.000    0.002   -0.025    0.255    1.552    0.283
-0.000    0.000   -0.000    0.003   -0.036    0.283    1.624

In [21]:
print "expect_kin:\n"
print_matrix(expect_kin)

expect_kin:

0.000    4.531   -0.509    0.054   -0.010    0.003   -0.001
-4.531    0.000    6.280   -0.916    0.116   -0.024    0.007
0.509   -6.280   -0.000    7.514   -1.358    0.201   -0.045
-0.054    0.916   -7.514    0.000    8.434   -1.857    0.324
0.010   -0.116    1.358   -8.434    0.000    9.088   -2.448
-0.003    0.024   -0.201    1.857   -9.088    0.000    9.411
0.001   -0.007    0.045   -0.324    2.448   -9.411    0.000


### Software¶

In [24]:
%reload_ext version_information
%version_information numpy, scipy, matplotlib, wavefunction

Out[24]:
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:08 2013 JST
In [ ]: