  We present efficient methods for Brillouin zone integration with a non-zero
but possibly very small broadening factor $\eta$, focusing on cases in which
downfolded Hamiltonians can be evaluated efficiently using Wannier
interpolation. We describe robust, high-order accurate algorithms automating
convergence to a user-specified error tolerance $\varepsilon$, emphasizing an
efficient computational scaling with respect to $\eta$. After analyzing the
standard equispaced integration method, applicable in the case of large
broadening, we describe a simple iterated adaptive integration algorithm
effective in the small $\eta$ regime. Its computational cost scales as
$\mathcal{O}(\log^3(\eta^{-1}))$ as $\eta \to 0^+$ in three dimensions, as
opposed to $\mathcal{O}(\eta^{-3})$ for equispaced integration. We argue that,
by contrast, tree-based adaptive integration methods scale only as
$\mathcal{O}(\log(\eta^{-1})/\eta^{2})$ for typical Brillouin zone integrals.
In addition to its favorable scaling, the iterated adaptive algorithm is
straightforward to implement, particularly for integration on the irreducible
Brillouin zone, for which it avoids the tetrahedral meshes required for
tree-based schemes. We illustrate the algorithms by calculating the spectral
function of SrVO$_3$ with broadening on the meV scale.
