I am trying to derive the result of the TKNN formula but am experiencing difficulty in deriving the Kubo formula. The Kubo formula used in the TKNN paper is, $$ \sigma_{xy}= \frac{ie^2}{\hbar} \sum_{E^a < E_F < E^b} \frac{\langle a|v_x|b \rangle \langle b|v_y|a \rangle - \langle a|v_y|b \rangle \langle b|v_x|a \rangle}{(E^a - E^b)^2} .$$
The following is my work so far. First of all, from linear response theory, for a general observable, $$ O=\langle\hat{O}(t)\rangle-i\int_{t_0}^{t}\mathrm{d} t' \langle[\hat{O}(t),H_{\mathrm{ext}}(t')]\rangle $$ where $H_{\mathrm{ext}}$ is the perturbing Hamiltonian. Since $H_{\mathrm{ext}}=-\int J\cdot A(r,t)\mathrm{d}^3r$, if we consider current density we have $$ \vec{J}_i(r,t)=\langle j_i(t)\rangle+i\int_{-\infty}^{t}\mathrm{d} t'\int\mathrm{d}^3{r'} \langle[j_i(r,t),j_j(r',t')]\rangle A_j(r',t'). $$ Ignoring the first term, we can rewrite this as $$ \vec{J}_i(r,t)=-\int_{-\infty}^\infty\mathrm{d}t\int \mathrm{d}^3r R_{ij}(r-r',t-t')A_j(r',t') $$ where $$ R_{ij}(r-r',t-t')=-i\theta(t-t')\langle[j_i(r,t),j_j(r',t')]\rangle. $$ We can assume that $A(r,t)=A(r)e^{-i\omega t}$ to show that $E(r,t)=i\omega A(r,t)$. Then by Fourier transforming the expression for the expectation value of the current density, we get $$ \vec{J}(k,\omega)=-\frac{R_{ij}(k,\omega)}{i\omega}E_j(k,\omega) $$ which means DC conductivity is $$ \sigma_{ij}=\lim_{\omega\rightarrow0}\frac{iR_{ij}(k,\omega)}{\omega} $$ The Fourier transform of the response function can be shown to be, $$ R_{ij}(k,\omega)=-i\int_0^\infty \mathrm{d} t e^{i\omega t}\langle[j_i(k,t),j_j(-k,0)]\rangle $$ If we evaluate expectation values using the grand canonical ensemble (ie $\langle \hat{O}\rangle = Tr(e^{-\beta H} \hat{O})/Z$ where $H$ includes the chem. pot.) and integrate over time and perform a first order $\omega$ expansion we end up with $$ R_{ij}(k,\omega)=\sum_{n,m}(e^{-E_m \beta}-e^{-E_n\beta})\left(\frac{\langle n|j_i(k,\omega)|m\rangle\langle m|j_j(-k,\omega)|m\rangle}{E_n-E_m+i\epsilon}-\omega\frac{\langle n|j_i(k,\omega)|m\rangle\langle m|j_j(-k,\omega)|m\rangle}{(E_n-E_m+i\epsilon)^2}\right)/Z $$ This is from where I don't know how to proceed. All derivations for the Kubo formula I have found on this site and other resources online do not consider the position dependence of the current density operator. Any advice that will point me in the right direction will be appreciated.