It can easily happen that the product of two eigenfunctions has an infinite eigenfunction expansion. Probably this is more _typical_ (in natural problems) than not. For example, the Laplace-Beltrami operator on compact Riemannian manifolds (or suitable values $(\Delta+c)^{-1}$ of its resolvent if one must have a _bounded_ operator) has discrete spectrum, but products of eigenfunctions are rarely eigenfunctions.

For example, on the sphere, (restrictions of) homogeneous harmonic polynomials are the eigenfunctions for the natural Laplacian. Products or squares of harmonic polynomials are rarely harmonic, but do have _finite_ expressions as sums of eigenfunctions. (From my viewpoint, this finiteness is predictable because the irreducibles of the rotation/orthogonal group are all finite-dimensional, because it is compact.)

Either from the a viewpoint of geometric analysis, or from a viewpoint of repn theory, in generality such re-expressions of products would rarely be finite, but should exist under mild hypotheses on the situation. Explicit examples where the decomposition/Fourier coefficients $\langle fg,h\rangle$ are explicitly expressible for triples of eigenfunctions are not so easy to manufacture, and the ones I know (having to do with automorphic forms, zonal spherical harmonics, etc.) I suspect are not in the direction of your interest.

It is also true that there seems to be fairly skimpy literature on such issues.

About your question 2: I don't have an easy example. About question 3: no, there are easy examples of the failure: on $[0,2\pi]$, with $f(x)=e^{ix}$ and $g(x)=e^{-ix}$ (or cosines and sines), and $E$ the space spanned by constants, $fg$ projects to $1$, while both products $f^2$ and $g^2$ project to $0$.

Edit: in response to the questioner's further comment/questions... First, the latter counter-example, with exponential functions, refers to $d^2/dx^2$.

Second, as @Matt Y. noted, someone like me is thinking of automorphic things, which may not be what everyone wants to see as examples... But, further, I think it is not so easy to make any such examples _explicit_. Even without the conjectured (and probable) typical-non-vanishing and so on, we can find definitive examples of infinite, discrete decompositions in the pseudo-Laplacians Y. Colin de Verdiere used c. 1981 to give another proof of meromorphic continuation of Eisenstein series: these are self-adjoint operators with compact resolvents (so discrete spectrum) on Hilbert spaces of automorphic functions which have, among their eigenfunctions _certain_ truncated Eisenstein series. The integral of three truncated Eisenstein series is not completely trivial to compute, but turns out (e.g., Zagier an others, circa 1978) to be a product/quotient of values of zeta at points that can be controlled directly... so can be made not to vanish. Thus, provably, an infinite-not-finite re-expression of product of two eigenfunctions (for a pseudo-Laplacian...) as sum of eigenfunctions. The "pseudo-Laplacian" is not wildly different from a/the genuine Laplacian...