I was trying to figure out how the electric potential and electric field are different in a 3D system versus in a 2D system (I take such a 2D 'world' to be the $xy$-plane, i.e. $z=0$, in a Cartesian coordinate system). The context/motivation has to do with the Coulomb interactions that electrons/holes experience in a 2D or a 3D metal, how this is different and how this results in the Coulomb screening being different between the two. (I want to emphasise that such a 2D material is a truly 2D system, i.e. the charges living in such a material do not know about a third dimension.)
I have been looking on the internet a lot, and I know that the answer has to be something like $$V_\text{2D}(r) \sim \log(r) \qquad E_\text{2D}(r) \sim \frac{1}{r},$$ and of course $$V_\text{3D}(r) \sim \frac{1}{r} \qquad E_\text{3D}(r) \sim \frac{1}{r^2},$$ but when trying to reproduce the 2D result I am struggling. (PS. I am using SI units)
Starting with Gauss' law,
$$\nabla \cdot \mathbf{E}\left(\mathbf{r}\right) = \frac{\rho\left(\mathbf{r}\right)}{\varepsilon_0}$$
I transformed it to its integral form in 2D and 3D using the divergence theorem in "$n$" dimensions, $$\int_U \text{d}^{n}{r} \, \nabla\cdot\mathbf{A} = \oint_{\partial U} \text{d}^{n-1}r \, \mathbf{A} \cdot \mathbf{\hat{n}}.$$
which gives
$$ \oint_{\partial S} \text{d}r \, \mathbf{E}\left(x,y\right) \cdot \mathbf{\hat{n}} = \frac{1}{\varepsilon_0} \int_S \text{d}^2 r \, \rho\left(x,y\right), \quad \text{(2D)}$$ $$\oint_{\partial V} \text{d}^2 r \, \mathbf{E}\left(x,y,z\right) \cdot \mathbf{\hat{n}} = \frac{1}{\varepsilon_0} \int_V \text{d}^3 r \, \rho\left(x,y,z\right). \quad \text{(3D)}$$
Here I explicitly show that $\mathbf{E}$ and $\rho$ are in 2D and 3D. PS. in 2D, $S$ is the surface you integrate over and $\partial S$ is its boundary 'line'; in 3D, $V$ is the volume you integrate over and $\partial V$ is its boundary surface
Now suppose we have a point charge $Q$ at the origin, i.e. $$\rho\left(\mathbf{r}\right) = Q \, \delta^{(3)}\left(\mathbf{r}\right)$$ (where $\delta^{(3)}$ is the Dirac delta function in 3 dimensions, i.e. with units $\text{m}^{-3}$). I am using this in 2D and in 3D, since $\rho\left(\mathbf{r}\right)$ must have dimensions of $\text{C} \, \text{m}^{-3}$ by virtue of Gauss' law, regardless the dimensionality.
In 3D, we use as the integration volume a sphere with radius $r$ centred at the origin (just like any elementary electromagnetism textbook does it).
In 2D, we use as the integration surface a circle with radius $r$ centred at the origin $(x,y)=(0,0)$ (remember, the 2D 'world' is the $xy$-plane, i.e. $z=0$ plane).
By symmetry considerations, it's clear that in 2D and in 3D the field only has a radial component. I will mark this as $E_r$ in both cases, although I want to emphasise that "$r$" in 2D means $\sqrt{x^2+y^2}$ while in 3D it means $\sqrt{x^2+y^2+z^2}$. This gives
$$ 2\pi r E_r = \frac{Q}{\varepsilon_0} \int_S \text{d}^2 r \, \delta^{(3)}\left(\mathbf{r}\right) = \frac{Q}{\varepsilon_0} \delta\left(z\right), \quad \text{(2D)}$$ $$ 4 \pi r^2 E_r = \frac{Q}{\varepsilon_0} \int_V \text{d}^3 r \, \delta^{(3)}\left(\mathbf{r}\right) = \frac{Q}{\varepsilon_0}, \quad \text{(3D)}$$
where I made use of $\delta^{(3)}\left(\mathbf{r}\right) = \delta(x)\delta(y)\delta(z)$. Rewriting gives my final result for the electric field:
$$ \mathbf{E}\left(x,y\right) = \frac{Q}{2\pi\varepsilon_0} \frac{1}{r} \delta\left(z\right) \mathbf{\hat{r}}, \quad \text{(2D)}$$ $$ \mathbf{E}\left(x,y,z\right) = \frac{Q}{4 \pi \varepsilon_0} \frac{1}{r^2} \mathbf{\hat{r}}, \quad \text{(3D)}$$
Disturbingly, I'm left with a Dirac delta function in my electric field, although I do get the $E_\text{2D} \sim 1/r$ dependence I expected. I guess it makes some sense that the $\delta$-function is there, since I'm only concerned with $z=0$ anyways, but it of course shouldn't be the case that a physical quantity like the electric field is merely a mathematical 'distribution' i.e. delta function. However, I cannot just leave out the $\delta$-function, since it carries a unit of $1/\text{meter}$ which is crucial to make the units consistent.
What am I doing wrong here? As you can see, I tried to do this derivation rigorously by starting from Gauss' law in coordinate-free form, and then converting using the divergence theorem; in this step, the difference between 2D and 3D should arise. I have seen other questions on this website where people would instead simply say that the right-hand side becomes the "enclosed charge" $$Q_\text{encl} = \int_V \text{d}^3 r \,\rho\left(\mathbf{r}\right)$$ even in 2D (e.g. this answer), but they end up (as far as I can see) with inconsistent units (i.e. they miss a $1/\text{meter}$ which I have in my $\delta\left(z\right)$). Another common method of getting rid of this I see is that people use a 'line charge density' $\lambda$ (units $\text{C} \, \text{m}^{-1}$) (e.g. this answer), but that is only if you have an infinitely long line with a charge per unit length $\lambda$. This shouldn't be valid for a point charge I think.
I hope anyone can help me here. If anything is unclear, please let me know and I will try to clarify it.