in the appendix that for an operator $O$ in the discrete DoFs that

$$d_t \langle \psi | O | \psi \rangle = 2E_a  \mathrm{Re} ( \eta_i^* \langle \epsilon_i|\partial_{k_a} \epsilon_j \rangle \langle \epsilon_j |O|\epsilon_k \rangle \eta_k$$

It is straightforward to show that this expression is $-i \langle \psi [ O, H_{MF}] |\psi\rangle$.

Recovery of the Kubo formula
----------------------------

Treating $H_{MF}$ as a perturbation to $H_0$ and applying the standard non-degenerate time-independent perturbation theory to the $i$th eigenstate of $H_0$, $|\epsilon_i\rangle$, the first order correction to $|\epsilon_i\rangle$ which we denote as $|\epsilon_i^{(1)}\rangle$ reads $$|\epsilon_i^{(1)}\rangle = \sum_j |\epsilon_j \rangle \frac{ \langle \epsilon_j|H_{MF}|\epsilon_i\rangle}{E_i - E_j}$$ so that to the correction to the expectation value of an observable $O$ for state $|\epsilon_i\rangle$ to first order in $\vec{E}$, $\delta \langle i| O | i \rangle$ is $$\begin{aligned}
    && \delta \langle i| O|i \rangle  \\ \nonumber
    &=& 2\mathrm{Re} (\langle \epsilon_i|O|\epsilon_i^{(1)}\rangle ) \\ \nonumber
    &=& 2\mathrm{Re} \sum_j \frac{  \langle \epsilon_i |O| \epsilon_j \rangle \langle \epsilon_j |H_{MF}|\epsilon_i\rangle}{E_i - E_j}. \label{c1oe}\end{aligned}$$

However, since $$\begin{aligned}
    && \partial_k \langle \epsilon_i | H_0 |\epsilon_i \rangle = 0 \\
    &\Rightarrow& \langle \partial_k \epsilon_i |\epsilon_j \rangle (E_i-E_j) = \langle \epsilon_i|\partial_k H_0 |\epsilon_j  \rangle \\
    &\Rightarrow& \langle \partial_k \epsilon_i|\epsilon_j \rangle = \frac{ \langle \epsilon_i| \partial_k H_0 