Solution (source code)

= Solution

Use the <positive Laplace-Beltrami operator> $\Delta=-\operatorname{div}\operatorname{grad}$. A <Riemannian submersion> is a surjective smooth <submersion> $\pi:(M,g)\to(N,h)$ for which, at every $p$, the restriction of $d\pi_p$ to $H_p=(\ker d\pi_p)^\perp$ is a <linear isometry> onto $T_{\pi(p)}N$. The spaces $V_p=\ker d\pi_p$ and $H_p$ are its vertical and horizontal spaces. Its <fibres> are <totally geodesic submanifolds> precisely when $\nabla_UV$ is vertical for vertical <vector fields> $U,V$: their <second fundamental form> vanishes. Equivalently, a <geodesic> initially tangent to a fibre remains in that fibre while defined.

For a smooth $f:N\to\mathbb R$, its <basic function> $F=f\circ\pi$ is constant along each fibre. The <Riemannian gradient> of $F$ is the <horizontal lift of a vector field through a submersion> of $\operatorname{grad}_Nf$, since
$$
g(\operatorname{grad}_MF,X)=dF(X)=df(d\pi X)=h(\operatorname{grad}_Nf,d\pi X)
$$
for horizontal $X$, and $dF(U)=0$ for vertical $U$. In particular no derivative of $f$ in a vertical direction occurs.

Here is the needed connection fact, which also follows directly from the <Koszul formula>: for horizontal lifts $X,Y$ of <vector fields> $\bar X,\bar Y$ on $N$, the horizontal component of $\nabla_XY$ projects to $\nabla^N_{\bar X}\bar Y$. To see this, pair the <Koszul formula> with a third horizontal lift $Z$. The horizontal <inner products> are pulled back from $N$, and the horizontal components of their <Lie brackets of vector fields> project to the brackets on $N$. Thus all six terms are the pullbacks of the corresponding terms on $N$.

Choose an adapted <Riemannian orthonormal frame> $X_1,\ldots,X_n,U_1,\ldots,U_r$, with the $X_a$ horizontal lifts. Using $\operatorname{Hess}F(A,A)=A(AF)-(\nabla_AA)F$, the connection fact gives
$$
\operatorname{Hess}_MF(X_a,X_a)=(\operatorname{Hess}_Nf)(\bar X_a,\bar X_a)\circ\pi.
$$
For a vertical $U_b$, both $U_bF=0$ and $(\nabla_{U_b}U_b)F=0$, the latter because the fibres are <totally geodesic submanifolds>. Taking the negative <metric trace> of the <Riemannian Hessian> therefore proves the <basic-function Laplacian identity>
$$
\boxed{\Delta_M(f\circ\pi)=(\Delta_N f)\circ\pi.}
$$
This identity is local and does not require <compactness>. Vanishing <mean curvature> of the <fibres> would already suffice; total geodesicity makes each vertical summand vanish separately.

For the discrete <eigenspace> assertion, assume the two <Riemannian manifolds> are <closed manifolds>. Without a discrete spectral realization, an unrestricted noncompact version need not have an <eigenbasis>. The projections of the <Riemannian product> $M\times N$ are <Riemannian submersions> with <totally geodesic submanifolds> as fibres. Its <Levi-Civita connection> splits into the two factor connections. Consequently its <positive Laplace-Beltrami operator> is
$$
\Delta_{M\times N}=\Delta_M\otimes I+I\otimes\Delta_N,
\qquad
\Delta_{M\times N}(u(x)v(y))=(\Delta_Mu)(x)v(y)+u(x)(\Delta_Nv)(y).
$$
The cross term in the <product rule for the positive Laplace-Beltrami operator> is zero because the two factor <Riemannian gradients> are <orthogonal>.

We use the standard compact elliptic <compact elliptic spectral theorem>: the <positive Laplace-Beltrami operator> on a <closed manifold> is <self-adjoint>, has <compact resolvent>, and has a complete <orthonormal eigenbasis> of smooth <eigenfunctions>, with finite-dimensional <eigenspaces> and <eigenvalues> tending to infinity. Let $\Delta_Mu_i=\mu_i u_i$ and $\Delta_Nv_j=\nu_jv_j$. <Fubini's theorem> and completeness on each factor show that $u_i(x)v_j(y)$ form a complete <orthonormal basis> of $L^2(M\times N)$. For example, a function orthogonal to all these products has, for each $i$, zero $u_i$-coefficient as an $L^2(N)$ function, hence is zero.

The displayed operator identity makes $u_i v_j$ an <eigenfunction> with <eigenvalue> $\mu_i+\nu_j$. Conversely, if $\Delta w=\lambda w$, self-adjointness of the <positive Laplace-Beltrami operator> gives
$$
(\mu_i+\nu_j-\lambda)\langle w,u_i v_j\rangle=0.
$$
All other coefficients vanish. Only finitely many pairs can have $\mu_i+\nu_j=\lambda$, since both <spectra> are nonnegative and have finitely many <eigenvalues> below any fixed bound. Thus the <product Laplacian eigenspace decomposition> is
$$
\boxed{E_\lambda(M\times N)=\bigoplus_{\mu+\nu=\lambda} E_\mu(M)\otimes E_\nu(N).}
$$
The <tensor product> summands are mutually <orthogonal>; their elements are actual smooth <eigenfunctions>, so this is an equality of <eigenspaces>, not just a formal expansion.