I have recently been delving into tight-binding modeling of Josephson junctions. As I have been learning to solve this system, and being rusty on my tight-binding basics (being an experimentalist), I have found it remarkably difficult to find examples of what the Hamiltonian matrices for these systems should look like, either in publications or through online lectures/notes. While there are several papers which use this method for solving all sorts of Josephson junctions (SIS, SNS, SFS, SFIS) and write the in-line Hamiltonian quite clearly, the whole host of symbols in these equations (\( \dagger , *, i, j, \langle\rangle, \sigma \)) can leave a newbie doubting their own construction of the matrices. So, here I am writing things as explicitly as possible to aid someone else who may not be sure about what they are doing.
The BdG Hamiltonian, \( \hat{H}_{BdG} \), is used to model superconductivity and electron dynamics. The equation is written as:
\[
\hat{H}_{BdG} = \sum_{i,\sigma} \epsilon_i \hat{c}_{i\sigma}^\dagger \hat{c}_{i\sigma}
- \sum_{\langle i,j\rangle,\sigma} \left( t_{i,j} \hat{c}_{i,\sigma}^\dagger \hat{c}_{j,\sigma} + \text{h.c.} \right)
+ \sum_{i} \left( \Delta_i \hat{c}_{i,\downarrow}^\dagger \hat{c}_{i,\uparrow}^\dagger +\Delta_i^* \hat{c}_{i,\uparrow} \hat{c}_{i,\downarrow}\right).
\]
This Hamiltonian has three key components:
- Onsite potentials (\( \epsilon_i \)): Represent energies like the chemical potential or external fields.
- Hopping term (\( t_{i,j} \)): Captures the kinetic energy contribution from electron hopping between neighboring sites.
- Pair potential (\( \Delta_i \)): Describes superconductivity within the mean-field approximation.
In most cases, the system is assumed to be homogeneous, with uniform hopping (\( t_{i,j} = t \)) and a constant pair potential (\( \Delta_i = \Delta_0 \)). This simplification helps in understanding superconducting systems while retaining the essential physics. Note that in the following discussion, I will drop the "hats" from the operator.
BdG Hamiltonian – Single Atom
The BdG hamiltonian is written in the basis of the Nambu spinor, wherein, each lattice site is represented through four elements - creation and annihilation of electrons, of spins up and down. The Nambu spinor is given by -
$$ \Psi = \begin{pmatrix} c_\uparrow \\ c_\downarrow \\ c^\dagger_\uparrow \\ c^\dagger_\downarrow \end{pmatrix} $$
where \( c_\uparrow \) and \( c_\downarrow \) are annihilation operators for spin-up and spin-down electrons, and \( c^\dagger_\uparrow \) and \( c^\dagger_\downarrow \) are their corresponding creation operators. In this basis, the full BdG hamiltonian can be written as
$$ H = \sum_{i,j} c^\dagger_i h_{ij} c_j $$
where \( h_{ij} \) are elements of the Hamiltonian matrix \( H_{BdG} \). For a single-site system, \( H_{BdG} \) is a 4 \( \times \) 4 matrix which can be written in a block form as
$$ H_{BdG} = \begin{pmatrix} H_0 & \Delta \\ \Delta^\dagger & -H_0^T \end{pmatrix} $$
Although the block form is not a major simplification for a single-site system, its form remains the same for a long chain of atoms too, and hence its good to get introduced to it already. Anyway, for the single-site system, the blocks in the Hamiltonian above are 2 \( \times \) 2 each, and given by
$$ H_0 = \begin{pmatrix} \epsilon_k & 0 \\ 0 & \epsilon_k \end{pmatrix}, \quad \Delta = \begin{pmatrix} 0 & \Delta_0 \\ -\Delta_0 & 0 \end{pmatrix}$$
Here, \( \epsilon_k \) is the onsite energy term and \( \Delta_0 \) is the superconducting gap magnitude. The anti-symmetry of the off-diagonal terms in \( \Delta \) ensures the spin-singlet pairing, i.e. \( \Delta_{\uparrow\downarrow} = -\Delta_{\downarrow\uparrow} \). To get a more intuitive picture of how the coefficients correspond to the Nambu spinor elements, the array representation below could be helpful
$$
\begin{array}
cc_{\uparrow}^\dagger (\epsilon_k) c_{\uparrow} & c_{\uparrow}^\dagger(0)c_{\downarrow} & c_{\uparrow}^\dagger(0)c_{\uparrow}^\dagger & c_{\uparrow}^\dagger (\Delta_0) c_{\downarrow}^\dagger \\
c_{\downarrow}^\dagger(0)c_{\uparrow} & c_{\downarrow}^\dagger (\epsilon_k) c_{\downarrow} & c_{\downarrow}^\dagger (-\Delta_0) c_{\uparrow}^\dagger & c_{\downarrow}^\dagger(0)c_{\downarrow}^\dagger \\
c_{\uparrow}(0)c_{\uparrow} & c_{\uparrow}(-\Delta_0)c_{\downarrow} & c_{\uparrow}(-\epsilon_k)c_{\uparrow}^\dagger & c_{\uparrow}(0)c_{\downarrow}^\dagger \\
c_{\downarrow}(\Delta_0)c_{\uparrow} & c_{\downarrow}(0)c_{\downarrow} & c_{\downarrow}(0)c_{\uparrow}^\dagger & c_{\downarrow}(-\epsilon_k)c_{\downarrow}^\dagger
\end{array}
$$
Two sites - Hopping term
When we go from a single-site to two sites, we need to additionally consider the energy (\( -t \)) associated with electrons hopping between the two sites. Here, we assume that electron hopping preserves the spin, i.e. the hopping energy is only associated with \( c_{i,\sigma}^\dagger c_{j,\sigma} \) for electrons and \( c_{i,\sigma} c_{j,\sigma}^\dagger \) for holes, but not with \( c_{i,\sigma}^\dagger c_{j,-\sigma} \) or \( c_{i,\sigma} c_{j,-\sigma}^\dagger \), where i, j are the site indices and \( \sigma \in \{\uparrow, \downarrow\} \). The Hamiltonian for the two-site system when represented as in (3) is now an 8 \( \times \) 8 matrix, with \( H_0 \) and \( \Delta \) taking the form
$$
H_0 = \begin{pmatrix}
\epsilon_k & 0 & -t & 0\\
0 & \epsilon_k & 0 & -t\\
-t & 0 & \epsilon_k & 0\\
0 & -t & 0 & \epsilon_k\\
\end{pmatrix}, \quad
\Delta = \begin{pmatrix}
0 & \Delta_0 & 0 & 0\\
-\Delta_0 & 0 & 0 & 0\\
0 & 0 & 0 & \Delta_0\\
0 & 0 & -\Delta_0 & 0\\
\end{pmatrix}.
$$
With the self-energy, superconducting gap and hopping terms accounted for, the tight-binding Hamiltonian for longer chains can be written by simply expanding the matrix and filling the three kinds of terms correctly. In most cases, the hopping term is limited to the nearest neighbor sites (usually denoted as \( \langle i,j\rangle \)). These three types of terms also happen to be all we need to model an S-N-S Josephson junction, as we will see next.
Modelling an SNS Josephson junction
We can model an SNS system by appropriately choosing the values of \( \epsilon_k, t \) and \( \Delta_0 \) for the superconducting and normal metal sites. In particular, \( \Delta_0 \) is set to zero for the normal metal region. For an insulating barrier (SIS junction), we should additionally choose \( \epsilon_{k,I} >> \epsilon_{k,S} \) and \( t_{I} << t_{S} \). Furthermore, the two superconducting regions can have a phase difference, \( \phi \). This can be accounted for by setting \( \Delta_{S1} = \Delta_0e^{-i\phi/2} \) and \( \Delta_{S2} = \Delta_0e^{i\phi/2} \).
As a toy example, I will try to write down the full Hamiltonian matrix for an SNS junction with a single-site in each region, which already ends up being 12 \( \times \) 12 in size. The basis for this Hamiltonian shall be
$$ \Psi = \begin{pmatrix}
c_{S1\uparrow} & c_{S1\downarrow} & c_{N\uparrow} & c_{N\downarrow} & c_{S2\uparrow} & c_{S2\downarrow} & c_{S1\uparrow}^\dagger & c_{S1\downarrow}^\dagger & c_{N\uparrow}^\dagger & c_{N\downarrow}^\dagger & c_{S2\uparrow}^\dagger & c_{S2\downarrow}^\dagger
\end{pmatrix}^T $$
Once again, the full Hamiltonian can be expressed in block form as shown in (3), with \( H_0 \) and \( \Delta \) taking the form
$$
H_0 = \begin{pmatrix}
\epsilon_S & 0 & -t & 0 & 0 & 0\\
0 & \epsilon_S & 0 & -t & 0 & 0\\
-t & 0 & \epsilon_N & 0 & -t & 0\\
0 & -t & 0 & \epsilon_N & 0 & -t\\
0 & 0 & -t & 0 & \epsilon_S & 0\\
0 & 0 & 0 & -t & 0 & \epsilon_S\\
\end{pmatrix},
$$
$$
\Delta = \begin{pmatrix}
0 & \Delta_{S1} & 0 & 0 & 0 & 0\\
-\Delta_{S1} & 0 & 0 & 0 & 0 & 0\\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & \Delta_{S2}\\
0 & 0 & 0 & 0 & -\Delta_{S2} & 0\\
\end{pmatrix}.
$$
So, finally, the full Hamiltonian matrix for the three-site SNS Josephson junction is (drum-roll ...)
$$
H^{SNS}_{BdG} = \begin{pmatrix}
\epsilon_S & 0 & -t & 0 & 0 & 0 & 0 & \Delta_{S1} & 0 & 0 & 0 & 0 \\
0 & \epsilon_S & 0 & -t & 0 & 0 & -\Delta_{S1} & 0 & 0 & 0 & 0 & 0 \\
-t & 0 & \epsilon_N & 0 & -t & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & -t & 0 & \epsilon_N & 0 & -t & 0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & -t & 0 & \epsilon_S & 0 & 0 & 0 & 0 & 0 & 0 & \Delta_{S2} \\
0 & 0 & 0 & -t & 0 & \epsilon_S & 0 & 0 & 0 & 0 & -\Delta_{S2} & 0 \\
0 & -\Delta_{S1}^* & 0 & 0 & 0 & 0 & -\epsilon_S & 0 & t & 0 & 0 & 0 \\
\Delta_{S1}^* & 0 & 0 & 0 & 0 & 0 & 0 & -\epsilon_S & 0 & t & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & t & 0 & -\epsilon_N & 0 & t & 0 \\
0 & 0 & 0 & 0 & 0 & 0 & 0 & t & 0 & -\epsilon_N & 0 & t \\
0 & 0 & 0 & 0 & 0 & -\Delta_{S2}^* & 0 & 0 & t & 0 & -\epsilon_S & 0 \\
0 & 0 & 0 & 0 & \Delta_{S2}^* & 0 & 0 & 0 & 0 & t & 0 & -\epsilon_S
\end{pmatrix}.
$$
After writing the MATLAB script for constructing this Hamiltonians, I also tried to visualise how the matrix looks for a 2-1-2 site SIS junction. MATLAB's syms framework and LaTeX scripting feature for symbolic expressions made it very easy to render this matrix. If interested, have a look here, noting that the symbols in it are a little different from the notation used above (though quite intuitive).
Current calculation
Now that we have built the Hamiltonian of our device, we should be able to calculate some observables of interest. And what better to look into than the device current? So, that is what we will try to figure out here. Also, in keeping up with this blog's spirit of being as explicit as possible, let's do this from the first principles, without skipping any steps!
The current operator for our superconducting device can be defined as the rate of change of the particle "number" operator with respect to time, i.e.
$$ \hat{I} = e \frac{d\hat{n}}{dt}, $$
where the electron charge is included to match the dimensions of the quantities. Using Heisenberg's equation of motion \( \left( d\hat{A}/dt = i/\hbar[H,A] \right) \), we can rewrite the current operator as
$$ \hat{I} = \frac{ie}{\hbar}[\hat{H},\hat{n}]. $$
So, calculating the current operator boils down to calculating the commutator of the superconductor Hamiltonian and the number operator, i.e. \( [\hat{H},\hat{n}] \), and that is what I will show step-by-step below. Let's first define the two operators:
$$ \hat{H} = \sum_{k,\sigma} \epsilon_k \hat{c}_{k\sigma}^\dagger \hat{c}_{k\sigma} + \sum_{k} \left( \Delta_k \hat{c}_{k\downarrow}^\dagger \hat{c}_{-k\uparrow}^\dagger +\Delta_k^* \hat{c}_{-k\uparrow} \hat{c}_{k\downarrow}\right), \quad \hat{n} = \sum_{m,\tau} \hat{c}_{m\tau}^\dagger \hat{c}_{m\tau} $$
To make the calculation managable, I will handle the commutation piece-wise, by taking the terms in the Hamiltonian one at a time. So, let's define \( \hat{H} = \hat{H_1}+\hat{H_2}+\hat{H_3} \), where
$$ \hat{H_1} = \sum_{k,\sigma} \epsilon_k \hat{c}_{k\sigma}^\dagger \hat{c}_{k\sigma}, \quad
\hat{H_2} = \sum_{k}\Delta_k \hat{c}_{k\downarrow}^\dagger \hat{c}_{-k\uparrow}^\dagger, \quad
\hat{H_3} = \sum_{k}\Delta_k^* \hat{c}_{-k\uparrow} \hat{c}_{k\downarrow}, $$
such that,
$$ [\hat{H},\hat{n}] = [\hat{H_1},\hat{n}] + [\hat{H_2},\hat{n}] + [\hat{H_3},\hat{n}] $$
For some of the gory steps below, I will drop the "hats" in the operator notations. Also, before we jump into the calculation, let us recall the fermionic anti-commutation relations which are needed for the algebra here -
$$ \{ c_i,c_j^\dagger\} = c_ic_j^\dagger + c_j^{\dagger}c_i = \delta_{ij}, $$
$$ \{ c_i^\dagger,c_j^\dagger\} = c_i^\dagger c_j^\dagger + c_j^{\dagger}c_i^\dagger = 0, $$
$$ \{ c_i,c_j \} = c_ic_j + c_jc_i = 0. $$
To make it easier to follow the simplifications below, I have highlighted in green the pairs of terms on which I apply an anticommutation relations. With that out of the way, let's start with the first term on the right hand side of (16).
$$ \begin{aligned}\
[H_1,n] &= \sum_{k,m,\sigma,\tau}\epsilon_k \left( c_{k\sigma}^\dagger \textcolor{green}{c_{k\sigma} c_{m\tau}^\dagger} c_{m\tau} - c_{m\tau}^\dagger c_{m\tau}c_{k\sigma}^\dagger c_{k\sigma} \right) \\
&= \sum_{k,m,\sigma,\tau}\epsilon_k \left( c_{k\sigma}^\dagger \left( -c_{m\tau}^\dagger c_{k\sigma} + \delta_{km}\delta_{\sigma\tau} \right) c_{m\tau} - c_{m\tau}^\dagger c_{m\tau}c_{k\sigma}^\dagger c_{k\sigma} \right) \\
&= \sum_{k,m,\sigma,\tau}\epsilon_k \left( \textcolor{green}{-c_{k\sigma}^\dagger c_{m\tau}^\dagger} c_{k\sigma} c_{m\tau} + \delta_{km}\delta_{\sigma\tau} c_{k\sigma}^\dagger c_{m\tau} - c_{m\tau}^\dagger c_{m\tau}c_{k\sigma}^\dagger c_{k\sigma} \right) \\
&= \sum_{k,m,\sigma,\tau}\epsilon_k \left( c_{m\tau}^\dagger c_{k\sigma}^\dagger \textcolor{green}{c_{k\sigma} c_{m\tau}} +
\delta_{km}\delta_{\sigma\tau} c_{k\sigma}^\dagger c_{m\tau} - c_{m\tau}^\dagger c_{m\tau}c_{k\sigma}^\dagger c_{k\sigma} \right) \\
&= \sum_{k,m,\sigma,\tau}\epsilon_k \left( -c_{m\tau}^\dagger \textcolor{green}{c_{k\sigma}^\dagger c_{m\tau}} c_{k\sigma} +
\delta_{km}\delta_{\sigma\tau} c_{k\sigma}^\dagger c_{m\tau} - c_{m\tau}^\dagger c_{m\tau}c_{k\sigma}^\dagger c_{k\sigma} \right) \\
&= \sum_{k,m,\sigma,\tau}\epsilon_k \left( -c_{m\tau}^\dagger \left( -c_{m\tau} c_{k\sigma}^\dagger + \delta_{km}\delta_{\sigma\tau} \right) c_{k\sigma} +
\delta_{km}\delta_{\sigma\tau} c_{k\sigma}^\dagger c_{m\tau} - c_{m\tau}^\dagger c_{m\tau}c_{k\sigma}^\dagger c_{k\sigma} \right) \\
&= \sum_{k,m,\sigma,\tau}\epsilon_k \left( \textcolor{red}{c_{m\tau}^\dagger c_{m\tau} c_{k\sigma}^\dagger c_{k\sigma}} - \delta_{km}\delta_{\sigma\tau} c_{m\tau}^\dagger c_{k\sigma} + \delta_{km}\delta_{\sigma\tau} c_{k\sigma}^\dagger c_{m\tau} - \textcolor{red}{c_{m\tau}^\dagger c_{m\tau}c_{k\sigma}^\dagger c_{k\sigma}} \right) \\
&= \sum_{k,m,\sigma,\tau}\epsilon_k \left( - \delta_{km}\delta_{\sigma\tau} c_{m\tau}^\dagger c_{k\sigma} + \delta_{km}\delta_{\sigma\tau} c_{k\sigma}^\dagger c_{m\tau} \right) \\
&= \sum_{k,\sigma} \epsilon_k \left( - c_{k\sigma}^\dagger c_{k\sigma} + c_{k\sigma}^\dagger c_{k\sigma} \right) \\
\therefore [H_1,n] &= 0
\end{aligned} $$
Neat! There's something satisfying about getting a zero at the end of a long algebriac ordeal. Now let's go for the second term on the right hand side of (16), i.e.
$$ \begin{aligned}\
[H_2,n] &= \sum_{k,m,\tau}\Delta_k \left( c_{k\downarrow}^\dagger \textcolor{green}{c_{-k\uparrow}^\dagger c_{m\tau}^\dagger} c_{m\tau} - c_{m\tau}^\dagger c_{m\tau} c_{k\downarrow}^\dagger c_{-k\uparrow}^\dagger \right) \\
&= \sum_{k,m,\tau}\Delta_k \left( -\textcolor{green}{c_{k\downarrow}^\dagger c_{m\tau}^\dagger} c_{-k\uparrow}^\dagger c_{m\tau} - c_{m\tau}^\dagger c_{m\tau} c_{k\downarrow}^\dagger c_{-k\uparrow}^\dagger \right) \\
&= \sum_{k,m,\tau}\Delta_k \left( c_{m\tau}^\dagger c_{k\downarrow}^\dagger \textcolor{green}{c_{-k\uparrow}^\dagger c_{m\tau}} - c_{m\tau}^\dagger c_{m\tau} c_{k\downarrow}^\dagger c_{-k\uparrow}^\dagger \right) \\
&= \sum_{k,m,\tau}\Delta_k \left( c_{m\tau}^\dagger c_{k\downarrow}^\dagger \left( -c_{m\tau} c_{-k\uparrow}^\dagger + \delta_{-km}\delta_{\uparrow\tau} \right) - c_{m\tau}^\dagger c_{m\tau} c_{k\downarrow}^\dagger c_{-k\uparrow}^\dagger \right) \\
&= \sum_{k,m,\tau}\Delta_k \left( -c_{m\tau}^\dagger \textcolor{green}{c_{k\downarrow}^\dagger c_{m\tau}} c_{-k\uparrow}^\dagger + \delta_{-km}\delta_{\tau\uparrow}c_{m\tau}^\dagger c_{k\downarrow}^\dagger - c_{m\tau}^\dagger c_{m\tau} c_{k\downarrow}^\dagger c_{-k\uparrow}^\dagger \right) \\
&= \sum_{k,m,\tau}\Delta_k \left( -c_{m\tau}^\dagger \left( -c_{m\tau} c_{k\downarrow}^\dagger +\delta_{km}\delta_{\downarrow\tau} \right) c_{-k\uparrow}^\dagger + \delta_{-km}\delta_{\tau\uparrow}c_{m\tau}^\dagger c_{k\downarrow}^\dagger - c_{m\tau}^\dagger c_{m\tau} c_{k\downarrow}^\dagger c_{-k\uparrow}^\dagger \right) \\
&= \sum_{k,m,\tau}\Delta_k \left( \textcolor{red}{c_{m\tau}^\dagger c_{m\tau} c_{k\downarrow}^\dagger c_{-k\uparrow}^\dagger} -\delta_{km}\delta_{\downarrow\tau} c_{m\tau}^\dagger c_{-k\uparrow}^\dagger+ \delta_{-km}\delta_{\tau\uparrow}c_{m\tau}^\dagger c_{k\downarrow}^\dagger - \textcolor{red}{c_{m\tau}^\dagger c_{m\tau} c_{k\downarrow}^\dagger c_{-k\uparrow}^\dagger} \right) \\
&= \sum_{k,m,\tau}\Delta_k \left( -\delta_{km}\delta_{\downarrow\tau} c_{m\tau}^\dagger c_{-k\uparrow}^\dagger+ \delta_{-km}\delta_{\tau\uparrow}c_{m\tau}^\dagger c_{k\downarrow}^\dagger \right) \\
&= \sum_{k}\Delta_k \left( -c_{k\downarrow}^\dagger c_{-k\uparrow}^\dagger + \textcolor{green}{c_{-k\uparrow}^\dagger c_{k\downarrow}^\dagger} \right) \\
&= \sum_{k}\Delta_k \left( -2c_{k\downarrow}^\dagger c_{-k\uparrow}^\dagger \right) \\
\therefore [H_2,n] &= -2H_2
\end{aligned} $$
Then, solving the third term in the right hand side of (16) quite similarly, we get \( [H_3,n] = 2H_3 \). So, we can finally write the commutator of the Hamiltonian and the number operator as
$$ \begin{align}
[\hat{H},\hat{n}] &= -2\hat{H_2}+2\hat{H_3} \\
&= -2\sum_{k} \left( \Delta_k \hat{c}_{k\downarrow}^\dagger \hat{c}_{-k\uparrow}^\dagger -\Delta_k^* \hat{c}_{-k\uparrow} \hat{c}_{k\downarrow}\right)
\end{align}$$
Expressing the superconducting order parameter as \( \Delta e^{i\phi} \) in (14), we can see that
$$ \frac{d\hat{H}}{d\phi} = i\sum_{k}\left( \Delta e^{i\phi} \hat{c}_{k\downarrow}^\dagger \hat{c}_{-k\uparrow}^\dagger -\Delta e^{-i\phi} \hat{c}_{-k\uparrow} \hat{c}_{k\downarrow}\right)
$$
So, from (23) and (24), we find that
$$ [\hat{H},\hat{n}] = 2i \frac{d\hat{H}}{d\phi}
$$
Consequently, we finally have the expression for the current operator that we've been seeking (another drum-roll ...)
$$ \boxed{
\hat{I} = -\frac{2e}{\hbar} \frac{d\hat{H}}{d\phi}
} $$
Quite interesting to see that the current turned out to be the differential of the Hamiltonian with respect to the superconducting phase parameter! Also, the factor of \( 2e \) accounts well for the fact that we are looking at a current carried by Cooper pairs. With the current operator well defined, we can calculate the current (its expectation value) using
$$ I = \expval{\hat{I}}{\Psi}. $$
This can be simplified by working some tricks on Schrödinger's eigenvalue equation, viz. differentiating it with respect to \( \phi \)
$$ \begin{aligned}
\frac{d}{d\phi} \hat{H}\ket{\Psi_n} &= \frac{d}{d\phi} E_n\ket{\Psi_n} \\
\therefore \frac{d\hat{H}}{d\phi}\ket{\Psi_n} + \hat{H}\frac{d\ket{\Psi_n}}{d\phi} &= \frac{dE_n}{d\phi}\ket{\Psi_n} + E_n\frac{d\ket{\Psi_n}}{d\phi}
\end{aligned}
$$
Now, taking the inner product with \( \bra{\Psi_n} \) and using the orthogonality of the eigenvectors we get,
$$ \begin{aligned}
\bra{\Psi_n} \frac{d\hat{H}}{d\phi}\ket{\Psi_n} + \bra{\Psi_n} \hat{H}\frac{d\ket{\Psi_n}}{d\phi} &= \bra{\Psi_n} \frac{dE_n}{d\phi}\ket{\Psi_n} + \bra{\Psi_n} E_n\frac{d\ket{\Psi_n}}{d\phi} \\
\therefore \bra{\Psi_n} \frac{d\hat{H}}{d\phi}\ket{\Psi_n} + \textcolor{red}{E_n\bra{\Psi_n}\frac{d\ket{\Psi_n}}{d\phi}} &= \braket{\Psi_n}{\Psi_n} \frac{dE_n}{d\phi} + \textcolor{red}{\bra{\Psi_n} E_n\frac{d\ket{\Psi_n}}{d\phi}} \\
\therefore \bra{\Psi_n} \frac{d\hat{H}}{d\phi}\ket{\Psi_n} &= \frac{dE_n}{d\phi}
\end{aligned}
$$
Thus, using (26), (27) and (29), the current through the device can be calculated as
$$ \boxed{
\begin{aligned}
I = -\frac{2e}{\hbar} \sum_{n}\frac{dE_n}{d\phi}
\end{aligned}
}
$$
MATLAB/Python scripts for solving the Hamiltonian to calculate the Andreev bound states spectrum and current-phase relation (CPR) of the junctions will soon be posted, along with some example results. Next, I will also delve into replacing the normal metal/insulator with a ferromagnet, and perhaps reproduce the 0-\( \pi \) transitions expected for the SFS Josephson junctions.
In the meantime, feedback and questions are welcome!