  We consider a general class of symmetric or Hermitian random band matrices
$H=(h_{xy})_{x,y \in \llbracket 1,N\rrbracket^d}$ in any dimension $d\ge 1$,
where the entries are independent, centered random variables with variances
$s_{xy}=\mathbb E|h_{xy}|^2$. We assume that $s_{xy}$ vanishes if $|x-y|$
exceeds the band width $W$, and we are interested in the mesoscopic scale with
$1\ll W\ll N$. Define the {\it{generalized resolvent}} of $H$ as $G(H,Z):=(H -
Z)^{-1}$, where $Z$ is a deterministic diagonal matrix with entries $Z_{xx}\in
\mathbb C_+$ for all $x$. Then we establish a precise high-probability bound on
certain averages of polynomials of the resolvent entries. As an application of
this fluctuation averaging result, we give a self-contained proof for the
delocalization of random band matrices in dimensions $d\ge 2$. More precisely,
for any fixed $d\ge 2$, we prove that the bulk eigenvectors of $H$ are
delocalized in certain averaged sense if $N\le W^{1+\frac{d}{2}}$. This
improves the corresponding results in \cite{HeMa2018} under the assumption
$N\ll W^{1+\frac{d}{d+1}}$, and in \cite{ErdKno2013,ErdKno2011} under the
assumption $N\ll W^{1+\frac{d}{6}}$. For 1D random band matrices, our
fluctuation averaging result was used in \cite{PartII,PartI} to prove the
delocalization conjecture and bulk universality for random band matrices with
$N\ll W^{4/3}$.
