Homework: Feynman diagram loops, confusion about natural units I've got another homework problem that I have trouble with:

4.3 Exercise 3
Use the Feynman rules to show that the expansion in number of loops is an expansion in powers of $\hbar$. Hint: obviously you will have to reinstate the factors $\hbar$ into the expressions for the interaction vertex and the Feynman propagators! You will also need the rule we derived which states that the number of loops (which is the same as the number of undetermined momenta) is equal to the number of propagators in a diagram minus the number of interaction vertices plus one. [So $L = I - V + 1$]

So far I've got:
The factors for internal lines will be proportional to some power $\hbar^x$ (this must be the same power for fermion lines and photon lines, or otherwise the question isn't answerable). I think both factors
$$
\text{i}D_{\text{F}\alpha\beta}(k) = \text{i} \frac{-g_{\alpha\beta}}{k^2 + \text{i}\varepsilon}
$$
and
$$
\text{i}S_F(p) = \text{i}\frac{1}{\not{p}-m+\text{i}\varepsilon}
$$
need to be dimensionless, but I don't see how you could get them to be dimensionless (or even to be the same dimension) by just adding $c$'s and $\hbar$'s. Could someone actually write these expressions with the appropriate factors inserted?
The factors for vertices will be proportional to some power $\hbar^y$. But since the expression here is 
$$
\text{i}e \gamma^\alpha
$$
which carries units of charge, I'm lost again.
Once that's solved, the total S-matrix amplitude of a diagram should be proportional to $(\hbar^x)^I (\hbar^y)^V = \hbar^{xI+yV}$. The question, again, cannot be answered properly unless it happens that $x = -y$ (which I'm therefore assuming will turn out to be the case), in which case you can use $L = I - V + 1$ to see that $\hbar^{xI+yV} = \hbar^{x(I-V)} = \hbar^{x(L-1)}$.
So I'm stuck until I find out how to properly deal with those units
Side question to anyone who's actually working in theoretical physics: are natural units really worth it? So far, for all the slight simplification of written expressions they allow, the dimensional confusion they create has caused me nothing but annoyance.
 A: The easiest way to see how you can recover the powers of $\hbar$ is by looking at the path integral, which (in natural units) is given by
\begin{align}
\mathcal{Z}=\int \left[\mathcal{D}\phi\right]e^{iS}\end{align}
with the action $S$. Going back to dimensionful units, $S$ will have dimension of an action, so the same as $\hbar$. To make the argument of the exponential dimensionless, we therefore need the path integral to be
\begin{align}
\mathcal{Z} = \int \left[\mathcal{D}\phi\right] e^{\frac{i}{\hbar}S},\end{align}
so the action is rescaled by $\frac{1}{\hbar}$. This means that every internal propagator gets a factor of $\hbar$ (because it is the inverse of the kinetic term) and every vertex gets a factor of $\frac{1}{\hbar}$. You can probably continue from here on your own. 
One thing to note is, that $D_F$ and $S_F$ don't need to have the same dimensions, as one of those has the dimensions of a bosonic field squared and the other those of a fermion field squared (remember that they are given as expectation values of time-ordered field products).
A: To explicitly show the dependence on $\hbar$ in the perturbative expansion, we work in units $c = \epsilon_0 = 1$ which means that $\hbar$ has units of energy times length: $[\hbar] = E \cdot L$. From the anti-commutation relations, we can find that,
$$\{ \psi^\dagger(x), \psi(y)\} = \hbar \delta^{(3)}(x-y)$$
from which we deduce $[\psi] = E^{1/2} \cdot L^{-1}$. In the path integral we have $\exp \left(\frac{i}{\hbar}S\right)$ and thus $\frac{1}{\hbar}S$ must be dimensionless. We thus have,
$$Z = \int D[\tilde \psi ]D[\bar{\tilde{\psi}}]D[\tilde A] \exp \left\{ i \int d^4x \left[ \bar{\tilde{\psi}}(i\gamma^\mu \partial_\mu - \tilde e \sqrt{\hbar} \gamma^\mu \tilde A_\mu - \tilde m )\tilde \psi - \frac14 \tilde F_{\mu\nu}\tilde F^{\mu\nu}\right] \right\}$$
so all $\hbar$-dependence is in the coupling, illuminated by the rescaling $\tilde \psi := \psi/\sqrt{\hbar}$ and $\tilde A := A/\sqrt{\hbar}$, as well as $e = \hbar \tilde e$ and similarly for the mass.
It becomes evident in this form, since we known each order in the expansion brings a factor of $\tilde e^2$ into the perturbative expansion that every order has a power of $\hbar$. $n$ loops brings $\alpha^n \propto \hbar^n$.
For more information on this approach see The $\hbar$ Expansion in Quantum Field Theory. For the standard explanation with propagators contributing $\hbar$ and vertices $\frac{1}{\hbar}$ see:


*

*Quantum Field Theory, Zuber et al, section 6.2.1.

*Functional Methods and Perturbation Theory, Rev. Mod. Phys. 47, 165.



Addendum
The more standard approach of this result is to show vertices carry $\frac{1}{\hbar}$ and propagators $\hbar$. In summary, to see why propagators carry $\hbar$, note that if each $\psi$ carries a $\sqrt \hbar$ then the propagator in the free theory, written as a correlation function,
$$D_{\alpha \beta} \sim \langle 0 | \psi_\alpha(x) \bar{\psi}_\beta(y) | 0 \rangle$$
will carry a factor of $\hbar$.
