In section 18.7 of Bruus & Flensberg the authors provide a microscopic derivation of the Josephson effect.
The hamiltonian on both sides of the tunnelling junction is just the typical BCS hamiltonian, on one side (with fermion operators $c$) \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} and on the other side (with fermion operators $f$) \begin{equation} H_f = \sum_{k,\sigma} \epsilon_k f_{k,\sigma}^\dagger f_{k,\sigma} - \sum_{k}(\Delta e^{-i\phi_c}f^\dagger_{k,\uparrow}f^\dagger_{-k,\downarrow}+\Delta e^{i\phi_c} f_{-k,\downarrow}f_{k,\uparrow}) \end{equation} Here we let the gap parameter for both superconductors have the same magnitude $\Delta$ but different phases $\phi_c$ and $\phi_f$. We then introduce a tunnelling Hamiltonian coupling the two superconductors \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}
We can deal with the phases by introducing a gauge transformation $c\rightarrow e^{-i\phi_c/2} c$ and $f \rightarrow e^{-i\phi_f/2} f$ so that the tunnelling coefficients acquire a phase $t \rightarrow e^{-i\phi/2} t$ with $\phi=\phi_c-\phi_f$. Then we see that the Josephson current is \begin{equation} I_J = \langle I \rangle = -2e \langle \dot{N} \rangle = -2e \bigg\langle \frac{\partial H_t}{\partial \phi}\bigg\rangle \end{equation} where we used the Heisenberg equations of motion $\dot{N} = -i[N,H]$ and the fact that $N=-i\frac{\partial}{\partial \phi}$. The current can be calculated perturbatively using the Dyson series \begin{equation} \exp(-\beta H) = \exp(-\beta H_0)\bigg[1-\int_0^\beta d\tau \ \hat{H}_t(\tau)\bigg]+o((H_t)^2) \end{equation} where $\hat{H}(\tau)$ is in the interaction picture. Then \begin{equation} I_J = -2e\bigg[\bigg\langle \frac{\partial H_t}{\partial \phi}\bigg\rangle_0-\frac{1}{2}\bigg\langle \int_0^\beta d\tau \ \hat{H}_t(\tau) H_t\bigg\rangle_0+\ ...\bigg] \end{equation} where $\langle \rangle_0$ denotes thermal averaging with respect to $H_0=H_c+H_f$. The first order contribution is \begin{equation} \bigg\langle \frac{\partial H_t}{\partial \phi}\bigg\rangle_0 = \frac{\partial}{\partial \phi}\text{Tr}(e^{-\beta H_0} H_t) = \frac{\partial}{\partial \phi} \sum_m e^{-\beta E^0_m} \langle m|H_t|m\rangle \end{equation} which I presume (but I am really not too sure about this) vanishes since $H_t$ moves an odd number of fermions from one superconductor to the other, so its action on any eigenstate $|m\rangle$ of the BCS hamiltonian will produce a state orthogonal to it. The second order contribution on the other hand is given as \begin{equation} \frac{\partial}{\partial \phi} \int_0^\beta d\tau \ (t^2 e^{i\phi} \mathcal{F}_{\downarrow \uparrow}(\textbf{k},\tau) \mathcal{F}^*_{\downarrow \uparrow}(\textbf{p},-\tau)+c.c.) \end{equation} where we defined $\mathcal{F}(\textbf{k},\tau)=-\langle \mathcal{T} c^\dagger_{k,\downarrow}(\tau)c^\dagger_{-k,\uparrow}(0)\rangle$. However I have a hard time making sense of this equation. I find that \begin{equation} \frac{1}{2}\frac{\partial}{\partial \phi} \int_0^\beta d\tau \ \bigg \langle \sum_{k,p,\sigma}\sum_{k',p',\sigma'} (tc^\dagger_{k\sigma} f_{p,\sigma}+t^* f^\dagger_{p\sigma} c_{k,\sigma})(tc^\dagger_{k'\sigma'} f_{p',\sigma'}+t^* f^\dagger_{p'\sigma'} c_{k',\sigma'})\bigg \rangle \end{equation} but I'm not sure how this simplifies to the result by Bruus & Flensberg.
