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$ 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.

The result of Berard-Besson-Gallot corresponds to $w(r)=e^{-r^2}$. For details [see this paper][1].


  [1]: https://arxiv.org/abs/1209.0639