Let $V : \mathbb{R}^2 \to \mathbb{R}$ be compactly supported, bounded away from the origin, and obey $$ |V(x)| \lesssim r^{-\delta_0}, \qquad 0 < |x| \le 1, \qquad r : =|x|,$$ for some $0 < \delta_0 < 2$. Note that this encompasses potentials which are "Coulomb-like" near zero. In this situation, the quadratic form $$H^1(\mathbb{R}^n) \ni u \mapsto q_V(u,u) : = \int |u|^2 Vdx$$ is "relatively form bounded" with respect to the natural quadratic form associated to the free Laplacian: $$H^1(\mathbb{R}^2) \ni u \mapsto \| \nabla u \|^2_{L^2}.$$ More precisely, what this means that is that for each $a > 0$, there exists $b > 0$ so that $$ \| \sqrt{|V|} u\|^2_{L^2} \le a \| \nabla u \|^2_{L^2} + b \|u\|^2_{L^2}.$$ A proof of this estimate may be found in Faris, Pacific J. Math. 1967. A key tool in the proof is the Hausdorff-Young equality.
Combining the above bound with the KLMN Theorem implies that the Schrödinger operator $$H : = -\Delta + V : L^2(\mathbb{R}^2) \to L^2(\mathbb{R}^2)$$ is self-adjoint with respect to the domain $$ \mathcal{D} := \{u \in H^1(\mathbb{R}^2): Hu \in L^2(\mathbb{R}^2) \},$$ where the membership $Hu \in L^2(\mathbb{R}^2)$ is interpreted in the distributional sense (note $uV \in L^1_{\text{loc}}(\mathbb{R})$ for any $u \in H^1(\mathbb{R}^2)$, thanks to our above bounds on $|V(x)| \mathbf{1}_{0 < |x| \le 1}$ and $\| \sqrt{|V|} u\|_{L^2}$). See also Proposition 1.1, Chantelau, Lett. Math. Phys 1990.
I would like to know whether there exists some subspace $D_0 \subseteq C^\infty(\mathbb{R}^2 \setminus \{0\})$ which is dense in $\mathcal{D}$ with respect to the standard graph norm $u \mapsto \| \cdot \|_{\mathcal{D}} := (\|u\|^2_{L^2} + \|Hu\|^2_{L^2})^{1/2}$ on $\mathcal{D}$.
In general, to have the inclusion $D_0 \subseteq \mathcal{D}$, it seems some (at least implicit) condition has to be imposed on the behavior of $u \in D_0$ as $|x| \to 0$, in order to counteract the blow-up of $V$. I have tried the barehanded approach of simply mollifying an element $u \in \mathcal{D}$,
$$ u_\varepsilon(x) : = \varepsilon^{-2} \int \varphi(y/\varepsilon)u(x - y)dy \in C^\infty(\mathbb{R}^2), \qquad 0< \varepsilon \ll 1,$$
where $\varphi \in C^\infty_0(\mathbb{R}^2)$ is an approximate identity. However, I fail to see how $u_\varepsilon$ even belongs to $\mathcal{D}$ (multiplication by $V(x)$ doesn't "play nicely" with the mollification).
On the other hand, I doubt that $(H, \mathcal{D})$ is the closure of $(H, C_0^\infty(\mathbb{R}^2 \setminus \{0\})$, due to the usual issue with controlling derivatives of cutoffs near zero (the situation would not be as severe in higher dimensions).
If the result holds only for a smaller range of $\delta_0$, I would still find this interesting and would like to know what is the optimal range. In terms of more recent literature on this topic, I came across Oliveira-Verri, Ann. Physics 2009. This paper discusses self-adjoint extensions for Coulomb potentials in dimensions 1, 2, and 3, but doesn't seem to directly address my density question here.
Comments, hints, or pointers to the literature are greatly appreciated.
Edit: When $0 < \delta_0 < 1$, it then holds that $V \in L^2(\mathbb{R}^2)$, in which case $V$ is self-adjoint with respect to the standard Sobolev space $H^2(\mathbb{R}^2)$ (in which $C^\infty_0(\mathbb{R}^n)$ is dense). This follows from Theorem 8 from Nelson 1964, J. Mathematical Phys.
