Lorentz Covariant formula for Noether Charges in QFT

Question

I'm looking for a Lorentz covariant expression of Noether charges and I found this article: https://arxiv.org/abs/hep-th/0701268, section II-A in particular. 

Consider specifically eq. (20-21), they claim: $$Q_\mu=\frac{1}{2}(\phi, P_\mu \rhd \phi),$$ $Q$ is the conserved charge,  "$\rhd$" is just an "acting on" symbol and the inner product is defined by ($\varepsilon$ is the sign function): $$(\phi_1, \phi_2)=\int d^4p \, \delta(p^2-m^2) \varepsilon(p_0)\tilde{\phi}_1^*(-p)\tilde{\phi}_2(p).$$
$\tilde{\phi}$ is the fourier transform of $\phi$.

Hence the Noether Charge is $$\begin{equation} \label{1} \tag{1} Q_\mu = \frac{1}{2}\int d^4p \, \delta(p^2-m^2) \varepsilon(p_0)p_\mu\tilde{\phi}^*(-p)\tilde{\phi}(p).  \end{equation}$$

Now I'm struggling to get the well known quantised expression in QFT: $$Q_\mu= \int d^3p\,\, p_\mu \,\,a^\dagger (\vec{p}) \,a (\vec{p}), \,\,\,\,[a^\dagger (\vec{p}), \,a (\vec{q})]=-\delta^3 (\vec{p}-\vec{q})$$  
from (1), by plugging in usual scalar klein gordon field with creation and annihilation operators.. 

If I'm not wrong (1) in coordinate space looks like $$\int d^4x \,\,d^4y \,\, \phi (x) \Delta (x-y)\left( -i\frac{\partial \phi(y)}{\partial y^\mu}\right)=Q_\mu, \tag{2}$$
where $\Delta$ is the usual commutator function $$\Delta (x-y)=\int d^4p \,\, \varepsilon(p_0) \delta (p^2-m^2)e^{-ip\cdot (x-y)}.$$ Just by substituting in (2) $$\phi(x)= \int \frac{d^3p}{\sqrt{2\omega}_\vec{p}}(a(\vec{p}) e^{-ip\cdot x}+a^\dagger(\vec{p}) e^{ip\cdot x} )$$ I don't seem to be getting the right answer.

Maybe I'm doing some calculation wrong or misinterpreted the article. Any help would be greatly appreciated!

My try:

For instance writing $\phi(x)=\int d^4p \,\, a(p) \delta(p^2-m^2) e^{-ip \cdot x}$ then $\tilde{\phi}(p)=a(p)\delta(p^2-m^2)$ and (1) becomes $$Q_\mu= \frac{1}{2} \int d^4p \,\, \varepsilon(p_0) \, p_\mu a(-p)a(p)\,\,\delta(p^2-m^2)\delta(p^2-m^2)\delta(p^2-m^2).$$
Is this right? How to work out the three deltas? I could use the identity $\delta(x)f(x)=\delta(x)f(0)$ with  $f=\delta$ twice to get $$\delta(p^2-m^2)\delta(p^2-m^2)\delta(p^2-m^2)=\delta(p^2-m^2)\delta(0)\delta(0)=\delta(p^2-m^2)\cdot \mathcal{S},$$ where $\mathcal{S}$ is an (infinite) surface contribution which I currenlty fail to see how cancels out.. what am I missing? Maybe they're using a different convention of their fourier transforms?