Eigenspace of Graph Laplacian related to clusters? Given an undirected graph $G$ with adjacency matrix $A$ and Graph Laplacian $L=D-A$, where $D$ is the degree matrix.
I wonder whether the eigenspaces of $L$, let's call them $(V_i)_{i=1}^n$, contain information about the clusters on the graph. 
 A: Yes. This is exactly the principle behind spectral clustering. The basic idea is that your adjacency matrix (more generally, any affinity matrix) describes the local relation between points through $L$. The eigendecomposition of $L$, the $(\lambda_i, V_i)$ eigenpairs, describes a new vector space in which the nodes are embedded, such that clustering the points there is very effective.
More specifically, however, it can be shown that solving 
$Lv = \lambda D v$ for the eigenvector $u$ corresponding to the second-smallest eigenvalue is the solution to the normalized cut problem, partitioning the graph $G=(V,E)$ into two parts (clusters) $A$ and $B$ such that 
$$
\mathfrak{N}(A,B) = \frac{C(A,B)}{T(A)} + \frac{C(B,A)}{T(B)},
\;\;\;C(A,B)=\sum_{u\in A, v\in B} W(u,v),\;\;\; T(A) = \sum_{u\in A,t\in V} W(u,t)
$$
is minimal.
Subsequent eigenvectors optimally subpartition these clusters.
The usefulness of $L$ for clustering can be interpreted in many ways. One is to consider the graph as a system of masses connected by springs: when shaken, the resulting oscillations (described by the Laplacian harmonics) "group" the closely connected masses that move together.
Another interpretation is to consider the normalized affinity (or adjacency) matrix $P$, so that it represents the probability that a random Markov walker moves from one node to another. The matrix powers $P^t$ represent node transition probabilities after $t$ steps. The Random Walk Laplacian is then $L=I-P = I - D^{-1} A$. 
Notice that the eigenpairs of $L$ are then very closely related to those of $P$. Computing $P^t$ can be easily done using the eigendecomposition of $L$. The interpretation is that we are grouping nodes together that a random walker easily reaches, while leaving out those to which it is difficult to diffuse.
(Diffusion maps are closely related to this approach as well.)
A great reference for all this is Shi and Malik, Normalized Cuts and Image Segmentation.
In the continuous case, when $L$ is viewed as an approximation of the Laplace-Beltrami operator $\Delta_M$ and the graph is a triangluation (mesh) of samples from a manifold $M$, its eigenfunctions (the manifold harmonics) are known to cluster the graph nodes just as well, which has applications to mesh segmentation in computer graphics (see e.g., Rustamov, Laplace-Beltrami Eigenfunctions for Deformation Invariant Shape Representation). The relation to heat diffusion still holds here as well, since the heat equation is $\Delta_M f = \partial_t f$ with solution (heat kernel) given by $K_t(x,y)=\sum_i \phi_i(x)\phi_i(y) \exp(-\lambda_i t)$, where $(\phi_i,\lambda_i)$ are the eigenpairs of $\Delta_M$.
