Let's examine the one-dimensional three-point stencil case in detail,
because I think it's important to be clear just how this behaviour
arises, and what it means to set a point to a certain value in a
finite-difference grid when the underlying function is
discontinuous. The equation will be
$$ u''(x) = \rho(x). $$
Instead of using the interval $[-1,1]$ with discontinuity at
$\frac12$, I will use the interval $[-\frac12,\frac12]$ with the
discontinuity placed at $0$. The grid size will be $h$, and I will only have to consider the interval $[-\frac h2,\frac h2]$ around the grid point $0$.
First, in a finite-difference approximation, we approximate
$$ u''(0) = \frac{\hat u(-h) - 2\hat u(0) + \hat u(h)}{h^2} $$
and solve
$$ \frac{\hat u(-h) - 2\hat u(0) + \hat u(h)}{h^2} = \hat \rho(0). $$
(Here the variables with hats are the numerical approximations on the
grid to the variables without hats.)
But this is a very bad specification of the problem, because your
function $\rho(0) = 2$ at $x=0$, and $\rho(x) = 0$ everywhere else, is
discontinuous. In particular, if we shift the grid point $0$ to either
side by some tiny amount $\epsilon$, the numerical solution changes
entirely and becomes exactly zero.
This means that this is a terribly misspecified problem.
We can make sense of it by
converting it to an equivalent finite-volume formulation, where it
will make much more sense.
In a finite-volume method, we solve
$$ \int_{-h/2}^{h/2} u''(x)\,dx = \int_{-h/2}^{h/2}\rho(x)\,dx, $$
by picking $u(x)$ to be a suitable approximation to the unknown
function. Let's pick, on the interval $[-\frac h2,\frac h2]$ the approximation
$$ u(x) = \hat u(-h) \phi(x+h) + \hat u(0) \phi(x) + \hat u(h) \phi(x-h), $$
where $\phi(x)$ is the basis function
$$ \phi(x) = \max\left(1-\frac{|x|}{h}, 0\right), $$
(it's a piecewise linear function that goes from $0$ at $-h$ to $1$ at
$0$ to $0$ at $+h$, thus interpolating between grid points.). The approximation $u(x)$ is a weighted sum of three basis functions that look like this:
We can then compute
$$ \phi''(x) = \frac1h\delta(h-|x|) - \frac2h \delta(x), $$
so that the approximation to the integral is
$$ \int_{-h/2}^{h/2} u''(x) = \frac{\hat u(-h) -2\hat u(0)+\hat
u(h)}{h}, $$
and the finite-volume approximation to our equation becomes
$$ \frac{\hat u(-h) -2\hat u(0) + \hat u(h)}{h} = \int_{-h/2}^{h/2}
\rho(x)\,dx = h \hat \rho(0). $$
There are two important things here. First, this is equivalent to
the finite-difference formulation in that we end up solving the same
equations. Second, the discontinuity in $\rho$ is given a very precise
meaning: when we use the value $2$ for $\hat \rho(0)$, we are saying
that this is the average value of $\rho$ on the interval $[-h/2,h/2]$:
$$2 = \hat\rho(0) = \frac1h \int_{-h/2}^{h/2} \rho(x)\,dx.$$
This interpretation is not available in the finite-difference
formulation. It is also not so sensitive to the location of the discontinuity: if the discontinuity were at some small distance $\epsilon$ away from $0$, the average value would be almost the same, but the value at $0$ might be completely different.
But if we say that $2$ is the average value of $\rho$ near the grid
point, we can then go back and compute the exact solution of the
equation with the right-hand side given by
$$ \tilde \rho(x) = 2[-h/2 < x < h/2]. $$
(We pick a function of our own choice that gives the right average.)
In this case, the exact solution will be
$$ \tilde u(x) = 2 \int_{-h/2}^{h/2} G(x; u)\,du \approx 2h G(x;
0), $$
in terms of the Green function for the Poisson equation. In the two-dimensional case it will be
$$ \approx 2h^2 G(x, y; 0, 0), $$
as in the other (correct) answer on MSE.
Finally, the outcome of all this is that when you say that you compare
your numerical approximate solution with the exact solution $u(x)=0$,
this is wrong. The exact solution should not be zero, it should be
$$ \approx 2h^2 G(x, y; 0,0). $$
It therefore should make perfect sense that the fourth-order solution does not converge to zero with order $4$: it should converge to the correct solution, which is not zero but has magnitude of order $O(h^2)$.
If you do want to get zero as the numerical solution and compare with the mathematical solution, you should set $\hat\rho(0)$ to be the average value of $\rho$, which is $0$, not $2$.
The fact that the exact solution depends on the chosen grid size indicates that this is not a good way to check whether you implemented the method correctly. A very straightforward technique known as the method of manufactured solutions is better for this.
