Inverting the Laplacian

I've had a hard time looking for literature on this, so here's my question:

We take a look at the Laplacian $$-\Delta$$ as an unbounded operator on $$\mathrm{L}^2(\mathbb{R}^3)$$. We know that $$-\Delta$$ is unitary equivalent to the multiplication operator with $$|p|^2$$ in Fourier space, so $$-\Delta=\mathcal{F}^{-1} |p|^2 \mathcal{F}$$.

So one could now use functional calculus and define an inverse Laplacian by setting $$(-\Delta)^{-1}=\mathcal{F}^{-1} (1/|p|^2 )\mathcal{F}$$, which will also be unbounded of course.

It is also known from theory of PDEs that one can invert the Laplacian on Schwartz functions using the Green's function $$1/4\pi |x|$$, which is derived as a distributional Fourier transform of $$1/|p|^2$$. So define $$G$$ on $$\mathrm{L}^2(\mathbb{R}^3)$$ by convolution with $$1/4\pi |x|$$: $$(G\phi)(x)=\int \frac{\phi(y)}{4\pi |x-y|} dy$$ $$G$$ is also unbounded and coincides at least for the Schwartz functions with the above defined inverse Laplacian, $$G\phi=(-\Delta)^{-1}\phi$$ for all $$\phi \in\mathcal{S}(\mathbb{R}^3)$$.

Now my question is: Are both operators the same? Does $$G=(-\Delta)^{-1}$$ hold, i. e. is $$D(G)=D((-\Delta)^{-1})$$ and $$G\phi=(-\Delta)^{-1}\phi$$ for all $$\phi \in D(G)=D((-\Delta)^{-1})$$?

My guess is that $$\mathcal{S}(\mathbb{R}^3)$$ is a core of $$(-\Delta)^{-1}$$, and as $$(-\Delta)^{-1}|_{\mathcal{S}(\mathbb{R}^3)}=G|_{\mathcal{S}(\mathbb{R}^3)}$$ equality should follow by closing the restrictions.

Any comments, hints on how to proceed or references are welcome! Thank you.

The operator $$(-\Delta +\epsilon I)^{-1} : L^2(\mathbb{R}^3)\rightarrow W^{2,2}(\mathbb{R}^3)$$ is a bicontinuous bijection that is equivalently defined by the Fourier transform $$\mathscr{F}$$ as $$(-\Delta +\epsilon I)^{-1}f = \mathscr{F}^{-1}\frac{1}{|\xi|^2+\epsilon}(\mathscr{F}f)(\xi),\;\;\; f\in L^2(\mathbb{R}^3).$$ Suppose $$f\in L^2$$. If $$g(\xi)=|\xi|^{-2}(\mathscr{F}f)(\xi)$$ is also in $$L^2$$, then $$g=L^2\mbox{-}\lim_{\epsilon\downarrow 0}(-\Delta+\epsilon I)^{-1}f=\mathscr{F}^{-1}\frac{1}{|\xi|^2}\mathscr{F}f = \frac{1}{4\pi}\int_{\mathbb{R}^3}\frac{1}{|x-y|}f(y)dy.$$ It is also true that, for such an $$f$$, the right side of the following converges to $$f$$ in $$L^2$$ as $$\epsilon\downarrow 0$$: $$-\Delta(-\Delta+\epsilon I)^{-1}f=f-\epsilon(-\Delta+\epsilon I)^{-1}f.$$ And $$(-\Delta+\epsilon I)^{-1}f$$ converges to $$\mathscr{F}^{-1}(|\xi|^{-1}\mathscr{F}f)$$. Because $$-\Delta : W^{2,2}\subset L^2\rightarrow L^2$$ is selfadjoint, it follows that $$\mathscr{F}^{-1}(|\xi|^{-2}\mathscr{F}f)\in\mathcal{D}(-\Delta)$$ and $$-\Delta\left[\mathscr{F}^{-1}\frac{1}{|\xi|^2}\mathscr{F}f\right] = f,$$ which is equal to $$-\Delta \frac{1}{4\pi}\int_{\mathbb{R}^3}\frac{1}{|x-y|}f(y)dy=f.$$