Physical poles in QFT scattering amplitudes?

In QFT, for instance in $\phi^3$ theory, the scattering amplitudes are said to be constrained to feature so called "physical poles" only.

Consider generalized Mandelstam variables $$s_{ij},s_{ijk},s_{ijkl},...$$

defined as

$$s_{i_1i_2,...,i_m}=\left(\sum_{j=1}^m p_{i_j}\right)^2$$ where each $p_{i_j}^\mu$ is a four momentum corresponding to kinematics of external particles in the scattering process.

In case of $\phi^3$ theory the set of physical poles for $n$ point scattering is given by the generalized Mandelstam variables with indices strictly neighboring - i.e. $s_{12},s_{2,3,4}$ or wrapping around as $s_{n-1,n,1,2,3}$. What tells us that these indices should be neighboring?

Or, more generally:

I wonder where does the information come from for how exactly physical poles are supposed to look like? How would we go about finding a constraint for the shape of physical poles in a generic amplitude in a different theory?

The scattering amplitude can be obtained in QFT from the $n$-point Green's function by the LSZ reduction formula. The $n$-point Green's function are a correlation function of a string of fields :

$$G(q_1,..,q_n)=\int dx_1...dx_n\,e^{-iq_1x_1}...e^{-iq_nx_n}\langle 0|\mathcal{T}\{A_1(x_1)...A_n(x_n)\}|0\rangle$$

Now, we can factorize the time-ordering using the theta function, giving terms like:

$$\int dx_1...dx_n\,e^{-iq_1x_1}...e^{-iq_nx_n}\langle 0|\mathcal{T}\{A_1(x_1)...A_r(x_r)\}\mathcal{T}\{A_{r+1}(x_{r+1})...A_n(x_n)\}|0\rangle \times$$ $$\times\theta(\min[x_1^0...x_r^0]-\max[x_{r+1}^0...x_n^0])$$

The idea is that very old idea in quantum mechanics when we insert the identity operator $I$ in a given representation:

$$I=\sum_{p,\sigma}|p,\,\sigma\rangle\langle p,\,\sigma|+...$$

where the first term is one-particle states projector and $...$ are multi-particle states projectors. Keeping just the one particle terms, we have:

$$\int dx_1...dx_n\,e^{-iq_1x_1}...e^{-iq_nx_n}\int d^3p\sum_\sigma\langle 0|\mathcal{T}\{A_1(x_1)...A_r(x_r)\}|p,\,\sigma\rangle\times$$ $$\times\langle p,\,\sigma|\mathcal{T}\{A_{r+1}(x_{r+1})...A_n(x_n)\}|0\rangle \,\theta(\min[x_1^0...x_r^0]-\max[x_{r+1}^0...x_n^0])$$

we get a factorization of the the correlation function into two pieces, connected just by the theta function. To get ride of this we can use translation symmetry to make:

$$\langle 0|\mathcal{T}\{A_1(x_1)...A_r(x_r)\}|p,\,\sigma\rangle=e^{ip.x_1}\langle 0|\mathcal{T}\{A_1(0)...A_r(y_r)\}|p,\,\sigma\rangle$$ $$\langle p,\,\sigma|\mathcal{T}\{A_{r+1}(x_{r+1})...A_n(x_n)\}|0\rangle=e^{-ip.x_{r+1}}\langle p,\,\sigma|\mathcal{T}\{A_{r+1}(0)...A_n(y_{n})\}|0\rangle$$

and under this new variables the theta function becomes:

$$\theta(x_1^0-x_{r+1}^0+\min[0...y_r^0]-\max[0...y_n^0])$$

using the Fourier representation:

$$\theta (\tau)=-\frac{1}{2\pi i}\int_{-\infty}^{+\infty}\frac{d\omega e^{-i\omega \tau}}{\omega + i\varepsilon}$$

now we can perform the integration over $x_1$, $x_{r+1}$ and $p$. Some delta Dirac will shows up enforcing conservation of the momentum between the two blobs and an extra delta enforcing $\omega$ being equal to the energy transferred between the blobs minus the energy $E_p$ of the one-particle state. Then, the pole $(\omega +i\varepsilon)^{-1}$ that comes from the theta function will give rise to a pole $(q^{0}-E_p+i\varepsilon)^{-1}$ where $q^0$ is the energy transferred between the blobs.

Around the pole, we can make $(q^{0}-E_p+i\varepsilon)^{-1}\rightarrow 2E_p (q^2+m^2-i\varepsilon)^{-1}$, with $\vec{q}=\vec{p}$. The term $2E_P$ is absorbed by the integrals to form a relativistic invariant measure. This is how the pole $(q^2+m^2-i\varepsilon)^{-1}$ show up.

Now let us look at the residue. After the LSZ reduction formula, the residue will be precisely the product of two new amplitudes:

$$\lim_{q^2\rightarrow-m^2}(q^2+m^2-i\varepsilon)A(q_1,...,q_n)=A(q,q_2,...,q_r)\times A(q,q_{r+2},...,q_n)$$

where $q=q_1+...+q_r=-(q_{r+1}+...+q_{n})$. This means that we have a pole whatever the amplitudes $A(q,q_2,...,q_r,p)$ and $A(q,q_{r+2},...,q_n)$ are non-zero.

For more detailed explanation and calculation see Weinberg, QFT, volume 1, chapter 10.

Turns out the answer is rather simple in $\phi^3$ theory. Considering scattering at tree level for simplicity, poles can only appear in an amplitude whenever a denominator in a propagator goes to zero. In massless theory a propagator looks like $1/p^2$ where any $p$ is determined by momentum conservation at every vertex according to the Feynman rules. $\phi^3$ only has a 3-vertex, which ensures that $p=p_{i_1}+p_{i_2}+...$ is always a set of adjacent external momenta in a corresponding Feynman diagram.