The optical theorem, which results from the unitarity of the $S$-matrix, relates the imaginary part of the forward scattering amplitude to the total cross section. When using this theorem in practice, one often invokes the fact that in perturbation theory, $S$-matrix elements come out purely real, unless there is a Feynman diagram contribution where an intermediate particle goes on shell. According to Peskin and Schroeder, this is "easily checked" in QED.
This is true, but I had a hard time seeing why it held in more general theories. For example, consider $\phi^n$ theory. Keeping track of only phases mod $\pi/2$, each vertex comes with a factor of $i$, and simultaneously each vertex yields $n/2$ propagators, each with a factor of $i$, so a vertex gives a factor of $i^{n/2+1}$. When $n$ is odd, vertices have to come in pairs, so different orders in perturbation theory are related by a factor of $i^{n+1}$, which is real. But when $n$ is $0 \, (\text{mod} \, 4)$, different orders in perturbation theory contribute with relative factors of $i$, so it looks like the statement is false.
The only way out that I can see is to assign a factor of $i$ to every loop integral $\int d^4 k$. If such a factor exists, then it's simple to establish the result using Euler's formula.
Indeed, this is precisely what happens in dimensional regularization, where the master formula is $$\int \frac{d^dp}{(2\pi)^d} \frac{p^{2a}}{(p^2-\Delta)^b} = i(-1)^{a-b} \frac{1}{(4\pi)^{d/2}} \frac{1}{\Delta^{b-a-d/2}} \frac{ \Gamma(a+\frac{d}{2}) \Gamma(b-a-\frac{d}{2})}{\Gamma(b) \Gamma(\frac{d}{2})}$$ and the right hand side has the factor of $i$ from the Wick rotation. On one hand, this is very strange: a real integral is being regulated to an imaginary number! But on the other hand, dimensional regularization just is strange, e.g. it sets massless integrals to zero.
What I find more disturbing is the apparent requirement that loops each contribute a factor of $i$. This doesn't appear to be true of any regularization scheme I know of besides dimensional regularization. Pauli-Villars, a Wilsonian hard cutoff, and the lattice work by modifying the integrand in the purely real loop integral on the left at high energies, and so can't possibly turn it pure imaginary. That would seem to imply that all these regularization schemes violate unitarity, and in fact violate it maximally. But I've never seen anybody say that, and moreover a lattice theory in a finite box is finite-dimensional, and in this case unitarity is trivial to establish.
What's going on here?