Let $f(z) = \log|F(z)|$, where $F: \mathbb{C} \rightarrow \mathbb{C}$ is entire. Then $f$ defines a distribution on $\mathbb{R}^2$, and we want to show that its distributional Laplacian is
$$\Delta f = \sum_{n \in I} 2 \pi \, d_n \, \delta_{z_n}$$
where $\{z_n \mid n \in I\}$ are the zeros of $F$, which have respective degree (multiplicity) $d_n$.
This is what I have got so far (though I suspect that there could be a completely different approach which is easier, so solutions that don't use this work are more than welcome!):
Let $U_\varepsilon = \mathbb{R}^2 \backslash \left( \bigcup\limits_{n \in I} B(z_n , \varepsilon) \right)$.
Then, for a test function $\varphi \in D(\mathbb{R}^2)$,
$$\left\langle f,\Delta\varphi\right\rangle=\lim_{\varepsilon\rightarrow0^+}\int_{U_\varepsilon} f\Delta\varphi \ dx,$$ and by Green's second identity $$\int_{U_\varepsilon} f\Delta\varphi \ dx=\int_{U_\varepsilon} \varphi\Delta f \ dx+\int_{\partial U_\varepsilon} f\frac{\partial\varphi}{\partial n}ds-\int_{\partial U_\varepsilon}\varphi\frac{\partial f}{\partial n} \ ds.$$
The first integral on the RHS vanishes since $f = Re(\log(F))$, so is harmonic on $U_\varepsilon$ where $\log(F)$ is holomorphic.
We can obviously split $\partial U_\varepsilon$ into circles around each singularity of $f$.
This is where I get stuck. How can we deal with the terms $f\frac{\partial\varphi}{\partial n}$ and $\frac{\partial f}{\partial n}$.
Any suggestions would be appreciated!
About the term $\int_{\partial U_\varepsilon} f\frac{\partial\varphi}{\partial n}ds$, note that $\partial U_\epsilon = \bigcup_n \partial B(z_n,\epsilon)$. Fix $n$ and observe that $$ F(z) = (z-z_n)^d G(z) $$ where $G(z)$ does not vanish near $z=z_n$. This implies $$ f(z) = d\log |z-z_n| + \log|G(z)|=d\log \epsilon + \log|G(z)| $$ on $|z-z_n|=\epsilon$. Note that $\log |G(z)|$ is a continuous harmonic function and hence $|\log |G(z)||\leq C$ for some $C>0$. Now, we see that $$ \int_{|z-z_n|=\epsilon} f\frac{\partial\varphi}{\partial n}ds =d\log\epsilon \int_{|z-z_n|=\epsilon} \langle \nabla \varphi,n\rangle ds +\int_{|z-z_n|=\epsilon} \log|G(z)|\langle \nabla \varphi,n\rangle ds. $$ From $$ d\Big|\log\epsilon \int_{|z-z_n|=\epsilon} \langle \nabla \varphi,n\rangle ds\Big| = d \log\epsilon^{-1} \int_{|z-z_n|=\epsilon} \|\nabla \varphi\| ds\leq d \log\epsilon^{-1} \cdot 2\pi\epsilon \cdot\sup\|\nabla \varphi\|\to 0, $$ and $$ \Big|\int_{|z-z_n|=\epsilon} \log|G(z)|\langle \nabla \varphi,n\rangle ds\Big|\leq 2\pi C\epsilon \sup\|\nabla\varphi\|\to 0, $$ we get $$ \lim_{\epsilon\to 0}\int_{|z-z_n|=\epsilon} f\frac{\partial\varphi}{\partial n}ds =0. $$ This holds for all $n$, and we conclude $\lim_{\epsilon\to 0}\int_{\partial U_\epsilon} f\frac{\partial\varphi}{\partial n}ds=0.$
Next, note that $$\frac{\partial f(z_n + re^{i\theta})}{\partial n} = \frac{\partial}{\partial r} f(z_n + re^{i\theta})= \frac{\partial}{\partial r}[d\log r + \log |G(z_n+re^{i\theta})|] = \frac{d}{r}+\partial_r\log |G(z_n+re^{i\theta})|. $$ This is because the unit outward normal vector $n$ of the circle $|z-z_n|=r$ at $z = z_n +re^{i\theta}$ is just $e^{i\theta}$. Now observe that $$ \int_{|z-z_n|=\epsilon} \varphi\frac{\partial f}{\partial n}ds = \frac{d}{\epsilon}\int_{|z-z_n|=\epsilon} \varphi ds + \int_{|z-z_n|=\epsilon} \varphi \left(\partial_r|_{r=\epsilon}\log |G(z_n +re^{i\theta})|\right)ds . $$ The first term is $$ \frac{d}{\epsilon}\int_{|z-z_n|=\epsilon} \varphi ds= \frac{d}{\epsilon}\int_0^{2\pi} \varphi(z_n + \epsilon e^{i\theta} )\epsilon d\theta =d\int_0^{2\pi} \varphi(z_n + \epsilon e^{i\theta} )d\theta \to 2\pi d \varphi(z_n). $$ The integrand in the second term is a continuous function and thus $$ \Big|\int_{|z-z_n|=\epsilon} \varphi \left(\partial_r|_{r=\epsilon}\log |G(z_n +re^{i\theta})|\right)ds\Big| \leq \sup |\varphi|\cdot \sup|\partial_r \log |G||\cdot 2\pi\epsilon \to 0. $$ Summing up over all $n$, we see that $$ \lim_{\epsilon\to 0}\int_{U_\varepsilon} f\Delta\varphi \ dx = \sum_n 2\pi d_n \varphi(z_n), $$ as desired.