It is common in QM, even in mathematical physics, to consider Hamiltonians with a Dirac delta in the potential:
$$
  \hat H = -\frac{\mathrm{d}^2}{\mathrm{d}x^2} + f(x) + \delta(x-a) \: .
$$
The most common interpretation of this is the self-adjoint operator $H$, such that:
$$
  H: \mathrm{D}(H) \subset L^2(\mathbb{R}) \to L^2(\mathbb{R}) \\[5pt]
  \mathrm{D}(H) = \big\{\;
    f \in W^{1,2}(\mathbb{R}) \cap W^{2,2}(\mathbb{R} \setminus \{0\})
    \;\; \big| \;\;
    f'(a+) - f'(a-) = f(a)
  \;\big\} \\[5pt]
  H\psi(x) = -\psi''(x) + f(x)\psi(x) \quad \text{on } \mathbb{R} \setminus \{0\}
$$
This is a very practical way to interpret the problem, however I feel like this approach sweeps the distributional aspect „under the rug“. I would love to understand the problem **in terms of distributions**.

My gut feeling is that something like this should be possible:
> Let $\mathcal{H} = L^2(\mathbb{R})$ and $\Phi = C^\infty_0(\mathbb{R})$ be the space of smooth functions with compact support, then $\Phi \subset \mathcal{H} \subset \Phi'$ is a Gelfand triple. We define a map $h: \Phi \to \Phi'$, such that
$$
  h(\varphi) = T_{-\varphi'' + f\varphi} + \varphi(a) \, \delta_a \: ,
$$
where $T_g(\psi) := \int \! g \, \psi$ is the regular distribution of function $g$.  Now we define a sesquilinear form $\omega: \Phi \times \Phi \to \mathbb{C}$, such that
$$
  \omega(\varphi, \psi) = \big< h(\overline\varphi), \; \psi \big>_\Phi
  = \big( h(\overline\varphi) \big)(\psi) \: .
$$
Now, if we extend $\omega$ to a maximum domain (in some sense), it uniquely determines the self-adjoint operator $H$:
$$
  \Omega|_{\Phi\times\Phi} = \omega \\[10pt]
  \Omega(\varphi, \psi)
  = \big( H\varphi, \; \psi \big)_{\mathcal{H}}
  = \big( \varphi, \; H\psi \big)_{\mathcal{H}}
$$

Is there any author that proves this rigorously, or approaches the problem in a similar way? Does this approach extend to potentials with distributions other than $\delta$?