Alright …
So the Lagrangian density
\begin{eqnarray}
\mathcal{L} \propto (\partial \phi)^2
\end{eqnarray}
should be integrated over ##d^4x## to give units of action. The Lagrangian density has units
\begin{eqnarray}
[(\partial\phi)^2] = \frac{1}{l^2} [\phi]^2.
\end{eqnarray}
And
\begin{eqnarray}
[d^4x] = l^4.
\end{eqnarray}
So, this gives
\begin{eqnarray}
[d^4x \mathcal{L}] = l^4 \frac{1}{l^2} [\phi]^2 = l^2 [\phi]^2 \equiv \frac{m l^2}{t},
\end{eqnarray}
or
\begin{eqnarray}
[\phi]^2 = \frac{m}{t}.
\end{eqnarray}
So the field has units of the square root of mass per time ... ?
For ##a(\boldsymbol{p})##, I hadn't thought of the commutation relation - I'm actually only interested in in the classical (real) solutions at this point. But,
\begin{eqnarray}
[\delta(\boldsymbol{p})] = \frac{1}{[p]^3} = \bigg(\frac{t}{ml}\bigg)^3
\end{eqnarray}
so as you said,
\begin{eqnarray}
[a] = \bigg( \frac{t}{ml}\bigg)^{\frac{3}{2}}.
\end{eqnarray}
So, from the natural units expression above,
\begin{eqnarray}
\phi = \alpha \int \frac{d^3 p}{(2\pi)^3 2E_p} \big( a(\boldsymbol{p}) e^{\cdots} + \cdots),
\end{eqnarray}
where ##\alpha## captures whatever units are missing. This gives
\begin{eqnarray}
\sqrt{\frac{m}{t}} &=& [\alpha] \bigg( \frac{ml}{t}\bigg)^3 \frac{1}{ \frac{ml^2}{t^2}} \bigg( \frac{t}{ml}\bigg)^{\frac{3}{2}} \\
&=& [\alpha] \sqrt{ \frac{mt}{l} }.
\end{eqnarray}
So
\begin{eqnarray}
[\alpha] = \frac{ \sqrt{l}}{t}.
\end{eqnarray}
But when I try to find the appropriate factor of ##\hbar## and ##c## (assuming ##G## doesn't enter in - that would be weird), I try
\begin{eqnarray}
[\hbar^ac^b] = \frac{m^a l^{2a+b}}{t^{a+b}}.
\end{eqnarray}
Setting this equal to the units for ##\alpha##, ##a## must be ##0##, but then there is no solution for ##b##.
Not sure what I've done wrong.