Finding singular “solutions” to the Dirichlet problem for Schrödinger operators that do not admit smooth solutions

Suppose that $\Omega\subset \mathbb{R}^n$ is a smooth open region and that $V:\Omega\to \mathbb{R}^+$ be a positive smooth function.
Then we have a family of operators $$L_\epsilon =-\Delta -\epsilon V.$$ It is straightforward to see that there is an $\epsilon_{crit}(V,\Omega)>0$ so that for $\epsilon<\epsilon_{crit}$ we have $$\lambda_1(L_\epsilon)>0$$ and $$\lambda_1(L_{\epsilon_{crit}})=0.$$

Consider now the following Dirichlet problem $$L_\epsilon u= 0 \mbox{ where } u|_{\partial \Omega}=1.$$

When $\epsilon<\epsilon_{crit}$ there is a unique solution smooth solution whereas for $\epsilon=\epsilon_{crit}$ there is no smooth solution (which follows from the Hopf boundary maximum principle and Green's identity).

My questions are:

Is there a natural way to make the limit $$\lim_{\epsilon \nearrow \epsilon_{crit}}u_\epsilon:= u_{\epsilon_{crit}}$$ make sense -- in particular as some sort of singular solution to the Dirichlet problem? Clearly, both the convergence and $u_{\epsilon_{crit}}$ would need to be singular in an appropriate sense.

If this is the case can one "localize" the singularities for instance restrict where they are allowed to occur in terms of $V$ and $\Omega$?

Has this been considered in the literature at all?

$$-\Delta u =\lambda V u \textrm{ in }\Omega, \quad u=0 \textrm{ on } \partial\Omega.$$ Now, $N(u,v) = \int_\Omega V u v\, dx$ is a scalar product as $V$ is positive on $\Omega$ (equivalent to the usual one). The value $\lambda=\epsilon_{\textrm{crit}}$ is the first eigenvalue of the (inverse) laplacian for this scalar product , and as (inverse of the) Laplacian is compact and positive, you know that it is a simple and isolated eigenvalue, with an eigenvector $\psi_{\textrm{crit}}$, normalised such that $$\max_{\Omega} \psi_{\textrm{crit}}=1$$ for example, which does not change sign on $\Omega$. Consider now the problem you asked, that is $$-\Delta u =\epsilon_{\textrm{crit}} V u \textrm{ in }\Omega, \quad u=1 \textrm{ on } \partial\Omega.$$ Change unknowns to $w=u-1$, to obtain $$-\Delta w =\epsilon_{\textrm{crit}} V w + \epsilon_{\textrm{crit}} V \textrm{ in }\Omega, \quad w=0 \textrm{ on } \partial\Omega.$$ If integrate by parts against $\psi_{\textrm{crit}}$, you obtain $$0 = \int_\Omega \psi_{\textrm{crit}}\epsilon_{\textrm{crit}} V dx$$ But that is a contradiction, since the integrand is positive. Therefore there is no solution to that problem in $H^1(\Omega)$.
If you look at the family of solutions $u_{\epsilon}$ with $\epsilon<\epsilon_{\textrm{crit}}$, and test against $\psi_{\textrm{crit}}$ you obtain $$\int_\Omega V u_{\epsilon} \psi_{\textrm{crit}} = \frac{\epsilon_{\textrm{crit}}}{\epsilon_{\textrm{crit}}-\epsilon} \int_\Omega V \psi_{\textrm{crit}}$$ Since everyone is positive on any $\bar\omega\subset\Omega$, you have $$\|u_{\epsilon}\|_{L^1(\omega)} > \frac{C(\omega)}{\epsilon_{\textrm{crit}}-\epsilon}$$ for some constant depending on $\omega$, $V$ and $\psi_{\textrm{crit}}$ but independent of $\epsilon$. Therefore no hope for an $L^1$ solution as a limit either. So you would have to define a notion of solution which is not $L^1$ on any measurable set with closure included in $\Omega$...