This is, of course, a long story incorporating many strands but I will try to give a quick overview. Firstly, it is, as so often, convenient to skip to a more general framework. In your case, this would be that of an unbounded self-adjoint operator $T$ on Hilbert space (here that would be the Laplacian--more later).
One can associate with it a Frechet space $H^\infty(T)$ (the intersection of the domains of definitions of its powers) and a $DF$-space $H^{-\infty}(T)$ which are in duality. In the case of classical differential operators (the most frequent examples occur with the Laplacian and Schrödinger operators), the former is a space of test functions, the latter of distributions. This is a fairly direct consequence of the spectral theorem (in the form that any such operator can be represented as one of multiplication by a measurable function on an $L^2$-space).
If the spectrum of $T$ is discrete and consists of a sequence $(\lambda_n)$ of eigenvalues which are such that $|\lambda_n|$ is asymptotically like $n^\alpha$ for some positive $\alpha$, then the situation is particularly transparent. $H^\infty$ and $H^{-\infty}$ are a nuclear Frechet space and Silva space respectively. The eigenfunctions of $T$ form a basis for both spaces and they are identifiable with the sequence spaces $s$ and $s´$ of rapidly decreasing resp. slowly increasing sequences via coefficients.
The smoothing operators, i.e., continuous linear operators from $H^{-\infty}$ into
$H^\infty$ are then identified with the two variable version of $s$ in a standard way.
It is classical that the Laplacian satisfies these condition in many cases, e.g.,on a closed manifold, a compact manifold (with suitable boundary conditions--Dirichlet or Neumann), as does the Schrödinger operator for suitable potential functions. The requisite estimates on the eigenvalues go under the generic name of Weyl inequalities.
You cannot, however, expect such results in the non-compact case. Here the Schrödinger operator is, perhaps, more appropriate.