Why is the argument of the path integral a pure phase? [Edit: moved to front] For which $(H,b)$ pairs, where $H$ is a Hamiltonian and $b$ is a basis, can we write:
$$\langle b_f\vert e^{-iHt/\hbar}\vert b_i\rangle=\int_{b(0)=b_i}^{b(t)=b_f} \mathcal{D}b(t)\, e^{iS[b(t)]},\tag{2}$$
where $S$ is real?

Background
The derivation of the path integral notion I am familiar with proceeds as follows. We would like to evaluate the propagator:
$$G(m,n,t)=\langle n\vert U(t)\vert m\rangle. $$
For some large number $N$, we break up the time evolution from $0$ to $t$ into $N$ pieces. For notational clarity, define:
$$W:=U\left(\frac{t}{N}\right),$$
so that:
$$G(m,n,t)= \langle n\vert W^N\vert m\rangle.$$
Now, denoting the matrix elements of $W$ by $W_{ij}$, we may write:
$$G(m,n,t)=\sum_{j,k,\ell,\dots}\underbrace{W_{nj}W_{jk}W_{k\ell}\cdots W_{\ell m}}_{\text{N times}}.$$
Thus we have written $G$ by summing over "histories" $\{j,k,\ell,\cdots\}$, with a product of $W$'s matrix elements being assigned to each history. [edited to add following Qmechanic's answer:] Assuming the $W_{ij}$'s are all nonzero, we may define $s_{ij}=-i\log W_{ij}$, and write:
$$G(m,n,t)=\sum_{i(t)\in\{\text{hists}\}}\exp\left(i\sum_{t}s_{i(t+1),i(t)}\right),$$
where now $s_{ij}$ has the interpretation of the increment of action. Note that nothing in our analysis so far has implied that $s$ is real.
For the particular case of working in the position basis with a Hamiltonian of the form $$H=p^2/(2m)+V(q),$$ a gaussian integral may be used to show that (as $N\to \infty$):
$$W_{ij}=\mathcal{N}e^{iF_{ij}},\tag{1}$$ 
with $F_{ij}$ real, which means that if we absorb the normalization into the definition of the path integral, the $s$'s are real.
This is kind of bizarre. To emphasize, what we have found is not that $W$ is unitary or any nice condition like that. We have found that all of the entries $W_{ij}$ have the same absolute value. I would not expect this property to be preserved under a unitary change of basis.
It is unobvious to me how general this result is (my suspicion is that it will definitely not hold for "most" Hamiltonians/bases). For which Hamiltonians/choice of bases are these matrix elements pure phases, and is there a simple or physical way to see that this must be so?

Counterexamples
I do not believe that (1) and (2) need to be satisfied in any basis for any choice of $H$. I will try to argue this by presenting 3 examples, in decreasing order of my confidence in them.
First consider $H=0$. In this case, $W_{x,x'}$ is nonzero when $x=x'$, and zero when $x\not =x'$. This shows that (1) cannot be satisfied, and indeed the $s$'s may not even be defined. I thus believe no $S$ exists to make (2) true.
We can also write down something pathological like:
$$H=\lambda\vert x=1\rangle\langle x=0\vert +\text{h.c.},$$
which, to my inexperienced eye, looks like $W_{0,0}$, $W_{0,1}$, and $W_{2,3}$ will all have different magnitudes. Again (1) cannot be satisfied.
Alternatively, we may consider a normal-looking Hamiltonian:
$$H=\frac{p^2}{2m}+E_0\log\left(1+\frac{x^2}{x_0^2}\right),$$
but consider the matrix elements of $W$ in the momentum basis. Doing the same Trotter formula trick and evaluating the integral over $x$ this time gives a result which I believe is not of the form $$W_{p,p'}=\mathcal{N}e^{iF(p,p')}$$ with $F$ real. I believe this because plugging in different values of $p$ and $p'$ numerically changes $\vert W_{p,p'}\vert ^2$.
Generically, for a separable local Hamiltonian $H(q,p)=K(p)+V(q)$, a phase-space path integral may be defined, and the position space path integral may be recovered by integrating out $p(t)$. I expect that only in the case of $K$ quadratic do we recover a position space path integral whose argument is of the form $e^{iS}$ for $S\in\mathbb{R}$, and generically the result would have a varying absolute value. I believe this provides some evidence that "representable as a real-action path integral wrt a particular basis" is a fairly strong constraint on Hamiltonians.

 A: That the Boltzmann factor $\exp\left(\frac{i}{\hbar}S\right)$ is a pure phase follows essentially from the fact that the time-evolution operator $U=\exp\left(-\frac{i}{\hbar}H\right)$ is unitary via the general derivation of the path integral formalism from the operator formalism. See also this related Phys.SE post.
Comments to update:


*

*There seems to be no reason why OP's eq. (1) should hold.

*For OP's potential $V=E_0{\rm Ln}\left(1+\left(\frac{x}{x_0}\right)^2\right)$, the Mathematica calculation seems to be calculating the momentum Fourier transform of $\exp\left(\frac{i}{\hbar}V\right)$. There seems to be no reason why that should be a pure phase.
A: Above equation (6.9), Srednicki's QFT book says

If $H(p,q)$ is no more than quadratic in the momenta ... then the integral over $p$ is gaussian, and can be done in closed form. If the term that is quadratic in $p$ is independent of q ... then the prefactors generated by the gaussian integrals are all constants, and can be absorbed into the definition of $\mathcal{D}q$. The result of
  integrating out $p$ is then
$$ \langle q'',t''| q', t' \rangle = \int \mathcal{D}q \exp \left[ i \int_{t'}^{t''} dt\, \mathcal{L}(\dot{q}(t), q(t)) \right],$$

where $\mathcal{L}$ is calculated from $H$ via the usual classical procedure.  So basically, the Hamiltonian has to take the form $H = p^2/(2m) + p\, A(x) + V(x)$. (The second term can come up as a cross-term between the canonical momentum and a vector potential in gauge QFT.)
The key point in understanding where those phases in the matrix elements come from is this: if $x$ and $p$ are canonically congugate variables (which they are), then the Stone–von Neumann theorem roughly says that they can be represented as Fourier transforms of each other, so that $\langle x | p \rangle \propto e^{i p x/\hbar}$.  That's why all the matrix elements have the same magnitude: in any momentum eigenstate of $p$, the probability $P(x) = |\psi(x)|^2 = | \langle x | p \rangle|^2$ of finding a particle at any position should be the same, and vice versa for position eigenstates.
