In the complex plane, one can write $$\delta(x) \delta(y) = \dfrac{1}{\pi} \partial_z \dfrac{1}{\bar{z}}.$$
How to prove this relation has been answered before and identifies
$$ 2\pi \delta(x) \delta(y) = \partial_x \left(\dfrac{x}{x^2 + y^2} \right) + \partial_y \left(\dfrac{y}{x^2 + y^2} \right). $$
However, I have not been able to show that this relation indeed holds. How can I prove that it is indeed true?
Let me write $z = (x,y)$ and $|z|^2 = x^2+y^2$. Then, your second relation can be written $$ \nabla\cdot(z/|z|^2) = 2\pi\,\delta(z) $$ where $\delta(z) = \delta(x)\,\delta(y)$ is the Dirac delta in $\Bbb R^2$. This follows from the well-known fact that $-\ln |z|$ is a fundamental solution of the Poisson equation, since then $$ \nabla\cdot(z/|z|^2) = \nabla\cdot(\nabla(\ln|z|)) = \Delta \ln |z| = 2\pi\,\delta(z). $$
There are a lot of derivations of the Poisson fundamental solution. One can remark that taking the derivative in the classical sense $\Delta \ln |z| = 0$ for $z\neq 0$, use a scaling argument to deduce that $\Delta \ln |z| = C\,\delta(z)$. Another that I prefer is to use the Fourier transform: since (see e.g. here) $\mathcal{F}(\ln(z)) = C\,\delta(z) - \mathrm{pf}\frac{1}{2\pi\,|z|^2}$, one deduces that $$ \mathcal{F}(\Delta\ln(z)) = -C\,4\pi^2 |z|^2\,\delta(z) + 2\pi = 2\pi $$ because $|z|^2\,\delta(z) = 0$, and so $\Delta\ln(z) = 2\pi\,\delta(z)$.
The strategy is to show that the delta function at the origin in the plane can be approximated by a family of smooth density functions that depend on a parameter $\epsilon \to 0$ that ultimately concentrate their mass at the origin.
First make use of the approximation $\frac{ 1}{\overline z} = \frac{ z}{ z\overline z} \sim f(z,\overline z)=\frac{ z}{ z\overline z+\epsilon}$ and then we use the Wirtinger calculus and quotient rule to compute
$ \partial_z f(z,\overline z)=\partial_z \frac{ z}{z\overline z+ \epsilon} =\frac{ (\epsilon + z\overline z)-z\overline z }{(|z|^2+ \epsilon)^2} =\frac{\epsilon}{(|z|^2+ \epsilon)^2} :=\mu(z,\overline z)$.
Next check that this mass density expression always has the same total mass in the plane; that is, show that $I=\int \int \mu \ dx \ dy$ is always the same constant regardless of the choice of $\epsilon$. (This is easy to verify by making the scaled change of variables $x=\sqrt{\epsilon} x', y= \sqrt{\epsilon} y'$ for which $ dx dy = \epsilon dx' dy'$.)
Next you can evaluate this normalization constant by setting $\epsilon =1$ and then using polar coordinates:
$I= \int_0^{2\pi} \int_{r=0}^{\infty} \frac {r \ dr \ d\theta}{ (r^2+ 1)^2} = \pi$
Lastly, (handwaving somewhat) you can confirm by graphing the density functions that the mass density distributions concentrate their mass in zones whose radii ultimately tend to zero.
Footnote. The function $f$ and the measure $\mu$ are connected to the Laplacian of perturbations of the log function as follows: $ \mu=\frac{\partial^2 \ln(\epsilon+ z \overline z)}{\partial z \partial\overline z}$ and $f= \frac{ \partial}{\partial \overline z} \ln(\epsilon+ z \overline z)$.
The usual real-variable Laplacian can be written in Wirtinger notation as $ u_{xx}+u_{yy}= 4 u_{z\overline z}$.
Another strategy starts by applying Stokes' theorem to Cauchy's integral formula : $$ f(\zeta) = \oint\frac{\mathrm{d}z}{2\pi i} \frac{f(z)}{z-\zeta} = \int\frac{\mathrm{d}\bar{z}\wedge\mathrm{d}z}{2\pi i} \frac{\partial}{\partial\bar{z}}\frac{f(z)}{z-\zeta} $$ With Cartesian coordinates, i.e. $z = x + iy$, one has $\mathrm{d}\bar{z}\wedge\mathrm{d}z = 2i \,\mathrm{d}x\wedge\mathrm{d}y$, hence the following representation of the 2D Dirac delta on the (Cartesian) complex plane : $$ \delta^{(2)}(x,y) = \delta(x)\delta(y) = \partial_{\bar{z}}\frac{1}{\pi z} $$ and, by taking the complex conjugate of this relation, $$ \delta(x)\delta(y) = \partial_z\frac{1}{\pi\bar{z}}. $$ As for your second identity, it can be proven by interpreting $\partial_z$ as a Wirtinger derivative, as follows : $$ \begin{array}{rcl} 2\pi\delta(x)\delta(y) &=& \displaystyle 2\partial_z\frac{1}{\bar{z}} \quad\mathrm{with}\quad \frac{1}{\bar{z}} = \frac{z}{|z|^2} \\ &=& \displaystyle (\partial_x - i\partial_y) \frac{x+iy}{x^2+y^2} \\ &=& \displaystyle \partial_x\frac{x}{x^2+y^2} + \partial_y\frac{y}{x^2+y^2} + i\partial_x\frac{y}{x^2+y^2} - i\partial_y\frac{x}{x^2+y^2} \\ &=& \displaystyle \partial_x\frac{x}{x^2+y^2} + \partial_y\frac{y}{x^2+y^2} \end{array} $$ because the last two terms compensate each other.
