So I am looking at a statistical path integral, meaning that I work with an Euclidean action. The propagator of my (Wiener) path integral is given by: $$ K(x_T,T|x_0,0)=\int\limits_{x(0)=0}^{x(T)=x_T}\mathcal{D}x\exp\left(-\int\limits_0^T\left[\frac{m}{2}\left(\dot{x}\right)^2+fx\right]d{t}\right), $$ which is basically a free particle in a gravity potential. Since the action is quadratic, the WKB formula $$ K(x_T,T|x_0,0)\approx\sqrt{-\frac{1}{2\pi}\frac{\partial^2 S[x_\mathrm{kl}(t)]}{\partial x_0\partial x_T}}\exp\left(-S[x_\mathrm{kl}(t)]\right) $$ should be exact.
The equation of motion gives me that the path of the particle meeting the boundary conditions is given by: $$ x_\mathrm{kl}(t)=\frac{f}{2m}(t-T)t+\frac{x_T- x_0}{T}t + x_0. $$ Using this path, I can calculate the classical action which becomes equal to: $$ S_\mathrm{kl}=-\frac{f^2}{24m}T^3+\frac{f T}{2}(x_T+x_0)+\frac{m}{2}\frac{(x_T-x_0)^2}{T}. $$
Substituting all of my results in the WKB formula then yields that the propagator is now given by: $$ K(x_T,T|x_0,0)=\sqrt{\frac{m}{2\pi T}}\exp\left(-\frac{m}{2}\frac{(x_T-x_0)^2}{T}-\frac{f T}{2}(x_T+x_0)+\frac{f^2T^3}{24m}\right). $$ The problem however with this propagator is that it does not stay normalized. If I demand that the propagator should remain normalized at all times $T$, then my propagator is given by: $$ K(x_T,T|x_0,0)=\sqrt{\frac{m}{2\pi T}}\exp\left(-\frac{m}{2}\frac{(x_T-x_0)^2}{T}-\frac{f T}{2}(x_T+x_0)+\frac{f^2T^3}{24m}\color{blue}-\color{blue}{\frac{fT}{6}}\left[\color{blue}{\frac{fT^2}{m}-6x_0}\right]\right), $$ which yields an extra term compared to the first version.
Extra: Time sliced method
I think that I have found the source of my problem, and it can be seen by looking at the infinitesimal propagator given by $$ K(x_j,t_j|x_{j-1},t_{j-1})=\sqrt{\frac{m}{2\pi\Delta t_j}}\exp\left(x_{j-1} f \Delta t_j + \frac{f^2}{2m}(\Delta t_j)^3\right)\\\times\exp\left(-\frac{m}{2\Delta t_j}\left[x_j-\left(x_{j-1}+\frac{f}{m}(\Delta t)^2\right)\right]^2\right). $$ In the upper part we indeed see that the normalization gets an extra exponential factor which will cause the path integral propagator (in time) to diverge. Also note that also has (more or less) the same form as the needed normalization factor (supporting my claim above)!
Question (new): For computational simplicity I'd like my propagator to still stay normalized. Is it okay to just use the second normalized version for my expectation values, or is that just wrong? Answer to old question of course still welcome as it may become relevant for further exploration of the path integral.
Question (old): Is there some kind of extra theorem that puts limitations on the correctness of the WKB formulas, or did I miss an extra important term here? I have recalculated the result a couple of times and it all seems to be correct at first sight.
About the solution
I also checked my solution in literature (Dittrich & Reuter) but they found the same (divergent) solution without any explanation. So at least I know the found solution is correct. Unfortunately I still have no idea what this means for my physics.