How can the Bessel function of the second kind be in the radial eigenfunction? Let $0<a<b<\infty$. Consider the heat equations or wave equations on the annulus or the spherical layers $$\Omega:=\{x\in\mathbb{R}^d\mid a<\|x\|_2<b\},$$ $$\partial\Omega:=\{x\in\mathbb{R}^d\mid \|x\|_2=a \}\cup \{x\in\mathbb{R}^d\mid \|x\|_2=b \},$$ i.e., the annulus ($d=2$) or the spherical layers ($d=3$) with the Dirichlet boundary condition. 
More compactly, a Bessel function of the second kind in eigenfunctions of a Laplace operator with Dirichlet boundary conditions on a spherical domain.
We assume the initial condition is nice enough, such that no subtleties arise and we can use the separation of variables.
Question: How can the Bessel function of the second kind be in the radial eigenfunction?
Context: When we use the separation of variables, the problem boils down to solving two ODEs in the angular and the radial part. In the radial part, when we consider the disk or the ball ($a=0$), we impose the condition $u(0)<+\infty$ at the origin ($u$ is the solution), which eliminates the Bessel function of the second kind from the consideration for the eigenfunction in the radial part.
Suppose $a>0$. I have found, e.g., in http://arxiv.org/pdf/1206.1278.pdf, in the radial part people consider both the Bessel function of the first and second. But if you see the plot of the second kind, http://mathworld.wolfram.com/BesselFunctionoftheSecondKind.html, you see it is unbounded not only in the neighborhood of the origin but the "unboundedness-point" shifts rightward.
How can one consider non-trivial eigenfunctions with Dirichlet boundary condition, but with possibly unbounded functions?
Also, are there books that treats $d=3$ in detail? Feels like it would be something like "eigenvalues are taken large enough so that we can avoid the unboundedness part", but I would like to see e.g., how large, in relation to the unboundedness part etc.
 A: I think you have misunderstood the nature of the solutions to Bessel's equation: the differential equation
$$ y'' + \frac{1}{x}y' + \left(\lambda^2-\frac{n^2}{x^2}\right) y = 0 $$
(Bessel's equation, with eigenvalue $\lambda^2$, $\lambda>0$) has a regular singular point at $x=0$, an irregular singular point at $\infty$, and all other points are regular points. At a regular point, a differential equation has two linearly independent solutions, each analytic on a disk about the point. On the other hand, at a regular singular point, there are two linearly independent solutions which either look like $x^{\sigma}$ or $x^{\sigma}\log{x}$, analytic on a cut disk about the point.
The upshot of this is that, by analytic continuation, the solutions to Bessel's equation are always analytic on the the plane with a single branch cut, which is normally chosen to be the negative real axis. Therefore, any solution to Bessel's equation on the interval $[a,b]$, $b>a>0$, will be analytic, and you just have to find the $\lambda$ such that
$$ \begin{vmatrix} J_n(\lambda b) & J_n(\lambda a) \\ Y_n(\lambda b) & Y_n(\lambda a) \end{vmatrix} = 0 $$
for the eigenvalues. (Since the equation has a nonzero solution if and only if this determinant vanishes.)
