  This work is concerned with fractional Gaussian fields, i.e. Gaussian fields
whose covariance operator is given by the inverse fractional Laplacian
$(-\Delta)^{-s}$ (where, in particular, we include the case $s >1$). We define
a lattice discretization of these fields and show that their scaling limits --
with respect to the optimal Besov space topology (up to an endpoint case) --
are the original continuous fields. As a byproduct, in dimension $d<2s$, we
prove the convergence in distribution of the maximum of the fields. A key tool
in the proof is a sharp error estimate for the natural finite difference scheme
for $(-\Delta)^s$ under minimal regularity assumptions, which is also of
independent interest.
