Inconsistency in microscopic derivation of the Josephson effect?

Question

My question is related to, but different from, other questions in this platform (see this, this and this).

One of the few places I have found in literature that provide a microscopic derivation of the Josephson Effect is in section 18.7 of Bruus & Flensberg. In the book, the authors consider two superconducting leads at different phases $\phi_c$ and $\phi_f$. The unperturbed Hamiltonian for the lead $c$ is \begin{equation} H_c = \sum_{k,\sigma} \epsilon_k c_{k\sigma}^\dagger c_{k\sigma} - \sum_{k}(\Delta e^{-i\phi_c}c^\dagger_{k\uparrow}c^\dagger_{-k\downarrow}+\Delta e^{i\phi_c} c_{-k\downarrow}c_{k\uparrow}), \end{equation} with an equivalent one holding for the lead $f$. They then introduce term to allow for fermion tunneling between the leads \begin{equation} H_t = \sum_{k,p,\sigma} (tc^\dagger_{k\sigma} f_{p\sigma}+t^* f^\dagger_{p\sigma} c_{k\sigma}). \end{equation}

They then calculate the current through the leads as \begin{equation} I_J = \langle I \rangle = -e \langle \dot{N}_c \rangle = -2e \bigg\langle \frac{\partial H_t}{\partial \phi}\bigg\rangle. \end{equation}

As I understand, to get the last equation one uses the Heisenberg equation of motion $\dot{N}_c = -i[N_c,H]$, where $N_c$ is the occupation number operator in the superconducting lead $c$, followed by the use of the identification $N_c - N_f = -i \partial_\phi$ together with some constraint $N_c + N_f = N \in \mathbb{Z}$.

This already I find quite problematic. I have read in some sources that the identification $N = - i \partial_\phi$ is mathematically inconsistent. So, how can we trust any results that rely on it? In particular, what justifies the use of this identity here?

Even if we ignore this problematic step, there is another detail I would like some clarification on.

Let us start from the second equality in the current equation. We use the Heisenberg equation and arrive at \begin{equation} I_J = i e \langle [N_c,H_c + H_f + H_t] \rangle. \end{equation}

The expectation value above needs to be evaluated according to some state (or distribution). The authors do not specify what it should be, so I think the simplest assumption is that the state are thermal states of the unperturbed Hamiltonian $H_0 = H_c + H_f$. At zero temperature, this reduces to the BCS ground state $H_0 |GS(\phi_c, \phi_f)\rangle = E_0 |GS(\phi_c, \phi_f)\rangle$, where \begin{equation} |GS(\phi_c, \phi_f)\rangle = \tfrac{1}{\mathcal{N}} e^{\Lambda_c^\dagger(\phi_c)} e^{\Lambda_f^\dagger(\phi_f)} | 0 \rangle, \end{equation} where the exponentials create BCS condensates on each lead separately. With this assumtpion, the $H_c + H_f$ term in the commutator vanishes in the expression for the current, since it acts on the ground state to become a number. We are then left with having to calculate $[N_c,H_t]$. After some dull but straightforward calculation, we find \begin{align} [N_c,H_t] &= \sum_{k, p, \sigma} \sum_{q, \tau} [c^\dagger_{q\tau} c_{q \tau}, (tc^\dagger_{k\sigma} f_{p,\sigma}+t^* f^\dagger_{p\sigma} c_{k,\sigma})] \\ &= \sum_{k, p, \sigma} \sum_{q, \tau} \delta_{q\tau, k\sigma} \left( t c_{q\tau}^\dagger f_{p\sigma} - t^* f_{p\sigma}^\dagger c_{q\tau} \right) \\ &= \sum_{k, p, \sigma} \left( t c_{k\sigma}^\dagger f_{p\sigma} - t^* f_{p\sigma}^\dagger c_{k\sigma} \right). \end{align} Here comes the issue: although $|GS(\phi_c, \phi_f)\rangle$ has an ill-defined particle number, it should have a definite fermion parity. In fact, each condensate $c,f$ should have a well defined fermion parity. Moreover, the operator from the expression above flips the fermion parity in each condensate (since it hops a single fermion from one region to the other). This means that the state \begin{equation} [N_c,H_t]|GS(\phi_c, \phi_f)\rangle \end{equation} must have the opposite fermion parity in each condensate, which means it must be orthogonal to $|GS(\phi_c, \phi_f)\rangle$ itself. In other words, this calculation seems to imply \begin{equation} I_J = i e \langle GS(\phi_c, \phi_f) | [N_c,H_t]|GS(\phi_c, \phi_f)\rangle = 0. \end{equation}

The current must be identically zero! Even though the calcultion through the original method seemed to imply otherwise.

How can this apparent paradox be resolved?

Answer 1

Let us analyze this problem in the limit $|t|\ll |\Delta |$, the superconducting gap. The tunneling term that you wrote down tunnels individual electrons and (as you clearly realize) excites each superconducting island above the gap. This is not the relevant term if you want to see the Josephson current.

Instead, we want the low-energy process -- Cooper pair hopping. This is generated by the tunneling term you wrote down in second order perturbation theory with strength $\sim t^2 / \Delta$. You can check that this will give you non-zero $\langle \dot{N}\rangle$, as there is no parity violation.

This may be confusing because I use time-independent perturbation theory to describe dynamics. Here is an argument which hopefully helps you reconcile this:

You can think of similar dynamics in a simpler case. Consider the Hamiltonian of a single particle in three orbitals $H = -t(c^\dagger_1 c_2 + c^\dagger_2 c_3 + H.c.) + U c^\dagger_2 c_2$ (with $U\gg t$), and start with the electron in the wavefunction $(1/\sqrt{2},0,i/\sqrt{2})$ and looking at the current from $1\to3$. I recommend plotting the three probability amplitudes on Mathematica to get a feel for it! At short times $T \sim 1/t \equiv T_1$, you will see that your result is correct; the initial $\langle\dot{N}_3 - \dot{N}_1\rangle$ is indeed zero. But if you plot it for very long times $\sim U/t^2 \equiv T_2$ (that is, you coarse grain and study dynamics over timescales $\gg T_1$), you see that there is a non-zero current from site 1 to 3 at $T=0$. Try to guess an effective Hamiltonian with only sites 1 and 3 that describes this dynamics!

---

Edit: I am adding the plots for this simple model from Mathematica for clarity. I am plotting the probability for the particle to be at site 1 as a function of time for $t = 1, U = 10$.

Now, imagine the experimentalist does not have the ability to resolve timescales at the order of $1/t$ but only sees the coarse-grained charge (as a function of time) in their measurements. What would they conclude is the current at time $=0$? The point is that this is non-zero even though the true current in the model is zero, and this non-zero value is estimated correctly by using the two-site effective model you get from second order perturbation theory $$ H = -\frac{t^2}{U}(c^\dagger_1 c_3 + \textrm{h.c.}) $$

---

In short: if you are being precise the Hamiltonian that you wrote down you are right that there is no current between the two superconducting islands at time $T=0$, which makes sense precisely because of your argument. But if you care about long-time behavior and look at the coarse grained $N_c(T)-N_f(T)$, you will see that there is current flowing at $T=0$ described by an effective Hamiltonian containing a Cooper pair tunneling term. In the experiments we observe the coarse grained current as a function of time.