The answer to both questions is no. This is due to two facts:
The Klein-Gordon two-point distribution $\omega_2(x,y)=\langle\Omega,\phi(x)\phi(y)\Omega\rangle$ in $\mathbb{R}^4$, where $\Omega_1=\Omega_2=\Omega$ is the vacuum state vector, has the property that $$WF(\omega_2)\subset\{(x,y;\xi,-\xi)\subset\mathbb{R}^8\times(\mathbb{R}^8\smallsetminus\{(0,0)\})\ |\ \eta(x-y,x-y)=\eta(\xi,\xi)=0\ ,\,\xi_0>0\}\ ,$$ where $\eta$ is the four-dimensional Minkowski metric. This is the so-called microlocal spectrum condition (in Minkowski space-time) and can be traced to two properties of $\omega_2$: (i) $\omega_2$ is a (distributional) bisolution of the Klein-Gordon equation which is invariant by translations in $\mathbb{R}^4$, that is, $\omega_2(x+z,y+z)=\omega(x,y)$ for all $z\in\mathbb{R}^4$, hence by Hörmander's propagation of singularities theorem its wave front set is contained in the superset of the set in the rhs of the above inclusion given by forfeiting the last condition $\xi_0>0$, which implies that $\xi$ should belong the closed forward (i.e. future) light cone; (ii) the spectral condition satisfied by the Klein-Gordon field, which on its turn enforces this latter condition. This is basically Theorem 4.6 , pp. 644-645 of the paper by R. Brunetti, K. Fredenhagen and M. Köhler, The Microlocal Spectrum Condition and Wick Polynomials of Free Fields on Curved Spacetimes, Commun. Math. Phys. 180 (1996) 633-652. This entails that $$WF(\omega_2)\cap N(\Delta)\subset\{(x,x;\xi,-\xi)\ |\ x,\xi\in\mathbb{R}^4\ ,\,\eta(\xi,\xi)=0\}\ .$$ One can show that the same result holds for $\Omega_1,\Omega_2$ being finite linear combinations 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}$ (see e.g. Theorem 4.5, pp. 643-644 of the same paper quoted above), which are dense in the Fock space and form the common domain of the smeared field operators $\phi(f)$ (which, let us recall, are always unbounded).
Let $\delta_\Delta(x,y)=\delta(x-y)$. Then $WF(\delta_\Delta)=N(\Delta)$ by Theorem 8.1,5, pp. 256 of Hörmander's book quoted by the OP.
Both facts together show that both distributions in the rhs of the formula defining $F(x,y)$ cannot mutually cancel all their singular directions conormal to $\Delta$ so as to have $WF(F)\cap N(\Delta)=\varnothing$, so $F(x,y)$ cannot be restricted to $\Delta$.
Anyhow (because this is what I reckon you actually want to know), this is not how you define the matrix elements of the Wick square of $\phi$. A proper microlocal construction of Wick polynomials and their matrix elements can be found in Section 5 of the paper quoted in 1. above.
(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.
Physically, the above procedure is often called point-splitting regularization of Wick products, which displays in a rigorous manner (thanks to the toolset of microlocal analysis) how short-distance singularities of $\phi(x)$ cancel each other out by the Wick ordering prescription so that the pullback by $\iota_n$ may be taken.
