First, let me apologize if the premise of this question is incorrect - I have no formal training in PDE's and am only fitfully self-taught. I can read through and mostly understand the proofs via separation of variables for a given coordinate system, and have a rough understanding of generalized Fourier series.
Question: Intuitive explanation why solution space for Helmholtz equation in $\mathbb{R}^N$ is finite-dimensional for $N=2$ but infinite-dimensional for $N \ge 3$? (I am interested in "infinite domain".)
Guess: Is it "heuristically correct" to say that the solution space "has a dimension for every piece of information required to specify any particular solution"???????
Since the 1D Helmholtz equation is a 2nd order ODE, only "two pieces of information" are required to specify a solution uniquely, so for any choice of $k \not=0$ the solution space is 2D.
But for the ND Helmholtz equation, any solution of the corresponding PDE, regardless of the value of $k \in \mathbb{R}$, requires a "boundary condition" to specify, which in particular requires specifying the value of the solution and/or one of its (partial) derivatives at infinitely many points, so the solution space needs to be infinite dimensional??????
In particular, the fact that any solution for a given $k$ can be written in terms of a generalized Fourier series of particular solutions we get using separation of variables reflects the fact that we have "infinite degrees of freedom", since being able to specify the values of infinitely many generalized Fourier coefficients requires the "infinite pieces of information" contained in a boundary condition???????
Clarification of my limited understanding in case premise of question is incorrect:
Here is my understanding of the situation:
For $N=1$, the Helmholtz equation for a given $k \in \mathbb{R}$ is just the ODE $f''(x) = k f(x)$. When $k > 0$, any solution is a linear combination of the exponentials $e^{\sqrt{k}x}$ and $e^{-\sqrt{k}x}$ (a 2D vector space, $2 < \infty$), and when $k < 0$, any solution is a linear combination of the sinusoids $\sin(\sqrt{-k}x)$ and $\cos(\sqrt{-k}x)$, which again is also a 2D vector space. Initial conditions allow one to choose the relevant solution in the vector space.
For $N=2$, $f_{xx} + f_{yy} = k f(x,y)$ there is apparently an infinite-dimensional vector space of solutions for the "infinite domain" (i.e. all of $\mathbb{R}^2$) problem (I am not actually sure about this though), so no obvious canonical "basis" for displaying them, but according to Wikipedia seemingly a "basis" in terms of products of sinusoids and Bessel functions will do, with every solution having "coordinates" with respect to that basis given in terms of a generalized Fourier series, and evaluating those coordinates requiring the specification of a boundary condition. Analogously for $N = 3$, $f_{xx} + f_{yy} + f_{zz} = k f(x,y,z)$, if one replaces "sinusoids" with "spherical harmonics" and "Bessel functions" with "spherical Bessel functions". I have no idea for $N \ge 4$, although I guess I conjecture the solution could involve a Fourier basis with functions that have "higher-dimensional spherical harmonics" and "higher-dimensional Bessel functions" as factors.
To be honest, as the above paragraph probably demonstrates, I don't actually understand the PDE case very well. In particular I am assuming that the coefficients in the generalized Fourier series depend on the boundary condition, but I am not actually sure whether that is the case or whether there is a completely different vector space of possible solutions for every admissible boundary condition.
I guess I am just struggling to generalize the intuition from the 1D case, where everything is very straightforward calculus 101, to higher dimensions - what generalizes, what doesn't, why, etc.
In particular I am confused why the analogous eigenfunction problem has finite-dimensional solution spaces for spherical harmonics and their higher-dimensional analogues. In particular it seems like it should be possible to write the spherical Laplacian as a differential operator on $\mathbb{R} ^ 2$. Yet the resulting differential operator would only have a finite-dimensional (linear combinations of the corresponding $2 \ell +1$ spherical harmonics) solution space, while the 2D Laplacian seemingly has an infinite dimensional solution space. Maybe it is because of the special geometry of the sphere meaning that we aren't really considering arbitrary solutions of that operator, but only those that satisfy the required periodicity conditions? And maybe those periodicity conditions can be thought of as "infinitely many boundary conditions", with the result that the solution space is now finite-dimensional and not infinite-dimensional? Perhaps the same argument also explaining why only certain values of $\ell$ will have solutions, instead of arbitrary $k \in \mathbb{R}$? (This is probably something I should ask as a separate question, but I just wanted clarification on the above first since I'm not sure whether that's the part I'm misunderstanding.)