  We investigate the large-distance asymptotics of optimal Hardy weights on
$\mathbb Z^d$, $d\geq 3$, via the super solution construction. For the free
discrete Laplacian, the Hardy weight asymptotic is the familiar
$\frac{(d-2)^2}{4}|x|^{-2}$ as $|x|\to\infty$. Here we prove that the
inverse-square behavior of the optimal Hardy weight is robust for general
elliptic coefficients on $\mathbb Z^d$. Namely, (1) annular averages have
inverse-square scaling in general, (2), for ergodic coefficients there is an
almost-sure pointwise inverse-square upper bound, and (3), for weakly random
i.i.d. coefficients there is a probabilistic inverse-square lower bound. The
results have continuum analogs and also yield $|x|^{-4}$-scaling for Rellich
weights on $\mathbb Z^d$. The technique relies on Green's function estimates
rooted in homogenization theory. Along the way, we identify for the first time
an asymptotic expansion of the averaged Green's function and of its derivatives
up to order $d+1$ at small ellipticity contrast. This leverages a recent
perturbative approach to homogenization theory developed by Bourgain. The
result implies rather universal asymptotics of the random Green's function via
concentration bounds of Marahrens-Otto.
