When dealing with constraint systems, we use dirac bracket instead of poisson bracket. In that procedure, we first find constraints $\Lambda_i$ and gauge $\Omega_i$, then we calculate the matrix consisting with their commutation, and the inverse of the matrix. It might be complicated.
My question is about a specific example. When dealing with classical electromagnetic field, we need to calculate the inverse of $\nabla^2\delta(x-y)$, but I'm stuck here. How to calculate its inverse?
Let $x,x',y$ denote position vectors.
Let $a(x,x')$ be the inverse of $\nabla^2\delta(x-y)$ in the following sense:
$$\int_\mathcal{M} dx\, a(x,x')\nabla^2\delta(x-y)=\delta(y-x').\tag{1}$$
Integrating the LHS by parts twice and getting rid of boundary terms, as $a(x,x')$ is expected to vanish there, we are left with
$$\int_\mathcal{M} dx\, a(x,x')\nabla^2\delta(x-y)=\int_\mathcal{M} dx\, \nabla^2\Big(a(x,x')\Big)\delta(x-y)=\nabla_y^2 \,a(y,x'),\tag{2}$$
where $\nabla_y^2$ means Laplacian with respect to the $y$ variables. The equation $(1)$ becomes $$\nabla_y^2a(y,x')=\delta(y-x').\tag{3}$$
There is a well known result $$\nabla^2\left(\frac{1}{|x-x'|}\right)=-4\pi\delta(x-x'),\tag{4}$$
so the solution to $(3)$ is $$a(y,x')=\frac{-1}{4\pi|y-x'|}.$$
