Current operator, why is this form valid up to second order in $q$?
Context
In <PERSON> textbook Topological insulators and topological superconductors, an approximate form of the current operator in momentum space is derived. It is said to be valid to second order in momentum, but I don´t understand why.
They start with a generic lattice Hamiltonian with translational invariance, so that it can be written in momentum space:
$$ H=\sum_{ij} c^\dagger_{i}h_{ij}c_{j}, \quad H=\sum_{\mathbf{k}}c_{\mathbf{k}}^\dagger h_\mathbf{k} c_{\mathbf{k}}.\tag{1}$$ Here, $i,j$ correspond to the lattice sites and $\mathbf{k}$ to momentum.
To obtain the current, they use the continuity equation in momentum space: $$ \dot{\rho}(\mathbf{x})+\nabla\cdot \mathbf{J}(\mathbf{x})=0 \implies \dot{\rho}\mathbf{q}-i\mathbf{q}\cdot\mathbf{J}\mathbf{q}=0\tag{2} .$$
Next, they use the density operator in momentum space and the <PERSON> equation of motion to write: $$ \rho_\mathbf{q}=\frac{1}{\sqrt{N}}\sum_\mathbf{k}c^\dagger_{\mathbf{k+q}}c_\mathbf{k}\implies -i\mathbf{q}\cdot\mathbf{J}\mathbf{q}=-\dot{\rho}\mathbf{q}=i[\rho_\mathbf{q},H], \tag{3} $$ which after some algebra becomes $$ \mathbf{q}\cdot \mathbf{J}\mathbf{q}=\frac{1}{\sqrt{N}}\sum\mathbf{k} (h_{\mathbf{k}+\mathbf{q}}-h_\mathbf{k})c_{\mathbf{k}+\mathbf{q}}^\dagger c_\mathbf{k}. \tag{4}$$
Core of the question
After using $h_{\mathbf{k}+\mathbf{q}}-h_\mathbf{k} \approx \mathbf{q}\cdot \partial_\mathbf{k} h_\mathbf{k}$, it is possible to identify the current operator by comparing both sides of equation (4) as: $$ \mathbf{J}\mathbf{q}=\frac{1}{\sqrt{N}}\sum\mathbf{k} \frac{h_\mathbf{k}}{\partial \mathbf{k}} c^\dagger_{\mathbf{k}+\mathbf{q}}c_\mathbf{k}$$
I understand everything up to this point. Clearly, this form of the current operator is valid to first order in $\mathbf{q}$. 281