I am considering the behavior of $$\frac{1}{h}\|(1-P_h)S(h)v\|,\tag{1}$$ and $$\frac{1}{h}\|(1-P_h)S(h)P_hv\|,\tag{2}$$ as $h\to 0^+$ for a fixed good enough $v$. I hope to show one of them converges to $0$.
Here $\|\cdot\|$ denotes the norm of $L^2(\mathbb{R})$. And given $v\in L^2(\mathbb{R})$, $S(h)v$ denotes the convolution w.r.t Gaussinan kernel defined by $$S(h)v(x)=\int_\mathbb{R} p(h,x,y)v(y)\, dy=\int_\mathbb{R} \frac{1}{\sqrt{2\pi h}}e^{-\frac{|x-y|^2}{2h}}v(y)\, dy,\quad x\in \mathbb{R},$$ and fixing a uniform partition $\{x_j\mid j\in \mathbb{Z}\}$ of $\mathbb{R}$ with $x_{j+1}-x_j=h$, let $P_h$ denotes the $L^2$ projection of $f\in L^2(\mathbb{R})$ on piecewise constant function $$P_hf(x)=\sum_j a_j 1_{(x_j,x_{j+1}]}(x),\; x\in \mathbb{R},\quad a_j=\frac{1}{h}\int_{x_j}^{x_{j+1}}f(y)\,dy.$$
In my opinion, by the observation $$\|(1-P_h)S(h)P_hv\|=\|(1-P_h)(S(h)-I)P_hv\|,$$ (2) may have higher convergence order since intuitively $1-P_h$ and $(S(h)-I)P_hv$ both contribute to the convergence.
Using the following method, I have shown $$\|(1-P_h)S(h)v\|\le Ch^{\frac{1}{2}}\|v\|.$$
Using $$\|(1-P_h)f\|\le h\|f\|_{H^1(\mathbb{R})}$$ it suffices to estimate $\|S(h)v\|_{H^1(\mathbb{R})}$. It's easy to show $\|S(h)v\|\le \|v\|$ for $v\in L^2$, hence we estimate $\|DS(h)v\|$. For each $x\in \mathbb{R}$, we have $$|DS(h)v(x)|\le \int_\mathbb{R} \frac{1}{h\sqrt{2\pi h}}|x-y|e^{-\frac{|x-y|^2}{2h}}|v(y)|\, dy,$$ from which, using Cauchy Schwartz inequality, we derive \begin{align*} |DS(h)v(x)|^2&\le \bigg(\int_\mathbb{R} \frac{1}{h\sqrt{2\pi h}}|x-y|e^{-\frac{|x-y|^2}{2h}}|v(y)|\, dy\bigg)^2\\ &\le \bigg(\int_\mathbb{R} \frac{1}{h\sqrt{2\pi h}}|x-y|e^{-\frac{|x-y|^2}{2h}}|v(y)|^2\, dy\bigg)\bigg(\int_\mathbb{R} \frac{1}{h\sqrt{2\pi h}}|x-y|e^{-\frac{|x-y|^2}{2h}}\, dy\bigg). \end{align*} The above inequality, along with the fact that $\int_\mathbb{R} |x-y|e^{-\frac{|x-y|^2}{2h}}|\, dy=2h$ gives us $$|DS(h)v(x)|^2\lesssim \frac{1}{h^2}\int_\mathbb{R} |x-y|e^{-\frac{|x-y|^2}{2h}}|v(y)|^2\, dy.$$ Integrating w.r.t $x$ and using Fubini's theorem, we obtain $$\|DS(h)v\|^2\le \frac{1}{h}\|v\|^2,$$ which yields $\|S(h)v\|_{H^1(\mathbb{R})}\le \frac{1}{\sqrt{h}}\|v\|$. Using $(1)$, we get the $\frac{1}{2}$ order of convergence.
But this estimate is not good enough to establish the desired results. Moreover this technique can't apply to (2) since $P_hv$ is not $H^1$ function.
John, I apologize for my careless answer, I don't have much time now.
If you compute $P_hS(h)v$, you should get something like
$$ P_hS(h)v(t) =\int_{\mathbb{R}}v(y)\sum_{j} \frac{1}{\sqrt{2\pi h}}\frac{1}{h}\int_{I_j}e^{-\frac{|x-y|^2}{2h}}\,dx\chi_{I_j}(t)\,dy. $$
You may notice that $\sum_{j} \frac{1}{\sqrt{2\pi h}}\frac{1}{h}\int_{I_j}e^{-\frac{|x-y|^2}{2h}}\,dx\chi_{I_j}(t) \approx \frac{1}{\sqrt{2\pi h}}e^{-\frac{|t-y|^2}{2h}}$, then you can think of the operator $P_hS(h)$ as being essentially $S(h)$. Hence you may write $(I-P_h)S(h)v = (S(h)v-v)+(v-P_hS(h)v)$ and it remains to estimate $||(I-S(h))v||$, which is easier.
Notice that $$ |(P_hS(h)-I)v(t)|^2\le \int |v(y)-v(t)|^2\sum_{j} \frac{1}{\sqrt{2\pi h}}\frac{1}{h}\int_{I_j}e^{-\frac{|x-y|^2}{2h}}\,dx\chi_{I_j}(t)\,dy \lesssim \int |v(y)-v(t)|^2\frac{1}{\sqrt{2\pi h}}e^{-\frac{|t-y|^2}{2h}}\,dy. $$
This is a long comment, it may help and tell me if there is a mistake :)