In the question 2D Schrodinger equation in polar coordinates, @Qmechanic has shown that $\mathrm{Y}_0(r)\sim \ln r$ should be omitted since it's energy expectation value will diverge:
$$ \def\dd{\mathrm{d}} \def\bm{\boldsymbol} \int_0^{a} \psi^*\nabla^2\psi\, \dd^2 \bm{x} =\int_0^{a} -\frac{\ln r}{r^2} r \, \dd r \, \dd\theta \to\infty $$
However, if we calculate the Neumann function instead of logarithm, it is actually not divergency:
$$ \def\Y{\mathrm{Y}} \nabla^2_{r,\theta}\Y_0(r)=-\Y_0(r) \implies\int_0^{a} \psi^*\nabla^2\psi\, \dd^2 \bm{x} =\int_0^{a} -\Y_0^2(r) \, r \, \dd r \, \dd\theta <\infty $$
and the function itself is normalizable similarly.
So should we omit this solution? How can we explain the infinity at the origin?