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 eigenvalues appears as many times as its multiplicity. Fix an orthonormal eixgen-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