According to wikipedia, the Green's function of $\partial_x^2 + \partial_y^2$ is $\frac{1}{2\pi}\log{\sqrt{x^2+y^2}}$.
But how to prove this ?
For $x \neq 0$ or $y \neq 0$, $$\frac{d^2}{dx^2}\log(x^2 + y^2) = -\frac{2 (x^2 - y^2)}{(x^2 + y^2)^2}$$ Thus, $$\left(\partial_x^2 + \partial_y^2 \right)\frac{1}{2\pi}\log{\sqrt{x^2+y^2}} = 0$$
However, I cannot calculate the derivative for $x=0$ and $y=0$ because the numerator of $-\frac{2 (x^2 - y^2)}{(x^2 + y^2)^2}$ will expand.
So, my question is how to calculate $\left(\partial_x^2 + \partial_y^2 \right)\frac{1}{2\pi}\log{\sqrt{x^2+y^2}} $ for $x=0$ and $y=0$ to prove this function is Green's function of $\partial_x^2+\partial_y^2$.
====================Postscript=====================
I add this part according to the helpful comment of Cameron Williams.
In polar coordinates, $$\partial_x^2 + \partial_y^2 $$ $$= \partial_r^2 + \frac{1}{r}\partial_r+\frac{1}{r^2}\partial_\theta^2$$
and
$$\frac{1}{2\pi}\log(\sqrt{x^2+y^2})$$ $$=\frac{1}{\pi}\log{r}$$
Thus,
$$\left(\frac{d^2}{dx^2}+\frac{d^2}{dy^2}\right)\log(x^2 + y^2)$$ $$=\left(\partial_r^2 + \frac{1}{r}\partial_r\right)\frac{1}{\pi}\log{r}$$
Thing looks simpler. However, I still don't understand why $\delta(x)\delta(y)$ will appear for $x=0$ and $y=0$.
The Dirac delta is a distribution characterized by the property that $\newcommand{\rx}{x}$
$$ \int f(\rx) \delta(\rx) \, \mathrm{d}\rx = f(0) $$
holds for all nice test functions $f$. (Or more precisely, compactly supported smooth functions $f$.) And this is what we have to check. Now let $G$ be the 2D Green function. Then for each test function $f$, the distributional Laplacian $\Delta G$ (where of course $\Delta = \partial_1^2 + \partial_2^2$ is the Laplacian) satisfies
\begin{align*} \int f(\rx) \Delta G(\rx) \, \mathrm{d}\rx = \int G(\rx)\Delta f(\rx) \, \mathrm{d}\rx = \lim_{\epsilon \to 0^+} \int_{\|\rx\| > \epsilon} G(\rx)\Delta f(\rx) \, \mathrm{d}\rx. \tag{*} \end{align*}
Since $G(\rx)\Delta f(\rx)$ is smooth and $\Delta G(\rx) = 0$ on the region $U_{\epsilon} = \{\rx\in\mathbb{R}^2:\|\rx\|>\epsilon\}$, applying Green's second identity gives
\begin{align*} \int_{U_{\epsilon}} G(\rx)\Delta f(\rx) \, \mathrm{d}\rx &= \int_{U_{\epsilon}} \left( G(\rx)\Delta f(\rx) - f(\rx)\Delta G(\rx) \right) \, \mathrm{d}\rx \\ &= \int_{\partial U_{\epsilon}} \left( G(\rx)\nabla f(\rx) - f(\rx)\nabla G(\rx) \right) \cdot \mathrm{n}(\rx) \, \sigma(\mathrm{d}x), \end{align*}
where $\sigma$ is the surface measure of $\partial U_{\epsilon}$ and $\mathrm{n}(\rx)$ is the outward unit normal to $\partial U_{\epsilon}$ at $\rx$. But since
$$\nabla G(\rx) = \frac{\rx}{2\pi\|\rx\|^2} $$
and $\partial U_{\epsilon}$ is the inward-oriented circle of radius $\epsilon$ centered at $0$, we have $\mathrm{n}(\rx) = -\rx/\epsilon$ and hence
\begin{align*} &\int_{\partial U_{\epsilon}} \left( G(\rx)\nabla f(\rx) - f(\rx)\nabla G(\rx) \right) \cdot \mathrm{n}(\rx) \, \sigma(\mathrm{d}\rx) \\ &\quad = \int_{\partial U_{\epsilon}} \left( \frac{\log \epsilon}{2\pi}\nabla f(\rx) - f(\rx) \frac{\rx}{2\pi \epsilon^2} \right) \cdot \left( - \frac{\rx}{\epsilon} \right) \, \sigma(\mathrm{d}\rx) \\ &\qquad = \frac{1}{2\pi\epsilon} \int_{\partial U_{\epsilon}} f(\rx) \, \sigma(\mathrm{d}\rx) - \frac{\log \epsilon}{2\pi} \int_{\partial U_{\epsilon}} \nabla f(\rx)\cdot \frac{\rx}{\epsilon} \, \sigma(\mathrm{d}\rx). \end{align*}
Since the first integral in the last line is the average of $f$ over $\partial U_{\epsilon}$, by continuity of $f$,
$$ \frac{1}{2\pi\epsilon} \int_{\partial U_{\epsilon}} f(\rx) \, \sigma(\mathrm{d}\rx) \xrightarrow[\epsilon \to 0^+]{} f(0). $$
On the other hand,
\begin{align*} \left| \frac{\log \epsilon}{2\pi} \int_{\partial U_{\epsilon}} \nabla f(\rx)\cdot \frac{\rx}{\epsilon} \, \sigma(\mathrm{d}\rx) \right| &\leq \frac{\log (1/\epsilon)}{2\pi} \int_{\partial U_{\epsilon}} \|\nabla f\|_{\infty} \, \sigma(\mathrm{d}\rx) \\ &= (\epsilon \log (1/\epsilon)) \|\nabla f\|_{\infty}, \end{align*}
which converges to $0$ as $\epsilon \to 0^+$. Therefore the limit in $\text{(*)}$ is computed as
$$ \lim_{\epsilon \to 0^+} \int_{\|\rx\| > \epsilon} G(\rx)\Delta f(\rx) \, \mathrm{d}\rx = f(0), $$
which tells that the distribution $\Delta G$ is exactly the Dirac delta.