Jupyter notebook mode
3.6. Solutions¶
3.8. Linear variational calculus for the Cooperpair box.¶
Click above to activate the Jupyter notebook mode and display the source code.
Look up the problem text in the dropdown frame below.
Problem
The Cooperpair 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 , where 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 , which means that every electron contributes an amount to the energy. The Hamiltonian reads with and is some charge offset. The first term reflects the electrostatic charging, and the second reflects tunnelling (also called hopping). The parameter is the socalled Josephson energy. Crucially, note that is a constant (not an operator) under the control of the experimentalist. Note that the charge basis is orthonormal, .
We will now use linear variational calculus to estimate the ground and firstexcited 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.

Write the matrix in the charge basis. The matrix is of course infinite dimensional, show just a subset of it, revealing its basic structure.

Write the matrix also in the charge basis.

Consider . Restrict your trial functions to the subspace . Plot the ground and firstexcited state energies as a function of in the range . Do this for , and .

Repeat for .

For and , plot the ground and firstexcited state energies as a function of in the range . Confirm McDonald's theorem: show that the energies decrease monotonically as you increase .

For , plot the ground and firstexcited state wave functions in the charge basis. That is, plot the coefficient in the expansion Do this for and . What choice of would you say is accurate enough?
3.8.1.¶
The matrix in the charge basis has diagonal elements:
and offdiagonal elements:
# 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:
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[N1,i] = w.real[0] excited[N1,i] = w.real[1] i = i + 1 # plot the ground state and the exctied state for Ndimenstions 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 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[N1,i] = w.real[0] excited[N1,i] = w.real[1] i = i + 1 # plot the ground state and the exctied state for Ndimenstions 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[N1] = N groundN[N1] = n.real[0] excitedN[N1] = n.real[1] # plot the ground state and the excited state as function of Ndimensions 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 Ndimension 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 firstexcited 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 Ndimension 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 firstexcited 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 firstexcited) state are symmetric (antisymmetric) around .