In Hajime Urakawa's monograph The Spectral Geometry of the Laplacian on page 41, we make an assumption that I can't quite justify on my own. The following is our setup:

Let $(M^n,g)$ be a closed connected Riemannian manifold, with $\Delta_g$ the Laplacian (a negative operator), and $0=\lambda_0<\lambda_1\le\lambda_2\le\cdots$ the eigenvalues of $-\Delta_g$ with corresponding eigenfunctions $u_0=\hbox{const},u_1,u_2,u_3,\ldots$. Then for $N\gg0$, $$\iota:M\ni x\mapsto(u_1(x),\ldots,u_N(x))\in\mathbb R^N$$ defines an embedding.

But why is this true? I know that $(u_i)_{i\ge0}$ forms an orthonormal basis for $L^2(M,\mu_g)$, but this doesn't seem to be quite enough to prove that we can have an embedding. Presumably we need to show that we can also approximate functions in $H^1(M,g)$?