Consequently, it follows that there exist second order, elliptic differential operator of the form $$ L u = a^{ij}\,\frac{\partial^2u}{\partial x^i\partial x^j} + b^i\,\frac{\partial u}{\partial x^i} + c\,u $$ (with $(a^{ij})$ positive definite) such that every local eigenfunction of $g$ satisfies an equation of the form $Lu = \phi(u)u$ where $\phi(u)$ is a smooth function (that could depend on $u$), and, moreover, $L$ is unique up to scalar multiplication and the addition of a $0$-th order term.

Consequently, it follows that there exist second order, elliptic differential operator of the form $$ L u = a^{ij}\,\frac{\partial^2u}{\partial x^i\,\partial x^j} + b^i\,\frac{\partial u}{\partial x^i} + c\,u $$ (with $(a^{ij})$ positive definite) such that every local eigenfunction of $g$ satisfies an equation of the form $Lu = \phi(u)u$ where $\phi(u)$ is a smooth function (that could depend on $u$), and, moreover, $L$ is unique up to scalar multiplication and the addition of a $0$-th order term.

Here is a sketch of an idea of how to show that, given the set $\mathcal{E}(g)\subset C^\infty(M)$ of all the eigenfunctions of the metric $g$ on a compact manifold $M$ will determine $g$ up to a constant multiple. (Note that I'm do not assume that the corresponding eigenvalues are given, which would be easier.) Clearly, this is the best one can hope for, since for any constant $c>0$, $\mathcal{E}(cg) = \mathcal{E}(g)$. In fact, I think it's very likely that knowing a sufficiently 'large' finite subset of $\mathcal{E}(g)$ should be sufficient, but that remains to be seen.

The basic idea is this: Let $J^k(M)\to M$ denote the vector bundle consisting of the $k$-jets of smooth functions on $M$. When $M$ has dimension $n$, the bundle $J^k$ has rank ${n+k}\choose n$ over $M$. Now, given $g$, it is easy to show that there is a closed quadratic cone bundle $Q(g)\subset J^3(M)$ of codimension $n$ such that the $3$-jets of all the local eigenfunctions of $g$ lie in $Q$. In fact, it is easy to show that the $3$-jets of local eigenfunctions fill out $Q(g)$.

Note that $Q(g)$ is not a linear subbundle of $J^3(M)$ precisely because we are not specifying the eigenvalues of the eigenfunctions. Of course, the $3$-jets of local eigenfunctions with eigenvalue $\lambda$, fill out a linear subbundle $Q(g,\lambda)\subset J^3$, but the union of the $Q(g,\lambda)$ is not a linear subbundle. However, it is easy to show that $Q_x(g)\subset J^3_x$ is the zero locus of an ideal generated by $n$ polynomials homogeneous of degree $2$ on the vector space $J^3_x$ for all $x\in M$. $Q_x(g)$ is not a smooth variety in $J^3_x$ but the singular locus is quite small.

I claim that the subbundle $Q(g)\subset J^3$ determines $g$ up to a constant factor (assuming that $M$ is connected). Here is why: Let $Q^0(g)\subset Q(g)$ denote the subset consisting of those $3$-jets in $Q(g)$ whose $0$-jet vanishes. The projection $\pi^2_3(Q^0(g)_x)$ of $Q^0(g)_x$ into $J^2_x$ has codimension $n{+}1$ in $J^2_x$, cut out by the linear equation that says that the $0$-jet vanishes and $n$ quadratic equations that turn out to all be multiples of a single linear equation on the space of $2$-jets whose $0$-jet vanishes, as is easy to verify in local coordinates.

Finally, what one would expect is that, when $M$ is compact, if we now look at the $3$-jets of the elements of $\mathcal{E}(g)$, i.e., the global eigenfunctions of $g$, that a sufficiently large subset will determine sufficiently many points in $J^3_x$ that they will determine $Q(g)_x$, which, after all is known to be cut out by $n$ homogeneous quadratic polynomials for each $x$. (Of course, the number of points needed could be large, even for $n=2$, but it will be finite.) Assuming such a 'density' result, $\mathcal{E}(g)$ will determine $Q(g)\subset J^3$, which, as we have seen, will determine $g$ up to a constant multiple.

Here is a sketch of an idea of how to show that the set $\mathcal{E}(g)\subset C^\infty(M)$ of all the eigenfunctions of the metric $g$ on a compact manifold $M$ determines $g$ up to a constant multiple. (Note that I do not assume that the corresponding eigenvalues are given, which would be easier.) Clearly, this is the best one can hope for, since for any constant $c>0$, $\mathcal{E}(cg) = \mathcal{E}(g)$. In fact, I think it's very likely that knowing a sufficiently 'large' finite subset of $\mathcal{E}(g)$ should be sufficient, but that remains to be seen.

The basic idea is this: Let $J^k(M)\to M$ denote the vector bundle consisting of the $k$-jets of smooth functions on $M$. When $M$ has dimension $n$, the bundle $J^k$ has rank ${n+k}\choose n$ over $M$. Given $g$, it is easy to show that there is a closed quadratic cone bundle $Q(g)\subset J^3(M)$ of codimension $n$ such that the $3$-jets of all the local eigenfunctions of $g$ lie in $Q(g)$. In fact, the $3$-jets of local eigenfunctions fill out $Q(g)$.

Note that $Q(g)$ is not a linear subbundle of $J^3(M)$ precisely because we are not specifying the eigenvalues of the eigenfunctions. Of course, the $3$-jets of local eigenfunctions with eigenvalue $\lambda$, fill out a linear subbundle $Q(g,\lambda)\subset J^3$, but the union of the $Q(g,\lambda)$ as $\lambda$ varies is not a linear subbundle. However, it is easy to show that, for any $x\in M$, the subset $Q_x(g)\subset J^3_x$ is the zero locus of an ideal generated by $n$ polynomials homogeneous of degree $2$ on the vector space $J^3_x$. ($Q_x(g)$ is not a smooth variety in $J^3_x$ but the singular locus is quite small.)

I claim that the subbundle $Q(g)\subset J^3$ determines $g$ up to a constant factor (assuming that $M$ is connected). Here is why: Let $Q^0(g)\subset Q(g)$ denote the subset consisting of those $3$-jets in $Q(g)$ whose $0$-jet vanishes. The projection $\pi^2_3(Q^0(g)_x)$ of $Q^0(g)_x$ into $J^2_x$ has codimension $n{+}1$ in $J^2_x$, cut out by the linear equation that says that the $0$-jet vanishes and the $n$ quadratic equations that then turn out to all be multiples of a single linear equation on the space of $2$-jets whose $0$-jet vanishes, as is easy to verify in local coordinates.

Finally, what one would expect is that, when $M$ is compact, if we now look at the $3$-jets of the elements of $\mathcal{E}(g)$, i.e., the global eigenfunctions of $g$, that a sufficiently large subset will determine sufficiently many points in each $J^3_x$ that they will determine $Q(g)_x$, which, after all, is known to be cut out by $n$ homogeneous quadratic polynomials for each $x$. (Of course, the number of elements of $\mathcal{E}(g)$ needed to do this at each point could be large, even for $n=2$, but it will be finite.) Assuming such a 'density' result, $\mathcal{E}(g)$ will determine $Q(g)\subset J^3$, which, as we have seen, will determine $g$ up to a constant multiple.