Let $V \ge 0$ be a continuous function $\mathbb R \to \mathbb R$ such that $V(x) \to \infty$ as $x \to \pm \infty$. Let $M > 0$ and define $T_M u = -\Delta u + M V u$ on $\mathcal D = \{u \in H^2(\mathbb R) : V u \in L^2(\mathbb R)\}$, understood in $L^2(\mathbb R)$. We know that $T_M$ will be self-adjoint (at least after taking closure in $L^2(\mathbb R)$) and have a non-negative purely discrete spectrum for each $M > 0$, let's write this as $\{\lambda_1^{(M)}, \lambda_2^{(M)}, \ldots\}$. I want to know how $\sigma(T_M)$ scales varying $M$. Specifically, I completely intuitively suspect that $\lambda_j^{(M)} \to \infty$ as $M \to \infty$ for $j \ge 2$. It's plausible to think that we could have $\lambda_1^{(M)} = 0$ for each $M > 0$, but if this is not the case I think we should have $\lambda_1^{(M)} \to \infty$ as $M \to \infty$ as well. This guess holds for $V(x) = x^2$, where we have $\sigma(T_M) = \sqrt M \sigma(T_1)$.
I was not too sure where to start with this. We can write out: $$\lambda_1^{(M)} = \inf_{\phi \in \mathcal D, \, \lVert \phi \rVert = 1} \langle (-\Delta + M V)\phi, \phi\rangle$$ to see that $\lambda_1^{(M)} \ge \lambda_1^{(N)}$ if $M \ge N$, but I am not too sure how to do something a bit more sophisticated. The corresponding expressions for $\lambda_j^{(M)}$ for $j \ge 2$ seem less useful since they concern eigenspaces of $T_M$, so you will have dependence in both the inner product and range of $\phi$.
Are there any results that would be directly applicable here, or is an approach specialised to the specific $V$ needed?
The general behaviour depends on the way your potential “looks” close to its zeros. For instance, in the case where your zeros are all non-degenerate critical points this paper proves that the eigenvalues will scale like $\sqrt{M}$ (relating this back to your example of $V=x^2$, the idea is that one can replace $V$ by its Taylor approximation close to its zeros, and use the explicit scaling available in that case). In general, for critical points of order m, the proof can be generalized in a straightforward fashion to yield a scaling of $M^\frac{2}{2+m}$. For a non-quantitative result, it is easy to see that as long as the zero set of your potential is so that the only $H^1$ function whose support is contained in it is the zero function, then the spectrum must tend to infinity as $M$ does. Indeed, suppose not. Then one can find a sequence $u_M$ of $H^1$ functions of norm 1 so that $$ ||\nabla u_M||^2+M\int u_M^2Vdx\leq C $$ for some C. Rellich+the fact that the potential diverges for large $|x|$ implies that we can find a a converging subsequence in $L^2$ converging to some $u\in H^1$. By fatou, it follows that $uV=0$ a.e. Since $u$ is in $H^1$ with unit $L^2$ norm, this is a contradiction.