Fundamental solution of Discrete Laplace in the plane

For nearest neighbor Laplacian, you can find a short self-contained proof of the formula $$u(x)=\frac1{2\pi}\log |x|+c+O\left(\frac1{|x|^2}\right)$$ at these lecture notes, pages 6-7. It is a simple exercise in Laplace's method to work out further coefficients of the asymptotic expansion, if necessary. The proof there is not the shortest possible; the easiest way would be to compare the double integral formula defining $G_0$ on page 6 with its counterpart for continuous Laplacian. But then you do not get the beautiful formula (4.4) along the way.

"Random Walk: a modern introduction" by Lawler and Limic does the same in much greater generality (not necessarily nearest neighbour), using local CLT with error bounds.

Historically, the first sources usually cited for the above formula are

Stohr, A., [50]: "Uber einige lineare partielle Differenzengleichungen mit konstanten Koefizienten," I, II, and Ill, Mathematische Nachrichten 3 (1950),

or even

W. H. McCrea and F. J. W. Whipple, Random paths in two and three dimensions, Proc. Roy. Soc. Edinburgh 60 (1940),

although the latter apparently proves a slightly weaker result (I have no access to it at the moment to check).

Formula (71) in Basic Properties Of Discrete Analytic Functions, R. J. Duffin, Duke Math. J. Volume 23, Number 2 (1956), 335-363 (warning: it is not the fist paper where it appears) gives asymptotics $$ 2\pi u(x,y)=\log 4r+\gamma+\cos(4\varphi)/6r^2+O(r^{-4}), $$ where $(x,y)=(r\cos \varphi,r\sin \varphi)$. Remainder term remains the same for discrete derivatives (they are the same as discrete differences like $u(x+1,y)-u(x,y)$, yes?), which in turn are genuine derivatives in intermediate points. Taking derivative of these functions allows to get much more than you ask about.