Jupyter notebook mode

3.6. Solutions

3.8. Linear variational calculus for the Cooper-pair box.

Click above to activate the Jupyter notebook mode and display the source code.

Look up the problem text in the drop-down frame below.

Problem

The Cooper-pair box consists of a small piece of superconductor (called the island) coupled to a larger piece (called the reservoir) via a Josephson junction: a thin, insulating layer. Superconductivity will be addressed later in this course, and Josephson junctions not at all. Detailed knowledge however is not needed to do this exercise. In a superconductor, Cooper pairs, pairs of electrons with opposite spin and momentum, form a condensate which requires a finite energy to create excitations. The Josephson junction is thin enough for allowing Cooper pairs to tunnel through it, and this tunnelling generates a particular coupling between the two superconducting volumes connected by the junction.

We use the charge basis $\left| n \right\rangle$, where $n$ denotes the number of Cooper pairs that have tunnelled through the junction, and therefore the number of charges that have moved to the small superconductor.

This island has an electrostatic capacity $C$, which means that every electron contributes an amount $E_C = e^2/{2C}$ to the energy. The Hamiltonian reads $$H=4 E_C \sum_{n}(n-n_g)^2 \left| n \right\rangle\left\langle n \right| - \frac{E_J}{2}\sum_{n}\left(\left| n+1 \right\rangle\left\langle n \right|+\left| n \right\rangle\left\langle n+1 \right|\right),$$ with $n={...-1,0,1,...}$ and $n_g$ is some charge offset. The first term reflects the electrostatic charging, and the second reflects tunnelling (also called hopping). The parameter $E_J$ is the so-called Josephson energy. Crucially, note that $n_g$ is a constant (not an operator) under the control of the experimentalist. Note that the charge basis is orthonormal, $\langle n|m\rangle=\delta_{n,m}$.

We will now use linear variational calculus to estimate the ground and first-excited state energies and wave functions. Caution: please do not attempt to solve this problem analytically. Rather, use your favourite mathematical software: Mathematica, Matlab, anything you like! Please print out your code.

1. Write the matrix $H_{qp}$ in the charge basis. The matrix is of course infinite dimensional, show just a subset of it, revealing its basic structure.

2. Write the matrix $S_{qp}$ also in the charge basis.

3. Consider $E_C=10 E_J$. Restrict your trial functions to the subspace $n=-N,...,0,...+N$. Plot the ground and first-excited state energies as a function of $n_g$ in the range $n_g \in \left[-1,1\right]$. Do this for $N=1$, $N=2$ and $N=3$.

4. Repeat for $E_J=10 E_C$.

5. For $E_J=10 E_C$ and $n_g=0.5$, plot the ground and first-excited state energies as a function of $N$ in the range $N \in\left[1,5\right]$. Confirm McDonald's theorem: show that the energies decrease monotonically as you increase $N$.

6. For $E_J=10 E_C$, plot the ground and first-excited state wave functions in the charge basis. That is, plot the coefficient $c_n$ in the expansion $$\left| \Psi \right\rangle=\sum_{n=-N}^{N} c_n \left| n \right\rangle$$ Do this for $n_g=0$ and $n_g=0.5$. What choice of $N$ would you say is accurate enough?

3.8.1.

The matrix $H_{qp}$ in the charge basis has diagonal elements:

$$H_{qq}= 4 E_C\left(q-n_g\right)^2$$

and off-diagonal elements:

$$H_{qp} = -\frac{E_J}{2}\left(\delta_{q, p+1}+\delta_{q, p-1}\right).$$

# Solution 3.8
# Initialization cell (no output) 
import common
common.configure_plotting()
import numpy as np
from numpy.linalg import eig
import matplotlib.pyplot as plt

def build_hamiltonian(N, Ec, Ej, ng):

    ham = np.zeros(shape=(2*N+1, 2*N+1), dtype=complex) 
    ham_kin = np.zeros((2*N+1,), dtype=complex)
    for i in range(0, 2*N+1):
        ham_kin[i] = 4*Ec*(-N+i - ng)**2
    ham_pot = -Ej/2

    np.fill_diagonal(ham, ham_kin)

    offdiag = np.zeros(shape=(2*N,), dtype=complex)
    offdiag[:] = ham_pot
    np.fill_diagonal(ham[1:, :-1], offdiag)
    np.fill_diagonal(ham[:-1, 1:], offdiag)
    return ham

3.8.2.

As mentioned in the problem statement, the charge basis is orthonormal, thus:

$$S_{qp} = \unicode{120793}$$

3.8.3.

# Define constants
Ec = 10
Ej = 1
# Define N
dim = 3
# Define ng and its resolution
steps = 101
ng = np.linspace(-1, 1, steps)

# initialize vectors
ground = np.zeros((dim, steps))
excited = np.zeros((dim, steps))

# build Hamiltonians for N dimensions and solve the eigenvalue problem
for N in range(1, dim+1):
    i = 0;
    for x in ng:
        ham = build_hamiltonian(N,Ec,Ej,x)
        w, v = eig(ham); 
        w = np.msort(w)
        ground[N-1,i] = w.real[0]
        excited[N-1,i] = w.real[1]
        i = i + 1

# plot the ground state and the exctied state for N-dimenstions as function of the charge offset (in the basis of 2e)
fig, ax = plt.subplots()
gr01, = ax.plot(ng, ground[0], linestyle='-', label='Ground state N = 1', color='black')    
gr02, = ax.plot(ng, ground[1], linestyle='--', label='Ground state N = 2' , color='black')
gr02, = ax.plot(ng, ground[2], linestyle=':', label='Ground state N = 3', color='black')
ex01, = ax.plot(ng, excited[0], linestyle='-', label='Excited state N = 1', color='red')    
ex02, = ax.plot(ng, excited[1], linestyle='--', label='Excited state N = 2', color='red')
ex02, = ax.plot(ng, excited[2], linestyle=':', label='Excited state N = 3', color='red')
plt.xlabel(r'Charge offset $n_g$, units of [2e]', fontsize=15)
plt.ylabel(r'Energy [$E_j$]', fontsize=15)
plt.title('CPB ground state and excited state for $E_c = 10E_j$ as function of the charge offset and N', fontsize=12)
plt.legend(loc='best', fontsize=10)
fig.show()

In this case, the difference between the results for $N=1,2,3$ are barely noticeable.

3.8.4.

# Define constants
Ec = 0.1
Ej = 1
# Define N
dim = 3
# Define ng and its resolution
steps = 101
ng = np.linspace(-1, 1, steps)

# initialize vectors
ground = np.zeros((dim, steps))
excited = np.zeros((dim, steps))

# build Hamiltonians for N dimensions and solve the eigenvalue problem
for N in range(1, dim+1):
    i = 0;
    for x in ng:
        ham = build_hamiltonian(N,Ec,Ej,x)
        w, v = eig(ham); 
        w = np.msort(w)
        ground[N-1,i] = w.real[0]
        excited[N-1,i] = w.real[1]
        i = i + 1

# plot the ground state and the exctied state for N-dimenstions as function of the charge offset (in the basis of 2e)
fig, ax = plt.subplots()
gr01, = ax.plot(ng, ground[0], linestyle='-', label='Ground state N = 1', color='black')    
gr02, = ax.plot(ng, ground[1], linestyle='--', label='Ground state N = 2' , color='black')
gr02, = ax.plot(ng, ground[2], linestyle=':', label='Ground state N = 3', color='black')
ex01, = ax.plot(ng, excited[0], linestyle='-', label='Excited state N = 1', color='red')    
ex02, = ax.plot(ng, excited[1], linestyle='--', label='Excited state N = 2', color='red')
ex02, = ax.plot(ng, excited[2], linestyle=':', label='Excited state N = 3', color='red')
plt.xlabel(r'Charge offset $n_g$, units of [2e]', fontsize=15)
plt.ylabel(r'Energy [$E_j$]', fontsize=15)
plt.title('CPB ground state and excited state for $E_c = 0.1 E_j$ as function of the charge offset and N', fontsize=12)
plt.legend(loc='best', fontsize=10)
fig.show()

3.8.5.

# Define constants
Ec = 0.1
Ej = 1

# Define N
dim = 5
ng = 0.5
# initialize vectors
groundN = np.zeros(dim)
excitedN = np.zeros(dim)
Nx = np.zeros(dim)

# build Hamiltonians for N dimensions and solve the eigenvalue problem
for N in range(1, dim+1):
    ham = build_hamiltonian(N,Ec,Ej,ng)
    n, k = eig(ham); 
    n = np.msort(n)
    Nx[N-1] = N
    groundN[N-1] = n.real[0]
    excitedN[N-1] = n.real[1]

# plot the ground state and the excited state as function of N-dimensions
fig, ax = plt.subplots(figsize=(5, 5))
grN, = ax.plot(Nx, groundN, label='Ground state', color='black', marker= 'o')
exN, = ax.plot(Nx, excitedN, label='Excited state', color='red', marker= 'o')
plt.xlabel(r'N', fontsize=15)
plt.ylabel(r'Energy [$E_j$]', fontsize=15)
plt.legend(loc='best', fontsize=10)
plt.title('Energy approximation as f(N) for $E_c = 0.1 E_j$', fontsize=10)
fig.show()

Energy levels converge monotonically towards expected (lower) values as a function of N, confirming the McDonald's theorem.

3.8.6.

# Define constants
Ec = 0.1
Ej = 1
ng = 0
# Set number of dimensions N
dim = 8
# initialize the N-dimension axis
Nx = np.zeros(2*dim+1)
for i in range(0, 2*dim+1):
    Nx[i] = -dim + i

# build Hamiltonians for N dimensions and solve the eigenvalue problem
ham = build_hamiltonian(dim,Ec,Ej,ng)
evals, evecs = eig(ham)
sorted_indeces = np.argsort(evals.real)
evals_sorted = evals[sorted_indeces]
evecs_sorted = evecs[:,sorted_indeces]    

# plot the coefficients of the ground state and the excited state in the charge basis 
fig, ax = plt.subplots()
grN, = ax.plot(Nx, evecs_sorted.real[:,0], label='Ground state', color='black',  marker= 'o')
exN, = ax.plot(Nx, evecs_sorted.real[:,1], label='Excited state', color='red',  marker= 'o')
plt.title('Expansion of the ground and the first-excited states in the charge basis, for $n_g = 0$', fontsize=12)
plt.xlabel(r'N', fontsize=15)
plt.ylabel(r'$c_n$', fontsize=15)
plt.grid(axis='x', color='0.95')
plt.legend(loc='best', fontsize=10)
fig.show()

# Define constants
Ec = 0.1
Ej = 1
ng = 0.5
# Set number of dimensions N
dim = 8
# initialize the N-dimension axis
Nx = np.zeros(2*dim+1)
for i in range(0, 2*dim+1):
    Nx[i] = -dim + i

# build Hamiltonians for N dimensions and solve the eigenvalue problem
ham = build_hamiltonian(dim,Ec,Ej,ng)
evals, evecs = eig(ham)
sorted_indeces = np.argsort(evals.real)
evals_sorted = evals[sorted_indeces]
evecs_sorted = evecs[:,sorted_indeces]    

# plot the coefficients of the ground state and the excited state in the charge basis
fig, ax = plt.subplots()
grN, = ax.plot(Nx, evecs_sorted.real[:,0], label='Ground state', color='black',  marker= 'o')
exN, = ax.plot(Nx, evecs_sorted.real[:,1], label='Excited state', color='red',  marker= 'o')
plt.title('Expansion of the ground and the first-excited states in the charge basis, for $n_g = 0.5$', fontsize=12)
plt.xlabel(r'N', fontsize=15)
plt.ylabel(r'$c_n$', fontsize=15)
plt.grid(axis='x', color='0.95')
plt.legend(loc='best', fontsize=10)
fig.show()

Note that the expansion coefficients for the ground (the first-excited) state are symmetric (antisymmetric) around $N = n_g$.