I've been trying to rederive (13.37) from Ashcroft & Mermin:
$$ \sigma_{xy}(\omega+i\eta) = \frac{-ie^2}{\hbar^2(\omega + i\eta)} \sum_{n\neq n'} (f) \left( \frac{\langle n\mathbf k|\nabla_x|n'\mathbf k \rangle \langle n'\mathbf k|\nabla_y|n\mathbf k \rangle}{\hbar\omega + i\eta + \varepsilon_{n\mathbf k} - \varepsilon_{n'\mathbf k}} + \frac{\langle n\mathbf k|\nabla_y|n'\mathbf k \rangle \langle n'\mathbf k|\nabla_x|n\mathbf k \rangle}{-\hbar\omega - i\eta + \varepsilon_{n\mathbf k} - \varepsilon_{n'\mathbf k}} \right).$$
So I understand (pretty much) all derivation steps, but at some point the authors mention the use of time-dependent pertubation theory, p. 252, between (13.36) and (13.37).
I'm assuming the same applies as in appendix E (pp. 765-766), where a Hamiltonian about $\mathbf k + \mathbf q$ is considered, $$ H_{\mathbf k+\mathbf q} = H_{\mathbf k} + \frac{\hbar^2}m \mathbf q \cdot\left(\mathbf k + \frac\nabla i\right) + \frac{\hbar^2}{2m}q^2 $$ where the second and third terms are treated as pertubation $H_1$, and the eigenvalues can correspondingly be $$ \varepsilon_{\mathbf k+\mathbf q} = \varepsilon_{\mathbf k} + \sum_i \frac{\partial\varepsilon}{\partial k_i}q_i + \frac12 \sum_{ij} \frac{\partial^2\varepsilon}{\partial k_i\partial k_j} q_i q_j+ \dots $$
So far so good. I'm assuming I again want to match the orders of q. So since equation (13.34) requires the second order derivative (factors of q of order 2), I will only consider those in the pertubation.
Now I'm getting into territory where I'm not sure I understand. I have $$ A_0(t) = 1$$ and the other coefficient zero. First order approximation gives $$ A_1(t) = 1\quad\text{and}\quad B_1(t) = -\frac i\hbar\int \langle u_{n\mathbf k}| H_1 | u_{n'\mathbf k}\rangle e^{i\frac{\varepsilon_{n\mathbf k} -\varepsilon_{n'\mathbf k}}{\hbar}t}dt$$
This yields three terms for $B_1$, of which two I believe are zero (since the Bloch wave functions are orthogonal), leaving me with $$ B_1(t) = -q\frac{\hbar^2}{m}\frac{\langle u_{n\mathbf k}| \frac{\nabla}i | u_{n'\mathbf k}\rangle}{\varepsilon_{n\mathbf k} - \varepsilon_{n'\mathbf k}} \left[1-\exp\left(i\frac{\varepsilon_{n\mathbf k}-\varepsilon_{n'\mathbf k}}\hbar t\right)\right] $$ Repeating the procedure again, I get $$ A_2(t) = q^2\left(\frac{\hbar^2}{m}\right)^2\frac{\langle u_{n\mathbf k}| \frac{\nabla}i | u_{n'\mathbf k}\rangle\langle u_{n'\mathbf k}| \frac{\nabla}i | u_{n\mathbf k}\rangle}{\varepsilon_{n\mathbf k} - \varepsilon_{n'\mathbf k}} \int \left[\exp\left({i\frac{\varepsilon_{n\mathbf k} - \varepsilon_{n'\mathbf k}}{\hbar} t}\right)-\exp\left(2i\frac{\varepsilon_{n\mathbf k}-\varepsilon_{n'\mathbf k}}\hbar t\right)\right] dt $$ for the $q^2$ term, thus $$ A_2(t) = q^2\left(\frac{\hbar^2}{m}\right)^2\frac{\langle u_{n\mathbf k}| \frac{\nabla}i | u_{n'\mathbf k}\rangle\langle u_{n'\mathbf k}| \frac{\nabla}i | u_{n\mathbf k}\rangle}{\varepsilon_{n\mathbf k} - \varepsilon_{n'\mathbf k}} \left[\frac1{\varepsilon_{n\mathbf k}-\varepsilon_{n'\mathbf k}} - \frac2{\varepsilon_{n\mathbf k}-\varepsilon_{n'\mathbf k}}\right], $$ but I don't think this is correct. The expressions in the multi-level system don't differ very much, so I think the fault in my derivation is more fundamental, but I can't find it.
