How to prove that the quantity appearing in the exponent of the path integral is the Lagrangian? In Zee's Quantum Field Theory in a Nutshell, it is shown that, if $H = \frac{\hat{p}^2}{2m}$, then
\begin{equation}
\langle q_F | e^{-iHt} | q_I \rangle = \int e^{i\int \frac{1}{2}m\dot{q}^2 \, dt} Dq
\end{equation}
and there is a remark that if $H = \frac{\hat{p}^2}{2m} + V(q)$ then
\begin{equation}
\langle q_F | e^{-iHt} | q_I \rangle = \int e^{i\int \frac{1}{2}m\dot{q}^2 - V(q) \, dt} Dq
\end{equation}
This can be handily proven once we realize we are allowed to approximate $e^{(A + B)\delta t}$ as $e^{A \delta t} e^{B \delta t}$ in the limit where $\delta t \to 0$.
and Zee concludes that the quantity appearing in the integral in the exponent is just the Lagrangian.
I'm having trouble seeing how to derive this for more complex Hamiltonians (or in general). For example, let's generalize to 3 dimensions and take
\begin{equation}
H = \frac{1}{2m} (\hat{\mathbf{p}} - e\mathbf{A}(\hat{\mathbf{x}}))^2 + e \varphi(\hat{\mathbf{x}})
\end{equation}
If we follow Zee's derivation, we must then figure out how to evaluate the amplitude to propagate from $\mathbf{x}_j$ to $\mathbf{x}_{j+1}$ in time $\delta t$,
\begin{equation}
\left\langle \mathbf{x}_{j+1} \left| \ \exp\left(-i \, \delta t \, \left[\frac{1}{2m}(\hat{\mathbf{p}}^2 - e\mathbf{A}(\hat{\mathbf{x}}) \cdot \hat{\mathbf{p}} - e\hat{\mathbf{p}} \cdot \mathbf{A}(\hat{\mathbf{x}}) + e^2 \hat{\mathbf{A}}^2(\hat{\mathbf{x}})) + e \varphi(\hat{\mathbf{x}})\right]\right) \ \right| \mathbf{x}_j \right\rangle = \, ?
\end{equation}
Now I'm stuck. We can again consider the terms appearing in the exponent separately, but I have no idea how to evaluate something of the form $\langle \mathbf{x}_{j+1} | \exp(ik \hat{\mathbf{p}} \cdot \mathbf{A}(\hat{\mathbf{x}})) | \mathbf{x}_j\rangle$.
I know that in general multiplying the momentum operator by a constant then taking the exponential produces a translation operator, but I don't think that's valid here because instead of a constant we have $A(\hat{x})$,  an operator that doesn't commute with $\hat{p}$. Nevertheless, it turns out that if we pretend that this works, we end up getting the right answer. But this must be luck... right?
More generally, it's hard to see how evaluating these amplitudes can, in general, perform a Legendre transformation on the Hamiltonian. Why should this procedure be able to extract the needed derivative of $H$ with respect to $p$?
 A: *

*The derivation of the naive formal Lagrangian path integral
$$\begin{align}\langle q_f,t_f|q_i,t_i \rangle~\sim~& \int_{q(t_i)=q_i}^{q(t_f)=q_f} \!{\cal D}q \exp\left[\frac{i}{\hbar}S[q]  \right], \cr S[q]~=~&\int_{t_i}^{t_f}\!dt~L(q),\end{align} \tag{1}$$
and the naive formal Hamiltonian phase space path integral
$$\begin{align}\langle q_f,t_f|q_i,t_i \rangle~\sim~&\int_{q(t_i)=q_i}^{q(t_f)=q_f} \!{\cal D}q~{\cal D}p~\exp\left[ \frac{i}{\hbar}S_{H}[q,p]\right],\cr  S_H[q,p]~=~&\int_{t_i}^{t_f}\!dt~\left[ p\dot{q}- H(q,p)\right],\end{align}\tag{2}$$
is done in many textbooks. The trick is to insert infinitely many completeness relations of position (& momentum) eigenstates.


*The words naive & formal are here used to stress the fact that most derivations in textbooks do not discuss operator ordering ambiguities, and this seems to be OP's actual question. This is a huge topic in itself. See e.g. this Phys.SE post and links & references therein.


*Let us finally mention that in OP's concrete example from E&M, much progress can be done by working with the Hamiltonian phase space path integral (2) rather than the Lagrangian path integral (1), since position and momentum operators can then be projected onto their respective eigenstates (possibly after adequate operator commutations).
