This work is concerned with fractional Gaussian fields, i.e. Gaussian fields whose covariance operator is given by the inverse fractional Laplacian $(-Δ)^{-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 $(-Δ)^s$ under minimal regularity assumptions, which is also of independent interest.