The Green's function for the Laplacian in the plane is $$G(\mathbf{x},\mathbf{x_0}) = \frac{1}{4\pi}\log(||\mathbf{x}-\mathbf{x_0}||^2)$$
That is to say, $G$ satisfies $$\triangle_{\mathbf{x}}G = \delta(\mathbf{x}-\mathbf{x_0})$$
where $\triangle$ is the Laplacian operator, and $\delta$ is the Dirac function.
In an effort to convince myself of this, I computed by hand the Laplacian of $G$, and obtained $$\triangle_{\mathbf{x}}G = \frac{-1}{2\pi}\frac{2r^2-2r^2}{r^4}$$ where $r= ||\mathbf{x}-\mathbf{x_0}||^2$.
Two questions:
Is this really equal to the $\delta$ function? Certainly it's equal to $0$ for $r \neq 0$, but for $r=0$ can we really say that it blows up? I suppose since $r > 0$, the value tends to $\infty$.
The more important question: Why does the initial factor of $\frac{1}{4\pi}$ even matter? We could have started out with any arbitrary constant and obtained the same property for the Laplacian.
We can redefine the problem to circumvent use of the Dirac Delta. So, rather than describing the problem by the expression $\nabla^2 G(\vec \rho|\vec \rho')=\delta(\vec \rho-\vec \rho')$ we describe the function $G(\vec \rho|\vec \rho')$ to satisfy the equations
$$\begin{align} \nabla^2 G(\vec \rho|\vec \rho')&=0\,\,\ \dots\vec \rho\ne \vec \rho'\tag 1\\\\ \lim_{\epsilon\to 0}\oint_{|\vec \rho-\vec \rho'|=\epsilon}\nabla' G(\vec \rho|\vec \rho')\cdot \hat n'\,d\ell'&=1\tag2 \end{align}$$
As determined in the OP, $(1)$ is satisfied for $G(\vec \rho|\vec \rho')=\frac{1}{2\pi }\log(|\vec \rho-\vec \rho'|)$.
To see that this form of $G$ satisfies $(2)$ we see that on $|\vec \rho-\vec \rho'|=\epsilon$, we have $\color{blue}{\nabla' G(\vec \rho|\vec \rho')=\frac{1}{2\pi}\frac{\vec \rho'-\vec \rho}{|\vec \rho-\vec \rho'|^2}}$, $\color{red}{\hat n'=\frac{\vec \rho'-\vec \rho}{|\vec \rho-\vec \rho'|}}$, and $\color{green}{d\ell'=|\vec \rho-\vec \rho'|\,d\phi'}$, where $\phi'\in [0,2\pi]$. Putting everything together reveals
$$\oint_{|\vec \rho-\vec \rho'|=R} \color{blue}{\nabla'G(\vec \rho|\vec \rho')}\cdot \color{red}{\hat n'}\,\color{green}{d\ell'}=\int_0^{2\pi} \color{blue}{\left(\frac{1}{2\pi}\frac{\vec \rho'-\vec \rho}{|\vec \rho-\vec \rho'|^2}\right)}\cdot \color{red}{\left(\frac{\vec \rho'-\vec \rho}{|\vec \rho-\vec \rho'|}\right)}\,\color{green}{\left(|\vec \rho-\vec \rho'|\,d\phi'\right)}=1$$
as was to be shown! In fact, note that we didn't need to let $\epsilon\to 0$. We only require that $\epsilon$ is small enough so that the disk $|\vec \rho-\vec \rho'|=\epsilon$ is contained inside the open domain of interest.
We show now the motivation for describing the Green (or Green's) function to the problem $\nabla^2 \phi(\rho)=f(\rho)$, where $f$ is a smooth function with compact support. We write $\phi(\vec \rho)$ as the superposition integral
$$\phi(\vec \rho)=\int_S G(\vec\rho|\vec\rho')f(\vec \rho')\,dS'\tag 3$$
We have $\nabla^2 G(\vec \rho|\vec\rho')=0$ for $\vec\rho\ne\vec\rho'$. In addition, we require $\phi$ to satisfy $\nabla^2 \phi(\vec \rho)=f(\vec \rho)$.
Accordingly, we apply the Laplacian to $(3)$ and write
$$\begin{align} \nabla^2 \phi(\vec \rho)&=\nabla \cdot \int_S \nabla G(\vec\rho|\vec\rho')f(\vec \rho')\,dS'\\\\ &=-\nabla \cdot \int_S \nabla' G(\vec\rho|\vec\rho')f(\vec \rho')\,dS'\\\\ &=-\nabla \cdot \int_S \left(\nabla'\left(G(\vec\rho|\vec\rho')f(\vec \rho')\right)-G(\vec\rho|\vec\rho')\nabla'f(\vec \rho')\right)\,dS'\\\\ &=-\nabla \cdot \underbrace{\oint_{\partial S} G(\vec\rho|\vec\rho')f(\vec \rho') \,\hat n'\,d\ell'}_{=0\,\,\text{since}\,f\,\text{is of compact support}}+\nabla \cdot \int_S G(\vec\rho|\vec\rho')\nabla'f(\vec \rho')\,dS'\\\\ &=-\int_S \nabla' G(\vec\rho|\vec\rho')\cdot \nabla'f(\vec \rho')\,dS'\\\\ &=-\int_S \nabla' \cdot \left(\nabla' G(\vec\rho|\vec\rho')f(\vec \rho')-f(\vec \rho')\underbrace{\nabla'^2 G(\vec \rho|\vec \rho')}_{=0}\right)\,dS'\\\\ &=-\underbrace{\oint_{\partial S}\hat n'\cdot\nabla' G(\vec\rho|\vec\rho')f(\vec \rho')\,d\ell'}_{=0,\,\text{since}\,f\,\text{is of compact support}}+\oint_{|\vec \rho-\vec \rho'|=\epsilon}\hat n'\cdot \nabla' G(\vec\rho|\vec\rho')f(\vec \rho')\,d\ell'\\\\ &=\underbrace{\oint_{|\vec \rho-\vec \rho'|=\epsilon}\hat n'\cdot \nabla' G(\vec\rho|\vec\rho')(f(\vec \rho')-f(\vec \rho))\,d\ell'}_{\to 0\,\,\text{as}\,\,\epsilon \to 0}+f(\vec \rho)\oint_{|\vec \rho-\vec \rho'|=\epsilon}\nabla' G(\vec\rho|\vec\rho')\cdot \hat n'\,d\ell' \end{align}$$
Therefore, we see that $\nabla^2 \phi(\vec \rho)=f(\vec \rho)$ if $\nabla^2 G(\vec\rho|\vec\rho')=0$ for $\vec \rho\ne \vec \rho'$ and
$$\lim_{\epsilon \to 0}\oint_{|\vec \rho-\vec \rho'|=\epsilon}\nabla' G(\vec\rho|\vec\rho')\cdot \hat n'\,d\ell'=1$$
This is precisely the description given by $(1)$ and $(2)$ and establishes a way forward without using distributions.