Electric field and potential of a point charge in a (strictly) 2D 'world' 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.
 A: To model a 2D system in a 3D world, your quantities should be independent of $z$. In other words, every slice of constant $z$ should be identical. You use the density $\rho(x,y,z)=Q\delta(x)\delta(y)\delta(z)$, which does depend on $z$.
How do we fix this? The answer is simple. Just make the system translation invariant in the $z$ direction. For example $\rho(x,y,z)=\lambda \delta(x)\delta(y)$, which represents a wire with constant line density. The invariant variable is now $x^2+y^2$ (polar coordinates), instead of $x^2+y^2+z^2$. What do you get when you use this? What test surface would you use for a problem with cylindrical symmetry?
A: 
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.

This is your problem.  By the definition of density, you must have
$$
\int \rho(\mathbf{r}) \, d^n \mathbf{r} = Q_\text{enc}
$$
where $Q_\text{enc}$ has the units of charge and $n$ is the number of dimensions.  This implies that dimensionally you must have
$$
[\rho] [L^n] = [Q] \quad \Rightarrow \quad [\rho] = \frac{[Q]}{[L^n]}.
$$
In other words, the charge density in 2-D must have units of $\text{C/m}^2$, not $\text{C/m}^3$.
Note that what does change between dimensions is the dimensions of $\epsilon_0$;  we must also have from Gauss's Law
$$
\oint \mathbf{E} \cdot d^{n-1}\mathbf{a} = \frac{Q_\text{enc}}{\epsilon_0} \quad \Rightarrow \quad [\epsilon_0] = \frac{[Q]}{[E][L^{n-1}]}.
$$
If we want $\mathbf[E]$ to have the same units in all dimensions (which seems reasonable) then the easiest route is to say that the units of $\epsilon_0$ differ between the various dimensions.  So in 3 dimensions it has units of $\text{C}^2/\text{N}\cdot\text{m}^2$;  in 2 dimensions it would instead have units of $\text{C}^2/\text{N}\cdot\text{m}$.  If we modify $\epsilon_0$ in this way between the dimensions, then everything is nice & dimensionally consistent.
