Let $\phi(x)$ be a free Hermitian scalar field in $4D$ Minkowski spacetime with the metric $(1,-1,-1,-1)$.

Then, though I wrote it as $\phi(x)$, it is in fact an operator-valued tempered distribution on a Bosonic Fock space, requiring a test function to be smeared.

Moreover, there is a decompositon $\phi(x)=\phi^+(x)+\phi^-(x)$ into the creation part $\phi^+(x)$ and the annihilation part $\phi^-(x)$.

Let $D(x-y)=[\phi^-(x),\phi^+(y)]$ where $[,]$ is the commutator so that $D(x-y)$ is a version of the propagator. It is known that $D(x-y)$ is a numerical valued tempered distribution.

Then, the normal ordering is defined in all literatures as follows: \begin{equation} :\phi^2:(x)=\lim\limits_{y \to x} \Bigl(\phi(x)\phi(y)+iD(x-y)\Bigr) \end{equation}

However, I have great difficulty understanding the above formula rigorously, since every term is a distribution. How should I make sense of the limit $\lim\limits_{y \to x}$ in this case?

Could anyone please clarify?