Statement of exercise
On a page 11 of A.Zee's book QFT in a Nutshell, he derives Dirac's formulation of the path integral formulation of QM for a free particle. This starts with the free particle Hamiltonian
$$ H = \frac{\hat{p}^2}{2m} $$
The hat on $\hat{p}$ says that it is meant as an operator. Then he divides the time interval into $N$ parts of length $\delta t$, where each part is
$$ \begin{align} \tag{1}\langle q_{j+1}|e^{-iH \delta t}|q_j\rangle & = \langle q_{j+1}|e^{-i\delta t\,\hat{p}^2 /(2m)}|q_j\rangle \\ \tag{2} & = \int\frac{dp}{2\pi}\langle q_{j+1}|e^{-i\delta t\,\hat{p}^2 /(2m)}|p\rangle\langle p|q_j\rangle \\ \tag{3} & = \int\frac{dp}{2\pi}e^{-i\delta t\,p^2 /(2m)} \langle q_{j+1}|p\rangle\langle p|q_j\rangle \end{align} $$
In $(2)$, he used $1 = \int(dp/2\pi)|p\rangle\langle p|$ with momentum eigenstates $|p\rangle$, and in $(3)$ the fact that $\hat{p}|p\rangle = p|p\rangle$.
Then he goes on to derive the fact that the time integral of the Lagrangian emerges as the exponent in
$$ \langle q_F|e^{-iHt}|q_I\rangle = \int Dq\,e^{i\int dt S} $$
Now Zee poses the exercise to derive the same thing using the Hamiltonian
$$ H = \frac{\hat{p}^2}{2m}+V(\hat{q}) $$
What I tried
I was able to derive the result, but only IF the following is true
$$ \tag{4}\left\langle q_{j+1}\middle|e^{-i\delta t\left\{\hat{p}^2/(2m)+V(\hat{q})\right\}}\middle|q_j\right\rangle = \left\langle q_{j+1}\middle|e^{-i\delta t\left\{\hat{p}^2/(2m)+V(q_j)\right\}}\middle|q_j\right\rangle $$
in other words, if I can replace $V(\hat{q})$ with $V(q_j)$, then I can prove the result, because then, I can repeat the free particle calculation and the potential appears as an extra term in the integral that depends only on $q_j$.
However, I don't know how to prove that $(4)$ is true. $\hat{p}$ and $\hat{q}$ do not commute, so $e^{\hat{p}+V(\hat{q})}$ cannot simply be factorized into $e^{\hat{p}}e^{V(\hat{q})}$.
Why do I think that $(4)$ is true?
The following approximation led me to think that $(4)$ is true:
$$ \begin{align} \tag{5}\left\langle q_{j+1}\middle|e^{-it\left\{\hat{p}^2/(2m)+V(\hat{q})\right\}}\middle|q_j\right\rangle & \approx \langle q_{j+1}|1-it\,(\hat{p}^2/2m+V(\hat{q}))+O(\delta t^2)|q_j\rangle \\ & = \langle q_{j+1}|1-it\,\hat{p}^2/2m+O(\delta t^2)|q_j\rangle - i\delta t\langle q_{j+1}|V(\hat{q})|q_j\rangle \\ \tag{7}& = \langle q_{j+1}|1-it\,\hat{p}^2/2m+O(\delta t^2)|q_j\rangle - i\delta t\langle q_{j+1}|V(q_j)|q_j\rangle \\ & = \langle q_{j+1}|1-it\,(\hat{p}^2/2m+V(q_j))+O(\delta t^2)|q_j\rangle \\ \end{align} \\ $$
$(5)$ should be possible since we're only interested in the limit $\delta t \to 0$. But this seems to be a rather sloppy argument, because in the following steps, before we take the limit $\delta t \to 0$, a Gaussian integral over $p$ is used. So we really need the exponential function and not just some approximation.
Even $(7)$ seems a bit sloppy. Shouldn't I have some limits on the $V(q)$ function to know that $V(\hat{q})|q_j\rangle = V(q_j)|q_j\rangle$?
Any hints?
