There is a more general result. Fix an even Schwartz function $\newcommand{\bR}{\mathbb{R}}$ $w:\bR\to[0,\infty)$.
Let $\Delta$ be the Laplacian of the compact connected Riemann manifold $(M,g)$, $\dim M=m$. Its eigenvalues are
$$0=\lambda_0< \lambda_1\leq \lambda_2\leq \cdots$$
where each eigenvalue appears as many times as its multiplicity. Fix an orthonormal eigen-basis $(\Psi_k)_{k\geq 0}$ of $L^2(M,g)$,
$$\Delta\Psi_k=\lambda_k\Psi_k. $$
For each $\newcommand{\ve}{\varepsilon}$ $\ve >0$ define $\Xi_\ve: M\to L^2(M,g)$ by setting
$$\Xi_\ve(p)= \left(\frac{\ve^{m+2}}{d_m}\right)^{\frac{1}{2}}\sum_{k\geq 0}w\bigl(\,\ve \sqrt{\lambda_k}\,\bigr)^{\frac{1}{2}}\Psi_k(p)\Psi_k, $$
where
$$d_m:=\frac{2\pi^{\frac{m}{2}}}{m \Gamma(\frac{m}{2})}\int_0^\infty w(r) r^{m+1} dr.$$
Then for $\ve>0$ sufficiently small the map $\Xi_\ve$ is an embedding. Moreover, as $\ve\to 0$ the induced metric converges to the original metric. No assumption on the metric $g$ is required. Note that when $w$ is compactly supported the above sum consists of finitely many terms so $\Xi_\ve$ is actually an embedding into a finite dimensional space.
The result of Berard-Besson-Gallot corresponds to $w(r)=e^{-r^2}$. For details see this paper.