Quantum corrections in path integral I am working the following exercise:

Calculate the generating functional $$Z[j]=\int \mathcal{D}\Phi \exp\left(\frac{i}{\hbar}S[\Phi,j]\right),\quad S[\Phi,j]=\int d^4x(\mathcal{L}(\Phi)+j\Phi),$$ $$\mathcal{L}=\frac{1}{2}\partial_\mu\Phi\partial^\mu\Phi-\frac{1}{2}m^2\Phi^2+V(\Phi)$$
  using the method of stationary phase for fluctuations $\delta\Phi$ around the classical field configuration $\Phi_c$, for which $\delta S/\delta\Phi|_{\Phi_c}=0$. Show that this approximation is an asymptotic series in $\hbar$.

Now let me present you my calculations.
The method of stationary phase says that
$$\int d^nx e^{i\lambda S(x)}\simeq \sqrt{\frac{2\pi i}{\lambda^n \det S^{\prime\prime}_0}}e^{i\lambda S_0},$$
where $S_0$ denotes the extremum.
Generalising this to functional integrals, this gives
$$\int \mathcal{D}\Phi \exp\left(\frac{i}{\hbar}S[\Phi,j]\right)\simeq N(n,\hbar)\sqrt{\det S^{\prime\prime}_0}^{-1}e^{i\lambda S_0},$$
where $S^{\prime\prime}_0=\delta^2S/\delta\Phi^2|_{\Phi_c}=\left(V^{\prime\prime}(\Phi_c)-(\square_x+m^2)\right)\delta(x-y)$.
Now, I don't see how this gives us any asymptotic series in $\hbar$.
My first idea was to do the following instead. We substitute $\Phi=\Phi_c+\Delta$, which gives $\mathcal{D}\Phi=\mathcal{D}\Delta$. Then we insert the substitution into the action and obtain three contributions:
$$S[\Phi,j]=\int d^4x\left(\frac{1}{2}\partial_\mu\Phi_c\partial^\mu\Phi_c-\frac{1}{2}m^2\Phi^2_c+V(\Phi_c)+j\Phi_c\right)\\
+\int d^4x\left(V^\prime(\Phi_c)+j-(\square+m^2)\Phi_c\right)\Delta\\
+\int d^4x\left(\frac{1}{2}\partial_\mu\Delta\partial^\mu\Delta-\frac{1}{2}m^2\Delta^2\right).$$
The middle term vanishes, since $\delta S/\delta\Phi|_{\Phi_c}=0$ and thus $(\square+m^2)\Phi_c-V^\prime(\Phi_c)=j$, i.e. the classical field fulfills the EOM with source term.
This then gives us
$$Z[j]=\int \mathcal{D}\Phi \exp\left(\frac{i}{\hbar}S[\Phi,j]\right)\\
=\exp\left(\frac{i}{\hbar}S[\Phi_c,j]\right)\int\mathcal{D}\Delta\exp\left[\frac{i}{\hbar}\int d^4x\left(\frac{1}{2}\partial_\mu\Delta\partial^\mu\Delta-\frac{1}{2}m^2\Delta^2\right)\right],$$
which can then be solved as a Gaussian integral:
$$Z[j]\simeq\exp\left(\frac{i}{\hbar}S[\Phi_c,j]\right)\times N(n,\hbar)\sqrt{\det(\square+m^2)}^{-1}.$$
But again this is not an asymptotic series in $\hbar$.
To me it seems, we are only ever going to get an expansion in terms of $\hbar$, if we expand out the exponential. But then we cannot use the method of stationary phase anymore. However, if we use the method of stationary phase, we will only ever get $\hbar$ to appear in the divergent prefactor, which is divided out by normalisation anyway.
How do I make this work?
 A: We have 
$$S[\Phi_c +\sqrt\hbar\Delta]=S[\Phi_c]+\sqrt\hbar\int dx\ \overbrace{\frac{\delta S}{\delta\Phi}[\Phi_c(x)]}^{=0}\Delta(x)+\frac{\hbar}{2}\iint dxdy\ \frac{\delta^2 S}{\delta \Phi(x)\delta\Phi(y)}\Delta(x)\Delta(y)+\ldots$$
So in particular 
$$S[\Phi,j]=\int d^4x\left(\frac{1}{2}\partial_\mu\Phi_c\partial^\mu\Phi_c-\frac{1}{2}m^2\Phi^2_c+V(\Phi_c)+j\Phi_c\right)\\
+\sqrt\hbar\int d^4x\left(j-(\square+m^2)\Phi_c+\right)\Delta\\
+\frac{\hbar}{2}\int d^4x\left(\frac{1}{2}\partial_\mu\Delta\partial^\mu\Delta-\frac{1}{2}m^2\Delta^2\right)\\
+\sum_{n=1}^\infty\frac{\hbar^{n/2}}{n!}\int d^4x \iint\ldots \int d\xi_1d\xi_2\ldots d\xi_n \frac{\delta^n V(\Phi_c)}{\delta\Phi(\xi_1)\delta\Phi(\xi_2)\ldots\delta\Phi(\xi_n)}\Delta(\xi_1)\Delta(\xi_2)\ldots\Delta(\xi_n)$$
The $\hbar$ terms act like perturbations to a gaussian integral.
A: Briefly, 


*

*Identify $\lambda\equiv 1/\hbar$, 

*Divide the path integral measure ${\cal D}\Phi$ with the normalization factor $N(n,\hbar)$. More explicitly, this means replace ${\cal D}\Phi\to {\cal D}\frac{\Phi}{\sqrt{\hbar}}$.   

*Rescale fluctuations $\Phi=\Phi_c+\sqrt{\hbar}\Delta$ with a factor $\sqrt{\hbar}$. 


See e.g. this related Phys.SE post.
