[guided]We now justify that the formal calculation is not merely algebraic. The remainder amplitude has the form
\begin{align*}
R_N(x,y,\xi;h)=\sum_{|\alpha|=N}\frac{N}{\alpha!}\left(\frac{y-x}{2}\right)^\alpha\int_0^1(1-t)^{N-1}\partial_x^\alpha a\left(x+t\frac{y-x}{2},\xi;h\right)\,d\mathcal{L}^1(t).
\end{align*}
The factor $(y-x)^\alpha$ is the important part. Since $|\alpha|=N$, it gives exactly $N$ powers of $y_j-x_j$. For each such factor we use
\begin{align*}
(y_j-x_j)e^{i(x-y)\cdot \xi/h}=ih\,\partial_{\xi_j}e^{i(x-y)\cdot \xi/h}.
\end{align*}
After doing this $N$ times, we integrate by parts in $\xi$ with respect to $d\mathcal{L}^n(\xi)$. The boundary terms vanish in the oscillatory-integral sense after inserting the standard compact $\xi$-cutoff and passing to the limit; the symbol estimates are uniform in the cutoff, so the limit exists in $\mathcal{S}'(\mathbb{R}^n)$.
There is one point that must not be hidden: after the integrations by parts in a Taylor remainder, the amplitude can still depend on $y$. Thus a remainder calculation by itself does not automatically construct a left Kohn--Nirenberg symbol. We therefore apply the standard external amplitude-to-left-symbol reduction theorem directly to the full amplitude $A(x,y,\xi;h)=a((x+y)/2,\xi;h)$. In this theorem, an amplitude means a smooth function of $(x,y,\xi)$, while a left symbol is a smooth function only of $(x,\xi)$. The theorem says that if
\begin{align*}
|\partial_x^\beta\partial_y^\delta\partial_\xi^\gamma A(x,y,\xi;h)|\le C_{\beta\delta\gamma}\langle\xi\rangle^{m-|\gamma|}
\end{align*}
for all $x,y,\xi\in\mathbb{R}^n$ and all multi-indices $\beta,\delta,\gamma$, then the oscillatory integral with amplitude $A$ is $\operatorname{Op}_h(\widetilde b)+\operatorname{Op}_h(c)$ for one symbol $\widetilde b\in S^m_{1,0}(T^*\mathbb{R}^n)$ and one residual left symbol $c$. It also gives the asymptotic formula
\begin{align*}
\widetilde b(x,\xi;h)\sim \sum_\alpha \frac{1}{\alpha!}(-ih)^{|\alpha|}\partial_\xi^\alpha\partial_y^\alpha A(x,y,\xi;h)\big|_{y=x}.
\end{align*}
This is the step that produces a single residual symbol, rather than a family of remainders depending on $N$.
We verify those hypotheses. The differentiated amplitude is smooth in $(x,y,\xi)$ because $a$ is a smooth semiclassical symbol. If $\beta,\delta,\gamma\in\mathbb{N}_0^n$ are multi-indices, then differentiating $a((x+y)/2,\xi;h)$ in $x$ and $y$ only differentiates $a$ in its first variable and multiplies by powers of $1/2$. The $\xi$-derivatives lower the symbolic order exactly as in the definition of $S^{m}_{1,0}$. Therefore
\begin{align*}
|\partial_x^\beta\partial_y^\delta\partial_\xi^\gamma A(x,y,\xi;h)|\le C_{\beta\delta\gamma}\langle\xi\rangle^{m-|\gamma|}
\end{align*}
for constants depending on the corresponding $S^m_{1,0}$ seminorms of $a$. The reduction theorem therefore produces a single symbol $\widetilde b\in S^m_{1,0}(T^*\mathbb{R}^n)$ and a single residual-symbol remainder $c$.
Finally we compute the expansion of $\widetilde b$. Since $\partial_y^\alpha a((x+y)/2,\xi;h)=2^{-|\alpha|}\partial_x^\alpha a((x+y)/2,\xi;h)$, evaluating at $y=x$ gives
\begin{align*}
\partial_y^\alpha A(x,y,\xi;h)\big|_{y=x}=2^{-|\alpha|}\partial_x^\alpha a(x,\xi;h).
\end{align*}
Thus the theorem gives, for every $N\in\mathbb{N}$,
\begin{align*}
\widetilde b-\sum_{|\alpha|<N}\frac{1}{\alpha!}\left(-\frac{ih}{2}\right)^{|\alpha|}\partial_\xi^\alpha\partial_x^\alpha a\in h^N S^{m-N}_{1,0}(T^*\mathbb{R}^n).
\end{align*}[/guided]