If the S-Matrix is the only observable, that rules out both generalized free fields and Wick-ordered polynomials of generalized free fields as interesting Physical models, because both result in a unit S-matrix. Neither possibility has been developed since the 1960s when these results were proved, and when the S-Matrix was ascendant as the only observable in Particle Physics.
If the S-matrix is not the only quantum field observable, which it certainly seems not to be in Condensed Matter Physics and in Quantum Optics, to name just two fields in which Wightman or correlation functions play a large part in modeling, does that encourage us to construct Wick-ordered polynomials of generalized free fields as algebraic, non-dynamical deformations of free fields? In such a construction, the point is to avoid deforming a Hamiltonian or Lagrangian evolution, instead working with deformations of the observables, requiring for example that all observables and all states will be constructed as Wick-ordered polynomials in a Wick-ordered operator-valued distribution such as $\hat\Phi(x)=:\!\!\hat\phi(x)+\lambda\hat\phi^\dagger\hat\phi(x)^2\!\!:$. The Hamiltonian is taken to be a derived quantity in this approach.
The S-Matrix seems remarkable, in that it requires us to set up a hyperplane at an initial time $t_I$ and another at a final time $t_F$, and another two hyperplanes between which the interaction will be nontrivial, at times $t_I<t_A<t_B<t_F$, all of which is not Lorentz-invariant, then take the limits $t_I\rightarrow-\infty,t_A\rightarrow-\infty,t_B\rightarrow+\infty,t_F\rightarrow+\infty$. It seems that in no other branch of Physics would we construct such an idealization and say it's the only possible way to report the experimental data. In other words, I wonder, going further, whether the S-matrix even is an observable? Particle Physics seems to be the only field in Physics that seems to think it is.
post-Answer addition. It's Greenberg who proves in JMathPhys 3, 31(1962) that Wick-ordered finite polynomials in the field have a trivial S-matrix, and I think, without looking it up, that it's also Greenberg who introduced and ruled out generalized free fields. I've resisted doing something else because Wick-ordered polynomials in a generalized free field seem to give such a lot of freedom to create models that surely they can be useful, which I know is not a good reason, but I had to feel how bad it is before I could move on. I suppose I was almost ready to give up this way of doing things before I asked this question. While Tim van Beeks' response is definitely of interest, and Matt's Answer is clearly the "right" answer, but doesn't go into the periphery of the question in the right way to help me, it was Marcel's response that particularly pushed me.
My comment to Marcel indicates the way I'm now going to take this, to functions of the field such as $\hat\Phi(x)=\tanh{(\hat\phi(x))}$. Insofar as we can say that the measured value of $\hat\phi(x)$ is almost always $\pm\infty$ in the vacuum state, because the expected value $\left<0\right|\hat\phi(x)\left|0\right>$ is zero and we can say, loosely, that the variance is $\infty$ in the vacuum state, presumably $\left<0\right|\delta(\hat\Phi(x)-\lambda)\left|0\right>$ is non-zero only for $\lambda=\pm 1$ (though with some worries about this construction). It's not clear that we even have to introduce normal-ordering. $\tanh{(\hat\phi(x))}$ is of course a bounded operator of the unbounded Wightman free field, whereas no nontrivial polynomial in the field $\hat\phi(x)$ can be a bounded operator. The particular choice of $\hat\Phi(x)=\tanh{(\hat\phi(x))}$ is clearly a particular coordinatization; if we change the coordinatization by taking a function of $\hat\Phi(x)$, we in general get a measurement operator that results in two discrete values of the field, mapping $\pm 1$ to $a,b$ respectively. If we take something different from the real Klein-Gordon field, my aim as I envisage it four hours after first imagining it is to map something like an $SU(3)$ invariant Wightman field, say, to a finite number of discrete values. If such an $SU(3)$ symmetry is unbroken, the relationships between each of the discrete values will all be the same, but if the $SU(3)$ symmetry is broken, there presumably has to be a coordinate-free way in which the relationships between the different discrete values are different. There will, I now suppose, have to be enough distinction between raising and lowering operators between different sectors of the theory and measurement operators to allow there to be a concept of an S-matrix.
This may look crazy, but I'll also put on the table here why the structures of the Feynman graphs formalism encourage me. We introduce connected Feynman graphs at different orders to calculate $n$-point connected Wightman functions. Although we typically expand the series in terms of the number of loops, we can alternatively expand the series in terms of the number of points we have introduced between the $n$ points at which we measure; the extra points ensure that there are infinite numbers of different paths between the $n$ points. With something like $\hat\Phi(x)$, we introduce an infinite number of paths directly between the $n$ points, without introducing any extra points, so we need, in effect, a transformation of the superposition of an infinite number of Feynman path integrals into a superposition of an infinite number of weighted direct paths between the $n$ points. Getting the weights on the direct paths right is of course rather important, and I also imagine the analytic structure has to be rather carefully done, particularly if I don't use normal-ordering. I'll do any amount of work to avoid renormalization, even in its modern gussied up form.
It seems to me significant that the quantum field $\hat\Phi(x)=\tanh{(\hat\phi(x))}$ is not reducible. A number of proofs concerning Wightman fields rely on this property.
If anyone else understands this (or reads this far) I'll be surprised. In any case I expect it will look very different a few years down the road if I ever manage to get it into a journal. Although I've worked in and around the Wightman axioms for the last few years, I find it interesting that I can now feel some pull towards something like the Haag-Kastler axioms. Lots of work to do! Thank you all! Good luck with your own crazy schemes!
That's the question (and the to me unexpected state of play a day later). Completely separately, as an example, to show the way I'm going with this, hopefully (which, a day later, looks as if it will be only a background concern for the next little while, but I think not likely to be completely forgotten by me), the real-space commutator of the creation and annihilation operator-valued distributions of the free field of mass $m$ is, in terms of Bessel functions, $$C_m(x)=\frac{m\theta(x^2)}{8\pi\sqrt{x^2}}\left[Y_1(m\sqrt{x^2})+i\varepsilon(x_0)J_1(m\sqrt{x^2})\right]$$ $$\qquad\qquad+\frac{m\theta(-x^2)}{4\pi^2\sqrt{-x^2}}K_1(m\sqrt{-x^2})-\frac{i}{4\pi}\varepsilon(x_0)\delta(x^2).$$ If we take a weighted average of this object with the normalized weight function $w_{\alpha,R}(m)=\theta(m)\frac{R(Rm)^{\alpha-1} {\rm e}^{-Rm}}{\Gamma(\alpha)}, 0<\alpha\in \mathbb{R},\ 0<R$, $\int w_{\alpha,R}(m)C_m(x)dm$, we obtain the commutator of a particular generalized free field, which can be computed exactly in terms of Hypergeometric functions and which at space-like separation $\mathsf{r}=\sqrt{-x^\mu x_\mu}$ is asymptotically $\frac{\Gamma(\alpha+1)R^{\alpha}}{\pi\Gamma(\frac{\alpha+1}{2})^2\bigl(2\mathsf{r})^{\alpha+2}}$, and at time-like separation $\mathsf{t}=\sqrt{x^\mu x_\mu}$ is asymptotically $$\frac{\cos{(\frac{\pi\alpha}{2})}\Gamma(\frac{\alpha}{2}+1)R^\alpha} {4\sqrt{\pi^3}\Gamma(\frac{\alpha+1}{2})\mathsf{t}^{\alpha+2}} -i\frac{R^\alpha}{4\Gamma\left(-\frac{\alpha}{2}\right)\Gamma\left(\frac{\alpha+1}{2}\right)\mathsf{t}^{\alpha+2}},$$ except for $\alpha$ an even integer. At small space-like or time-like separation, the real part of this generalized free field is $-\frac{1}{4\pi^2x^\mu x_\mu}$, identical to that of the massless or massive free particle, independent of mass, but we can tune the 2-point function at large distances to be any power of the separation smaller than an inverse square law. On the light-cone itself, the delta-function component is again identical to that of the massless or massive free particle, independent of mass.
There is of course an infinity of possible normalized weight functions, a half-dozen of which I have worked out exactly and asymptotically, and somewhat obsessively, by use of MAPLE and Gradshteyn & Ryzhik, though I've managed to stop myself at the moment. In a subsequent edit, I can't resist adding what we obtain if we use the weight function $w_{\mathsf{sm}[R]}(m)=\frac{\theta(m)R\exp{\left(-\frac{1}{mR}\right)}}{2(mR)^4},\ 0<R$ [using 6.591.1-3 from Gradshteyn & Ryzhik]. This function is smooth at $m=0$, and results at time-like separation in $$\left\{Y_2\!\left(\sqrt{\frac{2\mathsf{t}}{R}}\right)+iJ_2\!\left(\sqrt{\frac{2\mathsf{t}}{R}}\right)\right\} \frac{K_2\!\left(\sqrt{\frac{2\mathsf{t}}{R}}\right)}{8\pi R^2} \asymp \frac{\exp{\left(-i\sqrt{\frac{2\mathsf{t}}{R}}+\frac{\pi}{4}\right)} \exp{\left(-\sqrt{\frac{2\mathsf{t}}{R}}\right)}} {16\pi\sqrt{R^3\mathsf{t}}}$$ and at space-like separation in $$\frac{K_2\!\left((1+i)\sqrt{\frac{\mathsf{r}}{R}}\right) K_2\!\left((1-i)\sqrt{\frac{\mathsf{r}}{R}}\right)}{4\pi^2R^2} \asymp \frac{\exp{\left(-2\sqrt{\frac{\mathsf{r}}{R}}\right)}} {16\pi\sqrt{R^3\mathsf{r}}}.$$ I wish I could do this integral for more general parameters, but hey! This weight function has observable effects only at space-like and time-like separation $\mathsf{r}<R$ and $\mathsf{t}<R$, so this is essentially unobservable if $R$ is small enough. The kicker is that this decreases faster than polynomially in both space-like and time-like directions.
The generalized free field construction always results in a trivial $n$-point function for $n>2$, however by introducing also Wick-ordered polynomials of these generalized free fields, we can also tune the $n$-point connected correlation functions, which are generally non-trivial and to my knowledge finite for all $n$. All this is far too constructive, of course, to prove much. I think I propose this more as a way to report the $n$-point correlation functions in a manner comparable to the Kallen-Lehmann representation of the 2-point correlation function than as something truly fundamental, because I think it does not generalize well to curved space-time.