Products of distributions in QFT In Quantum Field Theory quantum fields are operator valued distributions. Namely, given the Schwartz space $\mathcal{S}(M)$ defined on Minkowski spacetime $M$, fields are continuous linear maps $\phi : \mathcal{S}(M)\to \mathcal{L}(\mathcal{H})$ which gives back for each function one operator in some Hilbert space.
When just free fields are considered this is just fine. Now, when one considers interactions, the issue that appears is that we have to deal with products of quantum fields.
For instance, the $\phi^4$ theory has the interaction given by $\lambda \phi^4$. The scalar Yukawa theory also has products of fields.
The field equations become non-linear in the fields. The issue now is that if quantum fields are distributions, then necessarily we can't deal with these terms, because products of distributions is not defined. Physicists just ignore this most of the time.
My question here is: is there a rigorous known way to deal with this? Can we make this rigorous someway, even if we need to go into approximate theories? Or there is no way currently known to make this rigorous?
Edit: For more information on this matter, I've asked on Physics.SE about the issues with interacting QFT and found out the main problems. This product problem seemed the worst, and since the point here is about whether or not a mathematical construct can be made rigorous I thought that Math.SE could be a good place to discuss the matter further.
 A: I'm definitely not a QFT expert, but I believe this is what Martin Hairer has been working on.  In particular, his Regularity Structures allow for a rigorous treatment of products of distributions that arise in QFT and other nonlinear stochastic PDE scenarios.  He gives a fairly easy to follow talk on the subject here, where I seem to recall he discusses $\Phi^4$ renormalization theory.  I don't recall if he addresses operator valued distributions, but it doesn't seem too difficult to believe that if you can work out the theory for scalar valued distributions, operators are not too far off via spectral theory.
A: The fields in Quantum Field Theory (QFT) are taken (in Wightman axiomatic framework) to be operator-valued distributions. I try to answer your question step by step.


*

*What is the best representation of a field?


the main idea of QFT is that the degree of freedom are infinitely large. For example, many harmonic oscillators (of different momentum) can be obtained by action of creation operator on vacuum. The suitable representation here is the  Fourier transform of free field. Similarly, in case of interacting QFT, one has to define some operators playing the role of creation and annihilation operators to include basic idea of QFT. 


*How to multiply such operators?


We multiply them without taking care of their distribution nature.
3.If so, why do we need operators to be distributions?
The distribution nature of fields comes to play when you do renormalization. As Dirac insisted, renormalization is nothing but subtracting infinities (that are denoted by cut-off in integrals).
The infinities appear not only in interacting QFT (and renormalization), but also in free field theory.  For example in the calculation of Hamiltonian operator of free field, an infinity appears and everyone is satisfied by saying that the difference of energy is physical. Here distribution comes to provide a reasonable foundation for such ill-defined processes!!


*Where does multiplication of distributions appear?


Perhaps the most serious situation where multiplication of distributions comes to play is the multiplication of two (or more) Dirac delta functions that eventually takes the form (in perturbative QFT) 
$$\int_{\mathbb{R}} \delta(x-y)\delta(x-z)dx = \delta(y-z)$$ 
That can be easily verified by taking $\delta(x-y)=\frac{d}{dx}H(x-y)$ where $H(x)$ is Heaviside step function and then integrating by part.  
A: In addition to the other useful answers: when two distributions have disjoint "singular support", they can be multiplied meaningfully. More delicately, when they have disjoint "wavefront sets", they can be multiplied meaningfully. This seems not entirely adequate for the apparent needs of QFT, though I am very inexpert about the latter, so may be wrong. In general, I've found Dirac's intuitions very good, even if non-trivial to parse in "strict" mathematics.
A: There is a formal approach to algebraic quantum field theory based on Epstein-Glaser causal perturbation theory combined with tools from microlocal analysis and the Steinmann scaling degree method. You may find the paper 'Causal Perturbation Theory in Terms of Retarded Products, and a Proof of the Action Ward Identity' by Dütsch and Fredenhagen  interesting, along with references therein.
The idea is to restrict the products by certain axioms, thereby allowing the inductive construction of higher order products from lower order products. In the construction, higher order products $\phi(x_1)\cdots\phi(x_n)$ are uniquely determined from lower order products, except on the 'total diagonal' $x_1=\ldots=x_n$. The extension across the total diagonal can be carried out by considering the Steinmann scaling degree, but this extension introduces some amount of non-uniqueness (corresponding to subtraction terms in more standard approaches).
The procedure is formal, because the products one obtains are formal power series in $\hbar$, rather than proper convergent power series. In this sense, the approach is mathematically completely rigorous, although it leaves out an interesting component, namely the question of convergence of the power series.
Dütsch, Michael; Fredenhagen, Klaus, Causal perturbation theory in terms of retarded products, and a proof of the action ward identity, Rev. Math. Phys. 16, No. 10, 1291-1348 (2004). ZBL1084.81054.
A: There are rigorous ways to deal with this in some cases. As was mentioned, Hairer's work is probably related, although I don't know much about it. My understanding is that polynomial interaction terms are ok in 2 dimensions (you have to regularize the distributions e.g. by convolution and then do some analysis), and some polynomial interactions work in 3 dimensions too. Part II of the Glimm & Jaffe book deals with this (in 2 dimensions), although the book is quite advanced. See also this review by Jaffe.
More recently there has been work on an exponential interaction $e^\phi$, see Liouville field theory, in two dimensions. This is an exponential of a distribution and it was defined by using the theory of Gaussian multiplicative chaos, developed initially by Jean-Pierre Kahane in the 1980's.
Related to all these things is this review. All in all the subject is very delicate and it is hard to give any of the details briefly.
A: It can be done rigorously using as input the operator product expansion for pointwise correlation with suitable remainder bounds. This is the object of my CMP article "A second-quantized Kolmogorov-Chentsov theorem via the operator product expansion" (see also here for the preprint version). A quick account can be found in the Oberwolfach report for this 2016 workshop.
