  We consider a discrete elliptic equation on the $d$-dimensional lattice
$\mathbb{Z}^d$ with random coefficients $A$ of the simplest type: they are
identically distributed and independent from edge to edge. On scales large
w.r.t. the lattice spacing (i.e., unity), the solution operator is known to
behave like the solution operator of a (continuous) elliptic equation with
constant deterministic coefficients. This symmetric ``homogenized'' matrix
$A_{\mathrm {hom}}=a_{\mathrm {hom}}\operatorname {Id}$ is characterized by
$\xi\cdot A_{\mathrm {hom}}\xi=\langle(\xi+\nabla\phi )\cdot
A(\xi+\nabla\phi)\rangle$ for any direction $\xi\in\mathbb {R}^d$, where the
random field $\phi$ (the ``corrector'') is the unique solution of
$-\nabla^*\cdot A(\xi+\nabla\phi)=0$ such that $\phi(0)=0$, $\nabla\phi$ is
stationary and $\langle\nabla\phi\rangle=0$, $\langle\cdot\rangle$ denoting the
ensemble average (or expectation). It is known (``by ergodicity'') that the
above ensemble average of the energy density $\mathcal
{E}=(\xi+\nabla\phi)\cdot A(\xi+\nabla\phi)$, which is a stationary random
field, can be recovered by a system average. We quantify this by proving that
the variance of a spatial average of $\mathcal {E}$ on length scales $L$
satisfies the optimal estimate, that is, $\operatorname {var}[\sum \mathcal
{E}\eta_L]\lesssim L^{-d}$, where the averaging function [i.e., $\sum\eta_L=1$,
$\operatorname {supp}(\eta_L)\subset\{|x|\le L\}$] has to be smooth in the
sense that $|\nabla\eta_L|\lesssim L^{-1-d}$. In two space dimensions (i.e.,
$d=2$), there is a logarithmic correction. This estimate is optimal since it
shows that smooth averages of the energy density $\mathcal {E}$ decay in $L$ as
if $\mathcal {E}$ would be independent from edge to edge (which it is not for
$d>1$). This result is of practical significance, since it allows to estimate
the dominant error when numerically computing $a_{\mathrm {hom}}$.
