Correct way of computing the conductivity from Kubo's formula in superconductivity Consider a full Hamiltonian $H=H_0 +H_1$ on a lattice with hopping terms $H_0$ of the form $t c^* c$ and attractive interaction $H_1= -V c^* c^* c c$ so that superconductivity can occur. We can then approximate $H$ with the quasi-free BCS Hamiltonian $H_\text{BCS}$ of the form $\xi c^* c + \Delta^* cc+ \Delta c^* c^*$ where $\xi, \Delta$ can be found using the BdG equations.
In order to calculate the conductivity $\sigma$ of such system in the frameworks of linear response theory (Kubo formalism), we need to calculate the current-current correlation $\chi_{ab}(t)$, given as
$$
\chi_{ab} (t) = -i\theta(t) \langle [J_a(t),J_b]\rangle
$$
I can think of 2 ways to compute $\chi_{ab}$. One way is to approximate the time evolution $J_a(t) =e^{iHt} J_a e^{-iHt}$ by substituting $H \mapsto H_\text{BCS}$ and similarly substitue $\langle \cdots \rangle \mapsto \langle \cdots \rangle_{\text{BCS}}$. Another way would be to calculate $\chi_{ab} (t)$ perturbatively using one-loop corrections, etc.
Which way would be the "more correct" way of doing such calculations? I have seen Schrieffer use the first method, but also seen papers using the latter. What are some of the reasons towards one or the other?
EDIT: To be more concrete, in the one-loop corrections, we use $\langle\cdots \rangle \mapsto \langle \cdots \rangle_0$ and apply Wick's theorem based on $\langle \cdots \rangle_0$. I'm wondering if we can circumvent the problem of calculating a lot of diagrams by simply calculating
$$
\chi_{ab} (t) = -i\theta(t) \langle [e^{iH_\text{BCS} t} J_a e^{-iH_\text{BCS} t} ,J_b]\rangle_\text{BCS}
$$
 A: Let me try to provide a sketch of derivation. Consider the normal metal, so interaction with external field is
$$H=\int d^3r\psi^{\dagger}\epsilon(p-eA)\psi,$$
where $\epsilon$ is dispersion law. Performing variation with respect to $A$, one find the expression for current,
$$j=e\psi^{\dagger}\nabla_p\epsilon(p)\psi-e^2\psi^{\dagger}\nabla_p[(A\nabla_p)\epsilon(p)]\psi. (*)$$
In order to find the contribution from the first term, we should consider
$$H_{\text{int}}=-\int d^3j(r)A(r),$$
which is the first order approximation (in $A$),

and it can be computed in terms of Matsubara Green functions followed by analytic continuation. The key point is the resulting expression has form (in momentum space)
$$j(k)=Q(k)A(k),$$
where for normal metal the gauge invariance guarantees that $Q(0)=0$ (cf. $(*)$).
In case of superconductor, we should consider not only normal averages (as in case of normal metal) but anomalous, too. It means that when one performs Wick contraction, the anomalous contractions should be taken into account, too. A little bit tedious derivation says that $Q(0)\neq 0$.
Finally, for me it seems that both derivations are the same, because linear-responce is nothing more than 1-loop approximation.
