Notice that the above procedure is not the same as pulling back the matrix elements $\langle\Omega_1,\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)\Omega_2\rangle$$\langle\Omega_1,:\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)\Omega_2\rangle$ of the auxiliary Wick monomial $:\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)$ between two vectors $\Omega_1,\Omega_2\in\mathscr{D}$ by $\iota_n$ as suggested by the OP in the case $n=2$. The difference is the following: one takes the pullback by $\iota_n$ in the way described above, by first mollifying the auxiliary Wick monomials with $\delta_\epsilon$ before smearing the arguments of $\Omega_1,\Omega_2$ with test functions. Only then can the necessary singularity cancellations take place as the limit $\epsilon\rightarrow 0$ is performed at the end. More precisely, this procedure guarantees that one never makes both arguments of a single factor $\omega_2$ in the Wick formula coincide in that limit, which (as argued by 1. above) is something that cannot be done. The way proposed by the OP falls into that very trap - in other words, restricting $G$ to $\Delta_2$ (more precisely, pulling back $G$ by $\iota_2$) is not a well defined procedure by the same reason as restricting $F$ to $\Delta_2$ is not.
Notice that the above procedure is not the same as pulling back the matrix elements $\langle\Omega_1,\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)\Omega_2\rangle$ of the auxiliary Wick monomial $:\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)$ between two vectors $\Omega_1,\Omega_2\in\mathscr{D}$ by $\iota_n$ as suggested by the OP in the case $n=2$. The difference is the following: one takes the pullback by $\iota_n$ in the way described above, by first mollifying the auxiliary Wick monomials with $\delta_\epsilon$ before smearing the arguments of $\Omega_1,\Omega_2$ with test functions. Only then can the necessary singularity cancellations take place as the limit $\epsilon\rightarrow 0$ is performed at the end. More precisely, this procedure guarantees that one never makes both arguments of a single factor $\omega_2$ in the Wick formula coincide in that limit, which (as argued by 1. above) is something that cannot be done. The way proposed by the OP falls into that very trap - in other words, restricting $G$ to $\Delta_2$ (more precisely, pulling back $G$ by $\iota_2$) is not a well defined procedure by the same reason as restricting $F$ to $\Delta_2$ is not.
Notice that the above procedure is not the same as pulling back the matrix elements $\langle\Omega_1,:\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)\Omega_2\rangle$ of the auxiliary Wick monomial $:\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)$ between two vectors $\Omega_1,\Omega_2\in\mathscr{D}$ by $\iota_n$ as suggested by the OP in the case $n=2$. The difference is the following: one takes the pullback by $\iota_n$ in the way described above, by first mollifying the auxiliary Wick monomials with $\delta_\epsilon$ before smearing the arguments of $\Omega_1,\Omega_2$ with test functions. Only then can the necessary singularity cancellations take place as the limit $\epsilon\rightarrow 0$ is performed at the end. More precisely, this procedure guarantees that one never makes both arguments of a single factor $\omega_2$ in the Wick formula coincide in that limit, which (as argued by 1. above) is something that cannot be done. The way proposed by the OP falls into that very trap - in other words, restricting $G$ to $\Delta_2$ (more precisely, pulling back $G$ by $\iota_2$) is not a well defined procedure by the same reason as restricting $F$ to $\Delta_2$ is not.
(Edit in response to updates in OP) I will provide some of the details below as guidance to the reader. Let $\mathscr{H}$ be the Klein-Gordon Fock (Hilbert) space, which is the completion of the linear span $\mathscr{D}$ of vectors of the form $\phi(f_1)\cdots\phi(f_m)\Omega$, where $f_1,\ldots,f_m$ are test functions on $\mathbb{R}^4$ and $m\in\mathbb{N}$ and $\Omega$ is the vaccum state vector, just like in 1. above. In other words, $\mathscr{D}$ is the common domain of definition of the smeared Klein-Gordon field operators $\phi(f)$. One first defines the so-called auxiliary Wick monomials recursively as the following sequence of distributions $:\!\phi^{\otimes n}\!:$$:\!\phi^{\otimes n}\!\!:$ in $\mathbb{R}^{4n}$ taking values in unbounded linear operators on $\mathscr{H}$ with common domain $\mathscr{D}$:
$$:\!\phi^{\otimes 0}\!:\,\,=\,\mathbb{1}\ ,\,\,:\!\phi^{\otimes 1}\!:\!\!(x)=\phi(x)\ ;$$$$:\!\phi^{\otimes 0}\!\!:\,\,=\,\mathbb{1}\ ,\,\,:\!\phi^{\otimes 1}\!\!:\!(x)=\phi(x)\ ;$$ $$:\!\phi^{\otimes(n+1)}\!:\!\!(x_1,\ldots,x_{n+1})=\,\,:\!\phi^{\otimes n}\!:\!\!(x_1,\ldots,x_n)\phi(x_{n+1})$$$$:\!\phi^{\otimes(n+1)}\!\!:\!(x_1,\ldots,x_{n+1})=\,\,:\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)\phi(x_{n+1})$$ $$\phantom{\!\phi^{\otimes(n+1)}\!:\!\!(x_1,\ldots,x_{n+1})=\,\,}-\sum^n_{k=1}\omega_2(x_k,x_{n+1})\,\!\!:\!\phi^{\otimes(n-1)}\!:\!\!(x_1,\ldots,\widehat{x_k},\ldots,x_n)\ ,\,n\in\mathbb{N}\ ,$$$$\phantom{\!\phi^{\otimes(n+1)}\!\!:\!(x_1,\ldots,x_{n+1})=\,\,}-\sum^n_{k=1}\omega_2(x_k,x_{n+1})\,\!\!:\!\phi^{\otimes(n-1)}\!\!:\!(x_1,\ldots,\widehat{x_k},\ldots,x_n)\ ,\,n\in\mathbb{N}\ ,$$
where the hat over an argument denotes its omission. This is exactly the standard formula for the Wick monomials $:\!\phi^n\!:\!\!(x)$ prior to the restriction to the small diagonal $\Delta_n=\{(x,\ldots,x)\ |\ x\in\mathbb{R}^4\}\subset\mathbb{R}^{4n}$ (so that $\Delta_2=\Delta$). To perform the latter operation, we define the mollified Wick monomials $:\!\phi^{\otimes n}_\epsilon\!:$ as $$:\!\phi^{\otimes n}_\epsilon\!:(f)=:\!\phi^{\otimes n}\!:(\delta_\epsilon(f))\ ,\,n>1\ ,f\in\mathscr{D}(\mathbb{R}^4)\ ,$$$$:\!\phi^{\otimes n}_\epsilon\!\!:\!(f)=:\!\phi^{\otimes n}\!\!:\!(\delta_\epsilon(f))\ ,\,n>1\ ,f\in\mathscr{D}(\mathbb{R}^4)\ ,$$ where $\delta_\epsilon:\mathscr{D}(\mathbb{R}^4)\rightarrow\mathscr{D}'(\mathbb{R}^{4n})$ is a continuous linear map taking values in $\mathscr{D}(\mathbb{R}^{4n})$ for each $\epsilon>0$ and converging in $\mathscr{D}'(\mathbb{R}^{4n})$ to the formal adjoint ( = transpose) of the pullback $\iota^*_n$ by the diagonal map $\iota_n:\mathbb{R}^4\rightarrow\mathbb{R}^{4n}$, $\iota_n(x)=(x,\ldots,x)$ (so that $\iota_n(\mathbb{R}^4)=\Delta_n$) as $\epsilon\rightarrow 0$. For instance, we may choose $$\delta_\epsilon(f)(x_1,\ldots,x_n)=g_\epsilon(x_1-x_2)\cdots g_\epsilon(x_{n-1}-x_n)f(x_n)\ ,$$ where $g_\epsilon(x)=\frac{1}{\epsilon^4}g\left(\frac{x}{\epsilon}\right)$ and $0\leq g\in\mathscr{D}(\mathbb{R}^4)$ is such that $g(0)=1$ and $\int_{\mathbb{R}^4}g(x)dx=1$ - indeed, we have for all $h\in\mathscr{D}(\mathbb{R}^{4n})$ that $$\delta_\epsilon(f)(h)=\int_{\mathbb{R}^{4n}}g_\epsilon(x_1-x_2)\cdots g_\epsilon(x_{n-1}-x_n)f(x_n)h(x_1,\ldots,x_n)dx_1\cdots dx_n$$ $$\stackrel{\epsilon\rightarrow 0}{\longrightarrow}\int_{\mathbb{R}^4}f(x_1)h(x_1,\ldots,x_1)dx_1\ .$$ What one shows is that the vacuum expectation values of products of auxiliary Wick monomials $\langle\Omega,:\!\phi^{\otimes n}\!:\!\!(x_{1,1},\ldots,x_{1,n})\cdots:\!\phi^{\otimes n}\!:\!\!(x_{k,1},\ldots,x_{k,n})\Omega\rangle$$\langle\Omega,:\!\phi^{\otimes n}\!\!:\!(x_{1,1},\ldots,x_{1,n})\cdots:\!\phi^{\otimes n}\!\!:\!(x_{k,1},\ldots,x_{k,n})\Omega\rangle$ is zero for $k$ odd and, for $k$ even, a linear combination of expressions of the form $\omega_2(x_{i_1,m_{i_1}},x_{j_1,n_{j_1}})\cdots\omega_2(x_{i_{\frac{k}{2}},m_{i_{\frac{k}{2}}}},x_{j_{\frac{k}{2}},n_{j_{\frac{k}{2}}}})$$\omega_2\left(x_{i_1,m_{i_1}},x_{j_1,n_{j_1}}\right)\cdots\omega_2\left(x_{i_{\frac{k}{2}},m_{i_{\frac{k}{2}}}},x_{j_{\frac{k}{2}},n_{j_{\frac{k}{2}}}}\right)$ as one runs over all partitions of $\{1,\ldots,k\}$ into two-element (disjoint) subsets $\{i_1,j_1\},\ldots,\{i_{\frac{k}{2}},j_{\frac{k}{2}}\}$ with $i_l<j_l$ for all $l=1,\ldots,\frac{k}{2}$ (i.e. all pairings of elements of $\{1,\ldots,k\}$ - notice that the set of all pairings of elements of $\{1,\ldots,k\}$ is empty if $k$ is odd). Notice that each such expression formally becomes $\omega_2(x_{i_1},x_{j_1})^n\cdots\omega_2(x_{i_{\frac{k}{2}}},x_{j_{\frac{k}{2}}})^n$$\omega_2\left(x_{i_1},x_{j_1}\right)^n\cdots\omega_2\left(x_{i_{\frac{k}{2}}},x_{j_{\frac{k}{2}}}\right)^n$ if $x_{i,1}=\cdots=x_{i,n}=x_i$ for all $i=1,\ldots,k$. The point is that the latter expression actually is well defined as a $k$-point distribution, thanks to the wave front set characterization of $\omega_2$ provided by 1. above and Theorem 8.2.10, pp. 267 ( = wave front set criterion of existence of the pointwise product of a pair of distributions) of Hörmander's book quoted by the OP, so one may take the limit $\epsilon\rightarrow 0$ to define the $k$-point vacuum expectation values of $:\!\phi^n\!:$ as $\lim_{\epsilon\rightarrow 0}\langle\Omega,:\!\phi^{\otimes n}_\epsilon\!:\!\!(f_1)\cdots:\!\phi^{\otimes n}_\epsilon\!:\!\!(f_k)\Omega\rangle$$\lim_{\epsilon\rightarrow 0}\langle\Omega,:\!\phi^{\otimes n}_\epsilon\!\!:\!(f_1)\cdots:\!\phi^{\otimes n}_\epsilon\!\!:\!(f_k)\Omega\rangle$. This allows us to define $:\!\phi^n\!:\!\!(f)$$:\!\phi^n\!\!:\!(f)$ as an operator in (at first sight) its own Hilbert space $\mathscr{H}_n$ thanks to the Wightman-GNS reconstruction theorem, so that we identify the above limit with $\langle\Omega_n,:\!\phi^n\!:\!\!(f_1),\cdots:\!\phi^n\!:\!\!(f_k)\Omega_n\rangle$$\langle\Omega_n,:\!\phi^n\!\!:\!(f_1),\cdots:\!\phi^n\!\!:\!(f_k)\Omega_n\rangle$ for some nonzero state vector $\Omega_n\in\mathscr{H}_n$. However, one can also show that we must have $\mathscr{H}_n=\mathscr{H}$ and $\Omega_n=\Omega$ so $:\!\phi^n\!:\!\!(f)$$:\!\phi^n\!\!:\!(f)$ exists as an (unbounded) operator on $\mathscr{H}$ with domain $\mathscr{D}$. All this - and more - is contained in Proposition 5.3, pp. 647 of the paper quoted in 1. above. In fact, one may construct in a similar way the $k$-point vacuum expectation values of different Wick powers of $\phi$ (see formula (16), pp. 647 of ibid.), which on their turn may be used to conclude that the matrix elements of $:\!\phi^n\!:\!\!(x)$$:\!\phi^n\!\!:\!(x)$ in $\mathscr{D}$ must be smooth in $x$ (possibly zero) - this is due to the fact that if this matrix element is not zero, then it must be a linear combination of products of copies of $\omega_2$ where all such factors have eachexpressions of their two variablesthe form either equal to $x$$\omega_2^n(x,f)\omega_2(f_1,g_1)\cdots\omega_2(f_m,g_m)$ or smeared with a test function$\omega_2^n(f,x)\omega_2(f_1,g_1)\cdots\omega_2(f_m,g_m)$, with the first case happeningwhere at most once per factor$f,f_1,\ldots,f_m,g_1,\ldots,g_m$ are test functions and precisely $n$ times overall$m\in\mathbb{N}\cup\{0\}$. By the wave front set characterization of $\omega_2$ provided in 1. above and Theorem 8.2.12, pp. 268 of Hörmander's book quoted by the OP, thiseach such product must be smooth in $x$.
Notice that the above procedure is not the same as pulling back the matrix elements $\langle\Omega_1,\!\phi^{\otimes n}\!:\!\!(x_1,\ldots,x_n)\Omega_2\rangle$$\langle\Omega_1,\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)\Omega_2\rangle$ of the auxiliary Wick monomial $:\!\phi^{\otimes n}\!:\!\!(x_1,\ldots,x_n)$$:\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)$ between two vectors $\Omega_1,\Omega_2\in\mathscr{D}$ by $\iota_n$ as suggested by the OP in the case $n=2$. The difference is the following: one takes the pullback by $\iota_n$ in the way described above, by first mollifying the auxiliary Wick monomials with $\delta_\epsilon$ before smearing the arguments of $\Omega_1,\Omega_2$ with test functions. Only then can the necessary singularity cancellations take place as the limit $\epsilon\rightarrow 0$ is performed at the end. More precisely, this procedure guarantees that one never makes both arguments of a single factor $\omega_2$ in the Wick formula coincide in that limit, which (as argued by 1. and 2. above) is something that cannot be done. The way proposed by the OP falls into that very trap - in other words, restricting $G$ to $\Delta_2$ (more precisely, pulling back $G$ by $\iota_2$) is not a well defined procedure by the same reason as restricting $F$ to $\Delta_2$ is not.
(Edit in response to updates in OP) I will provide some of the details below as guidance to the reader. Let $\mathscr{H}$ be the Klein-Gordon Fock (Hilbert) space, which is the completion of the linear span $\mathscr{D}$ of vectors of the form $\phi(f_1)\cdots\phi(f_m)\Omega$, where $f_1,\ldots,f_m$ are test functions on $\mathbb{R}^4$ and $m\in\mathbb{N}$ and $\Omega$ is the vaccum state vector, just like in 1. above. In other words, $\mathscr{D}$ is the common domain of definition of the smeared Klein-Gordon field operators $\phi(f)$. One first defines the so-called auxiliary Wick monomials recursively as the following sequence of distributions $:\!\phi^{\otimes n}\!:$ in $\mathbb{R}^{4n}$ taking values in unbounded linear operators on $\mathscr{H}$ with common domain $\mathscr{D}$:
$$:\!\phi^{\otimes 0}\!:\,\,=\,\mathbb{1}\ ,\,\,:\!\phi^{\otimes 1}\!:\!\!(x)=\phi(x)\ ;$$ $$:\!\phi^{\otimes(n+1)}\!:\!\!(x_1,\ldots,x_{n+1})=\,\,:\!\phi^{\otimes n}\!:\!\!(x_1,\ldots,x_n)\phi(x_{n+1})$$ $$\phantom{\!\phi^{\otimes(n+1)}\!:\!\!(x_1,\ldots,x_{n+1})=\,\,}-\sum^n_{k=1}\omega_2(x_k,x_{n+1})\,\!\!:\!\phi^{\otimes(n-1)}\!:\!\!(x_1,\ldots,\widehat{x_k},\ldots,x_n)\ ,\,n\in\mathbb{N}\ ,$$
where the hat over an argument denotes its omission. This is exactly the standard formula for the Wick monomials $:\!\phi^n\!:\!\!(x)$ prior to the restriction to the small diagonal $\Delta_n=\{(x,\ldots,x)\ |\ x\in\mathbb{R}^4\}\subset\mathbb{R}^{4n}$ (so that $\Delta_2=\Delta$). To perform the latter operation, we define the mollified Wick monomials $:\!\phi^{\otimes n}_\epsilon\!:$ as $$:\!\phi^{\otimes n}_\epsilon\!:(f)=:\!\phi^{\otimes n}\!:(\delta_\epsilon(f))\ ,\,n>1\ ,f\in\mathscr{D}(\mathbb{R}^4)\ ,$$ where $\delta_\epsilon:\mathscr{D}(\mathbb{R}^4)\rightarrow\mathscr{D}'(\mathbb{R}^{4n})$ is a continuous linear map taking values in $\mathscr{D}(\mathbb{R}^{4n})$ for each $\epsilon>0$ and converging in $\mathscr{D}'(\mathbb{R}^{4n})$ to the formal adjoint ( = transpose) of the pullback $\iota^*_n$ by the diagonal map $\iota_n:\mathbb{R}^4\rightarrow\mathbb{R}^{4n}$, $\iota_n(x)=(x,\ldots,x)$ (so that $\iota_n(\mathbb{R}^4)=\Delta_n$) as $\epsilon\rightarrow 0$. For instance, we may choose $$\delta_\epsilon(f)(x_1,\ldots,x_n)=g_\epsilon(x_1-x_2)\cdots g_\epsilon(x_{n-1}-x_n)f(x_n)\ ,$$ where $g_\epsilon(x)=\frac{1}{\epsilon^4}g\left(\frac{x}{\epsilon}\right)$ and $0\leq g\in\mathscr{D}(\mathbb{R}^4)$ is such that $g(0)=1$ and $\int_{\mathbb{R}^4}g(x)dx=1$ - indeed, we have for all $h\in\mathscr{D}(\mathbb{R}^{4n})$ that $$\delta_\epsilon(f)(h)=\int_{\mathbb{R}^{4n}}g_\epsilon(x_1-x_2)\cdots g_\epsilon(x_{n-1}-x_n)f(x_n)h(x_1,\ldots,x_n)dx_1\cdots dx_n$$ $$\stackrel{\epsilon\rightarrow 0}{\longrightarrow}\int_{\mathbb{R}^4}f(x_1)h(x_1,\ldots,x_1)dx_1\ .$$ What one shows is that the vacuum expectation values of products of auxiliary Wick monomials $\langle\Omega,:\!\phi^{\otimes n}\!:\!\!(x_{1,1},\ldots,x_{1,n})\cdots:\!\phi^{\otimes n}\!:\!\!(x_{k,1},\ldots,x_{k,n})\Omega\rangle$ is zero for $k$ odd and, for $k$ even, a linear combination of expressions of the form $\omega_2(x_{i_1,m_{i_1}},x_{j_1,n_{j_1}})\cdots\omega_2(x_{i_{\frac{k}{2}},m_{i_{\frac{k}{2}}}},x_{j_{\frac{k}{2}},n_{j_{\frac{k}{2}}}})$ as one runs over all partitions of $\{1,\ldots,k\}$ into two-element (disjoint) subsets $\{i_1,j_1\},\ldots,\{i_{\frac{k}{2}},j_{\frac{k}{2}}\}$ with $i_l<j_l$ for all $l=1,\ldots,\frac{k}{2}$ (i.e. all pairings of elements of $\{1,\ldots,k\}$ - notice that the set of all pairings of elements of $\{1,\ldots,k\}$ is empty if $k$ is odd). Notice that each such expression formally becomes $\omega_2(x_{i_1},x_{j_1})^n\cdots\omega_2(x_{i_{\frac{k}{2}}},x_{j_{\frac{k}{2}}})^n$ if $x_{i,1}=\cdots=x_{i,n}=x_i$ for all $i=1,\ldots,k$. The point is that the latter expression actually is well defined as a $k$-point distribution, thanks to the wave front set characterization of $\omega_2$ provided by 1. above and Theorem 8.2.10, pp. 267 ( = wave front set criterion of existence of the pointwise product of a pair of distributions) of Hörmander's book quoted by the OP, so one may take the limit $\epsilon\rightarrow 0$ to define the $k$-point vacuum expectation values of $:\!\phi^n\!:$ as $\lim_{\epsilon\rightarrow 0}\langle\Omega,:\!\phi^{\otimes n}_\epsilon\!:\!\!(f_1)\cdots:\!\phi^{\otimes n}_\epsilon\!:\!\!(f_k)\Omega\rangle$. This allows us to define $:\!\phi^n\!:\!\!(f)$ as an operator in (at first sight) its own Hilbert space $\mathscr{H}_n$ thanks to the Wightman-GNS reconstruction theorem, so that we identify the above limit with $\langle\Omega_n,:\!\phi^n\!:\!\!(f_1),\cdots:\!\phi^n\!:\!\!(f_k)\Omega_n\rangle$ for some nonzero state vector $\Omega_n\in\mathscr{H}_n$. However, one can also show that we must have $\mathscr{H}_n=\mathscr{H}$ and $\Omega_n=\Omega$ so $:\!\phi^n\!:\!\!(f)$ exists as an (unbounded) operator on $\mathscr{H}$ with domain $\mathscr{D}$. All this - and more - is contained in Proposition 5.3, pp. 647 of the paper quoted in 1. above. In fact, one may construct in a similar way the $k$-point vacuum expectation values of different Wick powers of $\phi$ (see formula (16), pp. 647 of ibid.), which on their turn may be used to conclude that the matrix elements of $:\!\phi^n\!:\!\!(x)$ in $\mathscr{D}$ must be smooth in $x$ (possibly zero) - this is due to the fact that if this matrix element is not zero, then it must be a linear combination of products of copies of $\omega_2$ where all such factors have each of their two variables either equal to $x$ or smeared with a test function, with the first case happening at most once per factor and precisely $n$ times overall. By the wave front set characterization of $\omega_2$ provided in 1. above and Theorem 8.2.12, pp. 268 of Hörmander's book quoted by the OP, this product must be smooth in $x$.
Notice that the above procedure is not the same as pulling back the matrix elements $\langle\Omega_1,\!\phi^{\otimes n}\!:\!\!(x_1,\ldots,x_n)\Omega_2\rangle$ of the auxiliary Wick monomial $:\!\phi^{\otimes n}\!:\!\!(x_1,\ldots,x_n)$ between two vectors $\Omega_1,\Omega_2\in\mathscr{D}$ by $\iota_n$ as suggested by the OP in the case $n=2$. The difference is the following: one takes the pullback by $\iota_n$ in the way described above, by first mollifying the auxiliary Wick monomials with $\delta_\epsilon$ before smearing the arguments of $\Omega_1,\Omega_2$ with test functions. Only then can the necessary singularity cancellations take place as the limit $\epsilon\rightarrow 0$ is performed at the end. More precisely, this procedure guarantees that one never makes both arguments of a single factor $\omega_2$ in the Wick formula coincide in that limit, which (as argued by 1. and 2. above) is something that cannot be done. The way proposed by the OP falls into that very trap - in other words, restricting $G$ to $\Delta_2$ (more precisely, pulling back $G$ by $\iota_2$) is not a well defined procedure by the same reason as restricting $F$ to $\Delta_2$ is not.
(Edit in response to updates in OP) I will provide some of the details below as guidance to the reader. Let $\mathscr{H}$ be the Klein-Gordon Fock (Hilbert) space, which is the completion of the linear span $\mathscr{D}$ of vectors of the form $\phi(f_1)\cdots\phi(f_m)\Omega$, where $f_1,\ldots,f_m$ are test functions on $\mathbb{R}^4$ and $m\in\mathbb{N}$ and $\Omega$ is the vaccum state vector, just like in 1. above. In other words, $\mathscr{D}$ is the common domain of definition of the smeared Klein-Gordon field operators $\phi(f)$. One first defines the so-called auxiliary Wick monomials recursively as the following sequence of distributions $:\!\phi^{\otimes n}\!\!:$ in $\mathbb{R}^{4n}$ taking values in unbounded linear operators on $\mathscr{H}$ with common domain $\mathscr{D}$:
$$:\!\phi^{\otimes 0}\!\!:\,\,=\,\mathbb{1}\ ,\,\,:\!\phi^{\otimes 1}\!\!:\!(x)=\phi(x)\ ;$$ $$:\!\phi^{\otimes(n+1)}\!\!:\!(x_1,\ldots,x_{n+1})=\,\,:\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)\phi(x_{n+1})$$ $$\phantom{\!\phi^{\otimes(n+1)}\!\!:\!(x_1,\ldots,x_{n+1})=\,\,}-\sum^n_{k=1}\omega_2(x_k,x_{n+1})\,\!\!:\!\phi^{\otimes(n-1)}\!\!:\!(x_1,\ldots,\widehat{x_k},\ldots,x_n)\ ,\,n\in\mathbb{N}\ ,$$
where the hat over an argument denotes its omission. This is exactly the standard formula for the Wick monomials $:\!\phi^n\!:\!\!(x)$ prior to the restriction to the small diagonal $\Delta_n=\{(x,\ldots,x)\ |\ x\in\mathbb{R}^4\}\subset\mathbb{R}^{4n}$ (so that $\Delta_2=\Delta$). To perform the latter operation, we define the mollified Wick monomials $:\!\phi^{\otimes n}_\epsilon\!:$ as $$:\!\phi^{\otimes n}_\epsilon\!\!:\!(f)=:\!\phi^{\otimes n}\!\!:\!(\delta_\epsilon(f))\ ,\,n>1\ ,f\in\mathscr{D}(\mathbb{R}^4)\ ,$$ where $\delta_\epsilon:\mathscr{D}(\mathbb{R}^4)\rightarrow\mathscr{D}'(\mathbb{R}^{4n})$ is a continuous linear map taking values in $\mathscr{D}(\mathbb{R}^{4n})$ for each $\epsilon>0$ and converging in $\mathscr{D}'(\mathbb{R}^{4n})$ to the formal adjoint ( = transpose) of the pullback $\iota^*_n$ by the diagonal map $\iota_n:\mathbb{R}^4\rightarrow\mathbb{R}^{4n}$, $\iota_n(x)=(x,\ldots,x)$ (so that $\iota_n(\mathbb{R}^4)=\Delta_n$) as $\epsilon\rightarrow 0$. For instance, we may choose $$\delta_\epsilon(f)(x_1,\ldots,x_n)=g_\epsilon(x_1-x_2)\cdots g_\epsilon(x_{n-1}-x_n)f(x_n)\ ,$$ where $g_\epsilon(x)=\frac{1}{\epsilon^4}g\left(\frac{x}{\epsilon}\right)$ and $0\leq g\in\mathscr{D}(\mathbb{R}^4)$ is such that $g(0)=1$ and $\int_{\mathbb{R}^4}g(x)dx=1$ - indeed, we have for all $h\in\mathscr{D}(\mathbb{R}^{4n})$ that $$\delta_\epsilon(f)(h)=\int_{\mathbb{R}^{4n}}g_\epsilon(x_1-x_2)\cdots g_\epsilon(x_{n-1}-x_n)f(x_n)h(x_1,\ldots,x_n)dx_1\cdots dx_n$$ $$\stackrel{\epsilon\rightarrow 0}{\longrightarrow}\int_{\mathbb{R}^4}f(x_1)h(x_1,\ldots,x_1)dx_1\ .$$ What one shows is that the vacuum expectation values of products of auxiliary Wick monomials $\langle\Omega,:\!\phi^{\otimes n}\!\!:\!(x_{1,1},\ldots,x_{1,n})\cdots:\!\phi^{\otimes n}\!\!:\!(x_{k,1},\ldots,x_{k,n})\Omega\rangle$ is zero for $k$ odd and, for $k$ even, a linear combination of expressions of the form $\omega_2\left(x_{i_1,m_{i_1}},x_{j_1,n_{j_1}}\right)\cdots\omega_2\left(x_{i_{\frac{k}{2}},m_{i_{\frac{k}{2}}}},x_{j_{\frac{k}{2}},n_{j_{\frac{k}{2}}}}\right)$ as one runs over all partitions of $\{1,\ldots,k\}$ into two-element (disjoint) subsets $\{i_1,j_1\},\ldots,\{i_{\frac{k}{2}},j_{\frac{k}{2}}\}$ with $i_l<j_l$ for all $l=1,\ldots,\frac{k}{2}$ (i.e. all pairings of elements of $\{1,\ldots,k\}$ - notice that the set of all pairings of elements of $\{1,\ldots,k\}$ is empty if $k$ is odd). Notice that each such expression formally becomes $\omega_2\left(x_{i_1},x_{j_1}\right)^n\cdots\omega_2\left(x_{i_{\frac{k}{2}}},x_{j_{\frac{k}{2}}}\right)^n$ if $x_{i,1}=\cdots=x_{i,n}=x_i$ for all $i=1,\ldots,k$. The point is that the latter expression actually is well defined as a $k$-point distribution, thanks to the wave front set characterization of $\omega_2$ provided by 1. above and Theorem 8.2.10, pp. 267 ( = wave front set criterion of existence of the pointwise product of a pair of distributions) of Hörmander's book quoted by the OP, so one may take the limit $\epsilon\rightarrow 0$ to define the $k$-point vacuum expectation values of $:\!\phi^n\!:$ as $\lim_{\epsilon\rightarrow 0}\langle\Omega,:\!\phi^{\otimes n}_\epsilon\!\!:\!(f_1)\cdots:\!\phi^{\otimes n}_\epsilon\!\!:\!(f_k)\Omega\rangle$. This allows us to define $:\!\phi^n\!\!:\!(f)$ as an operator in (at first sight) its own Hilbert space $\mathscr{H}_n$ thanks to the Wightman-GNS reconstruction theorem, so that we identify the above limit with $\langle\Omega_n,:\!\phi^n\!\!:\!(f_1),\cdots:\!\phi^n\!\!:\!(f_k)\Omega_n\rangle$ for some nonzero state vector $\Omega_n\in\mathscr{H}_n$. However, one can also show that we must have $\mathscr{H}_n=\mathscr{H}$ and $\Omega_n=\Omega$ so $:\!\phi^n\!\!:\!(f)$ exists as an (unbounded) operator on $\mathscr{H}$ with domain $\mathscr{D}$. All this - and more - is contained in Proposition 5.3, pp. 647 of the paper quoted in 1. above. In fact, one may construct in a similar way the $k$-point vacuum expectation values of different Wick powers of $\phi$ (see formula (16), pp. 647 of ibid.), which on their turn may be used to conclude that the matrix elements of $:\!\phi^n\!\!:\!(x)$ in $\mathscr{D}$ must be smooth in $x$ (possibly zero) - this is due to the fact that if this matrix element is not zero, then it must be a linear combination of expressions of the form either $\omega_2^n(x,f)\omega_2(f_1,g_1)\cdots\omega_2(f_m,g_m)$ or $\omega_2^n(f,x)\omega_2(f_1,g_1)\cdots\omega_2(f_m,g_m)$, where $f,f_1,\ldots,f_m,g_1,\ldots,g_m$ are test functions and $m\in\mathbb{N}\cup\{0\}$. By the wave front set characterization of $\omega_2$ provided in 1. above and Theorem 8.2.12, pp. 268 of Hörmander's book quoted by the OP, each such product must be smooth in $x$.
Notice that the above procedure is not the same as pulling back the matrix elements $\langle\Omega_1,\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)\Omega_2\rangle$ of the auxiliary Wick monomial $:\!\phi^{\otimes n}\!\!:\!(x_1,\ldots,x_n)$ between two vectors $\Omega_1,\Omega_2\in\mathscr{D}$ by $\iota_n$ as suggested by the OP in the case $n=2$. The difference is the following: one takes the pullback by $\iota_n$ in the way described above, by first mollifying the auxiliary Wick monomials with $\delta_\epsilon$ before smearing the arguments of $\Omega_1,\Omega_2$ with test functions. Only then can the necessary singularity cancellations take place as the limit $\epsilon\rightarrow 0$ is performed at the end. More precisely, this procedure guarantees that one never makes both arguments of a single factor $\omega_2$ in the Wick formula coincide in that limit, which (as argued by 1. above) is something that cannot be done. The way proposed by the OP falls into that very trap - in other words, restricting $G$ to $\Delta_2$ (more precisely, pulling back $G$ by $\iota_2$) is not a well defined procedure by the same reason as restricting $F$ to $\Delta_2$ is not.
where the hat over an argument denotes its omission. This is exactly the standard Wick formula for the Wick monomials $:\!\phi^n\!:\!\!(x)$ prior toprior to the restriction to the small diagonal $\Delta_n=\{(x,\ldots,x)\ |\ x\in\mathbb{R}^4\}\subset\mathbb{R}^{4n}$ (so that $\Delta_2=\Delta$). To perform the latter operation, we define the mollified Wick monomials $:\!\phi^{\otimes n}_\epsilon\!:$ as $$:\!\phi^{\otimes n}_\epsilon\!:(f)=:\!\phi^{\otimes n}\!:(\delta_\epsilon(f))\ ,\,n>1\ ,f\in\mathscr{D}(\mathbb{R}^4)\ ,$$ where $\delta_\epsilon:\mathscr{D}(\mathbb{R}^4)\rightarrow\mathscr{D}'(\mathbb{R}^{4n})$ is a continuous linear map taking values in $\mathscr{D}(\mathbb{R}^{4n})$ for each $\epsilon>0$ and converging in $\mathscr{D}'(\mathbb{R}^{4n})$ to the formal adjoint ( = transpose) of the pullback $\iota^*_n$ by the diagonal map $\iota_n:\mathbb{R}^4\rightarrow\mathbb{R}^{4n}$, $\iota_n(x)=(x,\ldots,x)$ (so that $\iota_n(\mathbb{R}^4)=\Delta_n$) as $\epsilon\rightarrow 0$. For instance, we may choose $$\delta_\epsilon(f)(x_1,\ldots,x_n)=g_\epsilon(x_1-x_2)\cdots g_\epsilon(x_{n-1}-x_n)f(x_n)\ ,$$ where $g_\epsilon(x)=\frac{1}{\epsilon^4}g\left(\frac{x}{\epsilon}\right)$ and $0\leq g\in\mathscr{D}(\mathbb{R}^4)$ is such that $g(0)=1$ and $\int_{\mathbb{R}^4}g(x)dx=1$ - indeed, we have for all $h\in\mathscr{D}(\mathbb{R}^{4n})$ that $$\delta_\epsilon(f)(h)=\int_{\mathbb{R}^{4n}}g_\epsilon(x_1-x_2)\cdots g_\epsilon(x_{n-1}-x_n)f(x_n)h(x_1,\ldots,x_n)dx_1\cdots dx_n$$ $$\stackrel{\epsilon\rightarrow 0}{\longrightarrow}\int_{\mathbb{R}^4}f(x_1)h(x_1,\ldots,x_1)dx_1\ .$$ What one shows is that the vacuum expectation values of products of auxiliary Wick monomials $\langle\Omega,:\!\phi^{\otimes n}\!:\!\!(x_{1,1},\ldots,x_{1,n})\cdots:\!\phi^{\otimes n}\!:\!\!(x_{k,1},\ldots,x_{k,n})\Omega\rangle$ is zero for $k$ odd and, for $k$ even, a linear combination of expressions of the form $\omega_2(x_{i_1,m_{i_1}},x_{j_1,n_{j_1}})\cdots\omega_2(x_{i_{\frac{k}{2}},m_{i_{\frac{k}{2}}}},x_{j_{\frac{k}{2}},n_{j_{\frac{k}{2}}}})$ as one runs over all partitions of $\{1,\ldots,k\}$ into two-element (disjoint) subsets $\{i_1,j_1\},\ldots,\{i_{\frac{k}{2}},j_{\frac{k}{2}}\}$ with $i_l<j_l$ for all $l=1,\ldots,\frac{k}{2}$ (i.e. all pairings of elements of $\{1,\ldots,k\}$ - notice that the set of all pairings of elements of $\{1,\ldots,k\}$ is empty if $k$ is odd). Notice that each such expression formally becomes $\omega_2(x_{i_1},x_{j_1})^n\cdots\omega_2(x_{i_{\frac{k}{2}}},x_{j_{\frac{k}{2}}})^n$ if $x_{i,1}=\cdots=x_{i,n}=x_i$ for all $i=1,\ldots,k$. The point is that the latter expression actually is well defined as a $k$-point distribution, thanks to the wave front set characterization of $\omega_2$ provided by 1. above and Theorem 8.2.10, pp. 267 ( = wave front set criterion of existence of the pointwise product of a pair of distributions) of Hörmander's book quoted by the OP, so one may take the limit $\epsilon\rightarrow 0$ to define the $k$-point vacuum expectation values of $:\!\phi^n\!:$ as $\lim_{\epsilon\rightarrow 0}\langle\Omega,:\!\phi^{\otimes n}_\epsilon\!:\!\!(f_1)\cdots:\!\phi^{\otimes n}_\epsilon\!:\!\!(f_k)\Omega\rangle$. This allows us to define $:\!\phi^n\!:\!\!(f)$ as an operator in (at first sight) its own Hilbert space $\mathscr{H}_n$ thanks to the Wightman-GNS reconstruction theorem, so that we identify the above limit with $\langle\Omega_n,:\!\phi^n\!:\!\!(f_1),\cdots:\!\phi^n\!:\!\!(f_k)\Omega_n\rangle$ for some nonzero state vector $\Omega_n\in\mathscr{H}_n$. However, one can also show that we must have $\mathscr{H}_n=\mathscr{H}$ and $\Omega_n=\Omega$ so $:\!\phi^n\!:\!\!(f)$ exists as an (unbounded) operator on $\mathscr{H}$ with domain $\mathscr{D}$. All this - and more - is contained in Proposition 5.3, pp. 647 of the paper quoted in 1. above. In fact, one may construct in a similar way the $k$-point vacuum expectation values of different Wick powers of $\phi$ (see formula (16), pp. 647 of ibid.), which on their turn may be used to conclude that the matrix elements of $:\!\phi^n\!:\!\!(x)$ in $\mathscr{D}$ must be smooth in $x$ (possibly zero) - this is due to the fact that if this matrix element is not zero, then it must be a linear combination of products of copies of $\omega_2$ where all such factors have each of their two variables either equal to $x$ or smeared with a test function, with the first case happening at most once per factor and precisely $n$ times overall. By the wave front set characterization of $\omega_2$ provided in 1. above and Theorem 8.2.12, pp. 268 of Hörmander's book quoted by the OP, this product must be smooth in $x$.
Notice that the above procedure is not the same as pulling back the matrix elements $\langle\Omega_1,\!\phi^{\otimes n}\!:\!\!(x_1,\ldots,x_n)\Omega_2\rangle$ of the auxiliary Wick monomial $:\!\phi^{\otimes n}\!:\!\!(x_1,\ldots,x_n)$ between two vectors $\Omega_1,\Omega_2\in\mathscr{D}$ by $\iota_n$ as suggested by the OP in the case $n=2$. The difference is the following: one takes the pullback by $\iota_n$ in the way described above, by first mollifying the auxiliary Wick monomials with $\delta_\epsilon$ before smearing the arguments of $\Omega_1,\Omega_2$ with test functions. Only then can the necessary singularity cancellations take place as the limit $\epsilon\rightarrow 0$ is performed at the end. More precisely, this procedure guarantees that one never makes both arguments of a single factor $\omega_2$ in the Wick formula coincide in that limit, which (as argued by 1. and 2. above) is something that cannot be done. The way proposed by the OP falls into that very trap - in other words, restricting $G$ to $\Delta_2$ (more precisely, pulling back $G$ by $\iota_2$) is not a well defined procedure by the same reason as restricting $F$ to $\Delta_2$ is not.
where the hat over an argument denotes its omission. This is exactly the standard Wick formula for the Wick monomials $:\!\phi^n\!:\!\!(x)$ prior to the restriction to the small diagonal $\Delta_n=\{(x,\ldots,x)\ |\ x\in\mathbb{R}^4\}\subset\mathbb{R}^{4n}$ (so that $\Delta_2=\Delta$). To perform the latter operation, we define the mollified Wick monomials $:\!\phi^{\otimes n}_\epsilon\!:$ as $$:\!\phi^{\otimes n}_\epsilon\!:(f)=:\!\phi^{\otimes n}\!:(\delta_\epsilon(f))\ ,\,n>1\ ,f\in\mathscr{D}(\mathbb{R}^4)\ ,$$ where $\delta_\epsilon:\mathscr{D}(\mathbb{R}^4)\rightarrow\mathscr{D}'(\mathbb{R}^{4n})$ is a continuous linear map taking values in $\mathscr{D}(\mathbb{R}^{4n})$ for each $\epsilon>0$ and converging in $\mathscr{D}'(\mathbb{R}^{4n})$ to the formal adjoint ( = transpose) of the pullback $\iota^*_n$ by the diagonal map $\iota_n:\mathbb{R}^4\rightarrow\mathbb{R}^{4n}$, $\iota_n(x)=(x,\ldots,x)$ (so that $\iota_n(\mathbb{R}^4)=\Delta_n$) as $\epsilon\rightarrow 0$. For instance, we may choose $$\delta_\epsilon(f)(x_1,\ldots,x_n)=g_\epsilon(x_1-x_2)\cdots g_\epsilon(x_{n-1}-x_n)f(x_n)\ ,$$ where $g_\epsilon(x)=\frac{1}{\epsilon^4}g\left(\frac{x}{\epsilon}\right)$ and $0\leq g\in\mathscr{D}(\mathbb{R}^4)$ is such that $g(0)=1$ and $\int_{\mathbb{R}^4}g(x)dx=1$ - indeed, we have for all $h\in\mathscr{D}(\mathbb{R}^{4n})$ that $$\delta_\epsilon(f)(h)=\int_{\mathbb{R}^{4n}}g_\epsilon(x_1-x_2)\cdots g_\epsilon(x_{n-1}-x_n)f(x_n)h(x_1,\ldots,x_n)dx_1\cdots dx_n$$ $$\stackrel{\epsilon\rightarrow 0}{\longrightarrow}\int_{\mathbb{R}^4}f(x_1)h(x_1,\ldots,x_1)dx_1\ .$$ What one shows is that the vacuum expectation values of products of auxiliary Wick monomials $\langle\Omega,:\!\phi^{\otimes n}\!:\!\!(x_{1,1},\ldots,x_{1,n})\cdots:\!\phi^{\otimes n}\!:\!\!(x_{k,1},\ldots,x_{k,n})\Omega\rangle$ is zero for $k$ odd and, for $k$ even, a linear combination of expressions of the form $\omega_2(x_{i_1,m_{i_1}},x_{j_1,n_{j_1}})\cdots\omega_2(x_{i_{\frac{k}{2}},m_{i_{\frac{k}{2}}}},x_{j_{\frac{k}{2}},n_{j_{\frac{k}{2}}}})$ as one runs over all partitions of $\{1,\ldots,k\}$ into two-element (disjoint) subsets $\{i_1,j_1\},\ldots,\{i_{\frac{k}{2}},j_{\frac{k}{2}}\}$ with $i_l<j_l$ for all $l=1,\ldots,\frac{k}{2}$ (i.e. all pairings of elements of $\{1,\ldots,k\}$ - notice that the set of all pairings of elements of $\{1,\ldots,k\}$ is empty if $k$ is odd). Notice that each such expression formally becomes $\omega_2(x_{i_1},x_{j_1})^n\cdots\omega_2(x_{i_{\frac{k}{2}}},x_{j_{\frac{k}{2}}})^n$ if $x_{i,1}=\cdots=x_{i,n}=x_i$ for all $i=1,\ldots,k$. The point is that the latter expression actually is well defined as a $k$-point distribution, thanks to the wave front set characterization of $\omega_2$ provided by 1. above and Theorem 8.2.10, pp. 267 ( = wave front set criterion of existence of the pointwise product of a pair of distributions) of Hörmander's book quoted by the OP, so one may take the limit $\epsilon\rightarrow 0$ to define the $k$-point vacuum expectation values of $:\!\phi^n\!:$ as $\lim_{\epsilon\rightarrow 0}\langle\Omega,:\!\phi^{\otimes n}_\epsilon\!:\!\!(f_1)\cdots:\!\phi^{\otimes n}_\epsilon\!:\!\!(f_k)\Omega\rangle$. This allows us to define $:\!\phi^n\!:\!\!(f)$ as an operator in (at first sight) its own Hilbert space $\mathscr{H}_n$ thanks to the Wightman-GNS reconstruction theorem, so that we identify the above limit with $\langle\Omega_n,:\!\phi^n\!:\!\!(f_1),\cdots:\!\phi^n\!:\!\!(f_k)\Omega_n\rangle$ for some nonzero state vector $\Omega_n\in\mathscr{H}_n$. However, one can also show that we must have $\mathscr{H}_n=\mathscr{H}$ and $\Omega_n=\Omega$ so $:\!\phi^n\!:\!\!(f)$ exists as an (unbounded) operator on $\mathscr{H}$ with domain $\mathscr{D}$. All this - and more - is contained in Proposition 5.3, pp. 647 of the paper quoted in 1. above. In fact, one may construct in a similar way the $k$-point vacuum expectation values of different Wick powers of $\phi$ (see formula (16), pp. 647 of ibid.), which on their turn may be used to conclude that the matrix elements of $:\!\phi^n\!:\!\!(x)$ in $\mathscr{D}$ must be smooth in $x$ (possibly zero) - this is due to the fact that if this matrix element is not zero, then it must be a linear combination of products of copies of $\omega_2$ where all such factors have each of their two variables either equal to $x$ or smeared with a test function, with the first case happening at most once per factor and precisely $n$ times overall. By the wave front set characterization of $\omega_2$ provided in 1. above and Theorem 8.2.12, pp. 268 of Hörmander's book quoted by the OP, this product must be smooth in $x$.
where the hat over an argument denotes its omission. This is exactly the standard formula for the Wick monomials $:\!\phi^n\!:\!\!(x)$ prior to the restriction to the small diagonal $\Delta_n=\{(x,\ldots,x)\ |\ x\in\mathbb{R}^4\}\subset\mathbb{R}^{4n}$ (so that $\Delta_2=\Delta$). To perform the latter operation, we define the mollified Wick monomials $:\!\phi^{\otimes n}_\epsilon\!:$ as $$:\!\phi^{\otimes n}_\epsilon\!:(f)=:\!\phi^{\otimes n}\!:(\delta_\epsilon(f))\ ,\,n>1\ ,f\in\mathscr{D}(\mathbb{R}^4)\ ,$$ where $\delta_\epsilon:\mathscr{D}(\mathbb{R}^4)\rightarrow\mathscr{D}'(\mathbb{R}^{4n})$ is a continuous linear map taking values in $\mathscr{D}(\mathbb{R}^{4n})$ for each $\epsilon>0$ and converging in $\mathscr{D}'(\mathbb{R}^{4n})$ to the formal adjoint ( = transpose) of the pullback $\iota^*_n$ by the diagonal map $\iota_n:\mathbb{R}^4\rightarrow\mathbb{R}^{4n}$, $\iota_n(x)=(x,\ldots,x)$ (so that $\iota_n(\mathbb{R}^4)=\Delta_n$) as $\epsilon\rightarrow 0$. For instance, we may choose $$\delta_\epsilon(f)(x_1,\ldots,x_n)=g_\epsilon(x_1-x_2)\cdots g_\epsilon(x_{n-1}-x_n)f(x_n)\ ,$$ where $g_\epsilon(x)=\frac{1}{\epsilon^4}g\left(\frac{x}{\epsilon}\right)$ and $0\leq g\in\mathscr{D}(\mathbb{R}^4)$ is such that $g(0)=1$ and $\int_{\mathbb{R}^4}g(x)dx=1$ - indeed, we have for all $h\in\mathscr{D}(\mathbb{R}^{4n})$ that $$\delta_\epsilon(f)(h)=\int_{\mathbb{R}^{4n}}g_\epsilon(x_1-x_2)\cdots g_\epsilon(x_{n-1}-x_n)f(x_n)h(x_1,\ldots,x_n)dx_1\cdots dx_n$$ $$\stackrel{\epsilon\rightarrow 0}{\longrightarrow}\int_{\mathbb{R}^4}f(x_1)h(x_1,\ldots,x_1)dx_1\ .$$ What one shows is that the vacuum expectation values of products of auxiliary Wick monomials $\langle\Omega,:\!\phi^{\otimes n}\!:\!\!(x_{1,1},\ldots,x_{1,n})\cdots:\!\phi^{\otimes n}\!:\!\!(x_{k,1},\ldots,x_{k,n})\Omega\rangle$ is zero for $k$ odd and, for $k$ even, a linear combination of expressions of the form $\omega_2(x_{i_1,m_{i_1}},x_{j_1,n_{j_1}})\cdots\omega_2(x_{i_{\frac{k}{2}},m_{i_{\frac{k}{2}}}},x_{j_{\frac{k}{2}},n_{j_{\frac{k}{2}}}})$ as one runs over all partitions of $\{1,\ldots,k\}$ into two-element (disjoint) subsets $\{i_1,j_1\},\ldots,\{i_{\frac{k}{2}},j_{\frac{k}{2}}\}$ with $i_l<j_l$ for all $l=1,\ldots,\frac{k}{2}$ (i.e. all pairings of elements of $\{1,\ldots,k\}$ - notice that the set of all pairings of elements of $\{1,\ldots,k\}$ is empty if $k$ is odd). Notice that each such expression formally becomes $\omega_2(x_{i_1},x_{j_1})^n\cdots\omega_2(x_{i_{\frac{k}{2}}},x_{j_{\frac{k}{2}}})^n$ if $x_{i,1}=\cdots=x_{i,n}=x_i$ for all $i=1,\ldots,k$. The point is that the latter expression actually is well defined as a $k$-point distribution, thanks to the wave front set characterization of $\omega_2$ provided by 1. above and Theorem 8.2.10, pp. 267 ( = wave front set criterion of existence of the pointwise product of a pair of distributions) of Hörmander's book quoted by the OP, so one may take the limit $\epsilon\rightarrow 0$ to define the $k$-point vacuum expectation values of $:\!\phi^n\!:$ as $\lim_{\epsilon\rightarrow 0}\langle\Omega,:\!\phi^{\otimes n}_\epsilon\!:\!\!(f_1)\cdots:\!\phi^{\otimes n}_\epsilon\!:\!\!(f_k)\Omega\rangle$. This allows us to define $:\!\phi^n\!:\!\!(f)$ as an operator in (at first sight) its own Hilbert space $\mathscr{H}_n$ thanks to the Wightman-GNS reconstruction theorem, so that we identify the above limit with $\langle\Omega_n,:\!\phi^n\!:\!\!(f_1),\cdots:\!\phi^n\!:\!\!(f_k)\Omega_n\rangle$ for some nonzero state vector $\Omega_n\in\mathscr{H}_n$. However, one can also show that we must have $\mathscr{H}_n=\mathscr{H}$ and $\Omega_n=\Omega$ so $:\!\phi^n\!:\!\!(f)$ exists as an (unbounded) operator on $\mathscr{H}$ with domain $\mathscr{D}$. All this - and more - is contained in Proposition 5.3, pp. 647 of the paper quoted in 1. above. In fact, one may construct in a similar way the $k$-point vacuum expectation values of different Wick powers of $\phi$ (see formula (16), pp. 647 of ibid.), which on their turn may be used to conclude that the matrix elements of $:\!\phi^n\!:\!\!(x)$ in $\mathscr{D}$ must be smooth in $x$ (possibly zero) - this is due to the fact that if this matrix element is not zero, then it must be a linear combination of products of copies of $\omega_2$ where all such factors have each of their two variables either equal to $x$ or smeared with a test function, with the first case happening at most once per factor and precisely $n$ times overall. By the wave front set characterization of $\omega_2$ provided in 1. above and Theorem 8.2.12, pp. 268 of Hörmander's book quoted by the OP, this product must be smooth in $x$.
Notice that the above procedure is not the same as pulling back the matrix elements $\langle\Omega_1,\!\phi^{\otimes n}\!:\!\!(x_1,\ldots,x_n)\Omega_2\rangle$ of the auxiliary Wick monomial $:\!\phi^{\otimes n}\!:\!\!(x_1,\ldots,x_n)$ between two vectors $\Omega_1,\Omega_2\in\mathscr{D}$ by $\iota_n$ as suggested by the OP in the case $n=2$. The difference is the following: one takes the pullback by $\iota_n$ in the way described above, by first mollifying the auxiliary Wick monomials with $\delta_\epsilon$ before smearing the arguments of $\Omega_1,\Omega_2$ with test functions. Only then can the necessary singularity cancellations take place as the limit $\epsilon\rightarrow 0$ is performed at the end. More precisely, this procedure guarantees that one never makes both arguments of a single factor $\omega_2$ in the Wick formula coincide in that limit, which (as argued by 1. and 2. above) is something that cannot be done. The way proposed by the OP falls into that very trap - in other words, restricting $G$ to $\Delta_2$ (more precisely, pulling back $G$ by $\iota_2$) is not a well defined procedure by the same reason as restricting $F$ to $\Delta_2$ is not.