The Laplace operator in cylindrical coordinates reads
\[\begin{eqnarray} \vec \nabla ^2= \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial }{\partial r}\right) +\frac{1}{r^2}\frac{\partial^2}{\partial\theta^2}+\frac{\partial^2}{\partial z^2} . \tag{1} \end{eqnarray}\] The Green’s function \(G\) is defined as the solution of the following differential equation: \[\begin{eqnarray} \vec \nabla ^2 G(\vec r)= \delta^3(\vec r) =\frac{1}{r}\delta(r)\delta(\theta)\delta(z). \tag{2} \end{eqnarray}\] We will consider problems with no \(\theta\) and \(z\) dependence, and that will leave behind the \(r\) dependence. Let’s solve this for \(r>0\): \[\begin{eqnarray} \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial G}{\partial r}\right) =0 \implies r\frac{\partial G}{\partial r}=C \implies G(r)=C \ln r+ D, \tag{3} \end{eqnarray}\] where \(C\) and \(D\) are constants. The term \(D\) is not that interesting since it will drop up on being acted on by \(\frac{d}{dr}\), and therefore we can set it to \(0\). However, we do need to figure out what \(C\) is. Although, \(\ln(r)\) diverges at \(r=0\), the overall expression with the derivatives vanishes. This will require a lot of care to handle the singularity properly. To this end, let’s protect the function by introducing a parameter \(\epsilon\), which we will set to zero when all is said and done: \[\begin{eqnarray} G_\epsilon(r)\equiv C \ln (r+\epsilon), \tag{4} \end{eqnarray}\] Inserting this back into Eq. (3), we get \[\begin{eqnarray} \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial G_\epsilon}{\partial r}\right) = \frac{C}{r}\frac{\partial}{\partial r}\left(\frac{r}{r+\epsilon}\right) = \frac{C}{r}\frac{\partial}{\partial r}\left(1-\frac{\epsilon}{r+\epsilon}\right) = \frac{\epsilon C}{r(r+\epsilon)^2}=\frac{1}{r}\delta(r)\delta(\theta)\delta(z) \tag{5}. \end{eqnarray}\] Integrating this over whole space:
\[\begin{eqnarray} \int_0^\infty dr \int_0^{2\pi} d\theta \int_{-h/2}^{h/2} dz\, r \frac{\epsilon C}{r(r+\epsilon)^2}= 2\pi h C=\int_0^\infty dr \int_0^{2\pi} d\theta \int_{-h/2}^{h/2} dz\frac{r}{r}\delta(r)\delta(\theta)\delta(z) =1 \implies C = \frac{1}{2\pi h}\tag{6}, \end{eqnarray}\] where \(h\) is an arbitrary length along the \(z\) axis. Note that we need \(h\) so that the units make sense. In pure mathematical expression, one wouldn’t care about it. However, imagine that \(r\) has the unit of length \(L\). \(\nabla^2\) has the unit of \(L^{-2}\). \(\delta^3(\vec r)\) has the unit of \(L^{-3}\). To match the units, \(G\) has the have the unit as \(L^{-1}\), which comes from the coefficient \(C\). \(h\) can be set to \(1\), a unitless value, and the unit can be absorbed into \(G\) to get:
\[\begin{eqnarray} G(r)=\frac{1}{2\pi } \ln(r)\tag{7}. \end{eqnarray}\]
The \(\epsilon\) business is a tricky one, and we could have avoided it if we invoked the Gauss theorem by integrating Eq. (2): \[\begin{eqnarray} \int_\mathcal{V} dV \vec \nabla \cdot \left(\vec \nabla G(\vec r)\right)= \int_\mathcal{\partial V} d\vec S \cdot \vec \nabla G(\vec r)= \int_0^{2\pi} d\theta \int_{-h/2}^{h/2} dz r \frac{d}{dr}(C ln(r))=2\pi h C = \int_\mathcal{V} dV \delta^3(\vec r)=1, \tag{8} \end{eqnarray}\] which would return the same value of \(C\).