This equation will always give you a volume charge density. One way to see this is that surface charge density and volume charge density have different units - $\mathrm{C/m^2}$ and $\mathrm{C/m^3}$ respectively - and in order for the units to be consistent, $\rho$ has to be the latter. The fact that the equation is written with $\rho$ is a helpful reminder that it is a volume charge density.
Of course, keep in mind that the potential is not $kx^{4/3}$ everywhere. That function only describes the potential within a certain region. You also have to think about what's happening outside that region, and on the boundaries of the region.
If you try solving Poisson's equation $\nabla^2\varphi = -\rho/\epsilon_0$ in region where the potential is not so nicely behaved (as you have to do here, if you think about the boundaries), you might get a solution that involves a delta function. Just to pull an example out of thin air, something like
$$\rho(x, y, z) = \delta(x - L) e^{-y^2 - z^2}$$
That is the signature of a surface charge density being expressed as a volume charge density. $\sigma$ is the part other than the delta function; in general:
$$\rho(x, y, z) = \delta(x - a)\sigma(y, z)$$
so in this purely hypothetical example you could discern that $\sigma(y, z) = e^{-y^2 - z^2}$.
This is consistent with the statement that surface charge densities correspond to discontinuities in the electric field, because remember you can write Poisson's equation as
$$\vec\nabla\cdot\vec{E} = \frac{\rho}{\epsilon_0}$$
When $\vec{E}$ is discontinuous, its derivative is "infinite", and therefore $\rho$ needs to be represented as a product involving a delta function.
I think you are on the right track. The idea is to consider the field of a complete, unpunctured sphere of surface charge density $\sigma$ and "add" a little patch of surface charge density $-\sigma$. By superposition, the patch and the full sphere is equivalent to a sphere with a small hole in it.
You are correct that the field right at the surface for the unpunctured sphere is $E=\sigma/\epsilon_0$. Now the trick I believe is that, for a point very close to the surface, the patch can be considered as an infinite plate of surface charge $-\sigma$. This is possible on the same grounds that one considers a plate infinite in extent if the point where we want to field is much closer to the plate than any of the physical dimensions of the plate. This seems justified here as you are at the centre of the "plate" and arbitrarily close to it. The field of such a plate is $-\sigma/2\epsilon_0$ so the net field at the centre is thus
$$
\frac{\sigma}{\epsilon_0}-\frac{\sigma}{2\epsilon_0}=\frac{\sigma}{2\epsilon_0}
$$
Best Answer
I have personally also contemplated this issue, and have come up with a simple solution that is satisfactory, to me at least. I'm sure this can also be found in many textbooks. In general, we have
$$\tag{$\star$} Q=\int \rho\ d\tau$$
because we are considering a three dimensional space. Intuitively, we feel it should be possible to talk about a three dimensional charge distribution in every case. The question is how to conceptualize this when discussing surface, line or point charges. The solution comes in the form of the Dirac delta distribution (or function, depending who you ask).
Let's take a look at an example: Consider a 2-sphere of radius $R$, with some charge distribution $\sigma(\theta,\phi)$ on it. What is the three dimensional charge distribution $\rho(r,\theta,\phi)$ corresponding to this situation? Like I said, we have to use the Dirac delta:
$$\rho(r,\theta,\phi)=\delta(r-R)\sigma(\theta,\phi)$$
Now, $(\star)$ gives us: $$ Q=\int\rho\ d\tau=\int_{0}^{2\pi}\int_{0}^\pi\int_0^\infty \rho\ r^2\sin\theta\ dr\ d\theta\ d\phi=R^2\int_{0}^{2\pi}\int_{0}^\pi\sigma\ \sin\theta\ d\theta\ d\phi$$
Similarly, when considering a line or point charge, one uses two or three Dirac delta's to describe the distribution in 3-space.