[proofplan]
We rewrite the Weyl kernel in left Kohn--Nirenberg form by expanding the midpoint amplitude $a((x+y)/2,\xi;h)$ around $x$. Powers of $y-x$ are converted into $\xi$-derivatives of the exponential, and integration by parts transfers those derivatives to the symbol with the correct minus sign. The finite Taylor remainders satisfy symbol estimates of lower order, and Borel summation for semiclassical symbols realizes the resulting formal expansion by an actual symbol; the remaining operator is the left quantization of a residual symbol in the sense stated in the theorem.
[/proofplan]
[step:Expand the midpoint amplitude around the left endpoint]
Fix $u\in \mathcal{S}(\mathbb{R}^n)$. For each $N\in\mathbb{N}$, Taylor's formula in the $x$-variable gives, for $x,y,\xi\in\mathbb{R}^n$,
\begin{align*}
a\left(\frac{x+y}{2},\xi;h\right)=\sum_{|\alpha|<N}\frac{1}{\alpha!}\left(\frac{y-x}{2}\right)^\alpha \partial_x^\alpha a(x,\xi;h)+R_N(x,y,\xi;h).
\end{align*}
Here $\alpha=(\alpha_1,\dots,\alpha_n)\in\mathbb{N}_0^n$, $\alpha!=\alpha_1!\cdots\alpha_n!$, $(y-x)^\alpha=\prod_{j=1}^n(y_j-x_j)^{\alpha_j}$, and the remainder is
\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*}
Substituting this expansion into the Weyl formula separates the operator into a finite sum of oscillatory integrals plus a remainder operator with amplitude $R_N$.
[guided]
The only difference between the two quantizations is the point at which the symbol is evaluated. Weyl quantization uses the midpoint $(x+y)/2$, while Kohn--Nirenberg quantization uses the left endpoint $x$. Therefore the natural first move is to Taylor expand the map
\begin{align*}
z\mapsto a(z,\xi;h)
\end{align*}
at $z=x$ and then set $z=(x+y)/2$.
For a multi-index $\alpha=(\alpha_1,\dots,\alpha_n)\in\mathbb{N}_0^n$, write $\partial_x^\alpha=\partial_{x_1}^{\alpha_1}\cdots\partial_{x_n}^{\alpha_n}$ and $(y-x)^\alpha=\prod_{j=1}^n(y_j-x_j)^{\alpha_j}$. Taylor's formula with integral remainder gives
\begin{align*}
a\left(\frac{x+y}{2},\xi;h\right)=\sum_{|\alpha|<N}\frac{1}{\alpha!}\left(\frac{y-x}{2}\right)^\alpha \partial_x^\alpha a(x,\xi;h)+R_N(x,y,\xi;h),
\end{align*}
where
\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 point of writing the remainder in this explicit form is that the factor $(y-x)^\alpha$ has degree $N$. Each factor of $y_j-x_j$ will later become one factor of $h\partial_{\xi_j}$, so the $N$th Taylor remainder contributes only at order $h^N$ after integration by parts.
[/guided]
[/step]
[step:Move powers of $y-x$ onto $\xi$-derivatives of the symbol]
For each $j\in\{1,\dots,n\}$, the phase satisfies
\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*}
Thus, for every multi-index $\alpha$,
\begin{align*}
\left(\frac{y-x}{2}\right)^\alpha e^{i(x-y)\cdot \xi/h}=\left(\frac{ih}{2}\right)^{|\alpha|}\partial_\xi^\alpha e^{i(x-y)\cdot \xi/h}.
\end{align*}
Integrating by parts in $\xi$ with respect to Lebesgue measure $d\mathcal{L}^n(\xi)$ transfers $\partial_\xi^\alpha$ from the exponential to $\partial_x^\alpha a(x,\xi;h)$ and contributes the factor $(-1)^{|\alpha|}$. Hence the coefficient becomes $(-ih/2)^{|\alpha|}$, and the finite Taylor terms become
\begin{align*}
(2\pi h)^{-n}\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}e^{i(x-y)\cdot \xi/h}\sum_{|\alpha|<N}\frac{1}{\alpha!}\left(-\frac{ih}{2}\right)^{|\alpha|}\partial_\xi^\alpha\partial_x^\alpha a(x,\xi;h)u(y)\,d\mathcal{L}^n(y)\,d\mathcal{L}^n(\xi).
\end{align*}
Therefore the finite-order Kohn--Nirenberg symbol is
\begin{align*}
b_N(x,\xi;h)=\sum_{|\alpha|<N}\frac{1}{\alpha!}\left(-\frac{ih}{2}\right)^{|\alpha|}\partial_\xi^\alpha\partial_x^\alpha a(x,\xi;h).
\end{align*}
[/step]
[step:Identify the formal exponential series]
Because
\begin{align*}
\left(\sum_{j=1}^n\partial_{x_j}\partial_{\xi_j}\right)^k=\sum_{|\alpha|=k}\frac{k!}{\alpha!}\partial_x^\alpha\partial_\xi^\alpha,
\end{align*}
the finite symbols $b_N$ are exactly the partial sums of
\begin{align*}
\exp\left(-\frac{ih}{2}\sum_{j=1}^n\partial_{x_j}\partial_{\xi_j}\right)a.
\end{align*}
Indeed,
\begin{align*}
\sum_{k=0}^{N-1}\frac{1}{k!}\left(-\frac{ih}{2}\right)^k\left(\sum_{j=1}^n\partial_{x_j}\partial_{\xi_j}\right)^k a=\sum_{|\alpha|<N}\frac{1}{\alpha!}\left(-\frac{ih}{2}\right)^{|\alpha|}\partial_x^\alpha\partial_\xi^\alpha a.
\end{align*}
Since each $\xi$-derivative lowers the order of a symbol by one and $x$-derivatives preserve the order in $S^m_{1,0}$, the term with $|\alpha|=k$ belongs to $h^kS^{m-k}_{1,0}(T^*\mathbb{R}^n)$.
[/step]
[step:Apply amplitude reduction once to the full Weyl amplitude]
Let $A:\mathbb{R}^n_x\times\mathbb{R}^n_y\times\mathbb{R}^n_\xi\times(0,h_0]\to\mathbb{C}$ denote the full Weyl amplitude
\begin{align*}
A(x,y,\xi;h)=a\left(\frac{x+y}{2},\xi;h\right).
\end{align*}
For later comparison with the Taylor expansion, also let $A_N:\mathbb{R}^n_x\times\mathbb{R}^n_y\times\mathbb{R}^n_\xi\times(0,h_0]\to\mathbb{C}$ denote the remainder amplitude $R_N$. Here the word amplitude means a smooth function of both endpoint variables $x$ and $y$ and the covariable $\xi$, rather than a left symbol depending only on $x$ and $\xi$. The explicit formula for $R_N$ shows that $A_N$ is a finite sum of terms with a factor $(y-x)^\alpha$ for $|\alpha|=N$ and a symbol factor of order $m$ differentiated $N$ times in $x$. Converting each factor $y_j-x_j$ to $ih\partial_{\xi_j}$ and integrating by parts $N$ times in $\xi$ produces an oscillatory integral whose remaining amplitude is still allowed to depend on $y$ through $x+t(y-x)/2$.
We use the following standard external theorem from the semiclassical pseudodifferential calculus, stated in the precise form needed here. Suppose $A(x,y,\xi;h)$ is a smooth amplitude satisfying, for every multi-indices $\beta,\delta,\gamma\in\mathbb{N}_0^n$, an estimate
\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 $0<h\le h_0$. Then its oscillatory integral is equal, as an operator $\mathcal{S}(\mathbb{R}^n)\to\mathcal{S}'(\mathbb{R}^n)$, to $\operatorname{Op}_h(\widetilde b)+\operatorname{Op}_h(c)$ for a single left symbol $\widetilde b\in S^m_{1,0}(T^*\mathbb{R}^n)$ and a single residual left symbol $c$ in the sense stated in the theorem. Moreover $\widetilde b$ has the asymptotic expansion
\begin{align*}
\widetilde b(x,\xi;h)\sim \sum_\alpha \frac{1}{\alpha!}\left(-ih\right)^{|\alpha|}\partial_\xi^\alpha\partial_y^\alpha A(x,y,\xi;h)\big|_{y=x},
\end{align*}
meaning that the difference between $\widetilde b$ and the sum over $|\alpha|<N$ lies in $h^N S^{m-N}_{1,0}(T^*\mathbb{R}^n)$ for every $N\in\mathbb{N}$. The residual symbol $c$ is part of this single reduction theorem; it is not obtained by choosing a different remainder for each $N$. The theorem is proved by inserting a compact cutoff near the diagonal $x=y$, Taylor expanding in $y-x$, integrating by parts in $\xi$, and using nonstationary phase away from the diagonal to produce the residual symbol.
We verify these hypotheses for the full Weyl amplitude $A$. For every multi-index triple $\beta,\delta,\gamma\in\mathbb{N}_0^n$, differentiating the expression $(x+y)/2$ only produces the bounded coefficients $1/2$. The symbol estimates for $a$ therefore give constants $C_{\beta\delta\gamma}>0$, depending only on the corresponding $S^m_{1,0}$ seminorms of $a$, such that
\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 $0<h\le h_0$. Hence the reduction theorem gives a single left Kohn--Nirenberg symbol $\widetilde b\in S^m_{1,0}(T^*\mathbb{R}^n)$ and a single residual left symbol $c$ such that
\begin{align*}
\operatorname{Op}_h^w(a)=\operatorname{Op}_h(\widetilde b)+\operatorname{Op}_h(c).
\end{align*}
The expansion supplied by the same theorem is computed by differentiating $A$ in the $y$ variable at $y=x$. Since
\begin{align*}
\partial_y^\alpha A(x,y,\xi;h)\big|_{y=x}=2^{-|\alpha|}\partial_x^\alpha a(x,\xi;h),
\end{align*}
we obtain, 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]
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]
[/step]
[step:Name the reduced symbol and rewrite its expansion in exponential form]
Let $b:=\widetilde b$, where $\widetilde b$ is the single left symbol produced by the amplitude-to-left-symbol reduction theorem in the previous step. That theorem already gives a single residual left symbol $c$ such that
\begin{align*}
\operatorname{Op}_h^w(a)=\operatorname{Op}_h(b)+\operatorname{Op}_h(c)
\end{align*}
as operators $\mathcal{S}(\mathbb{R}^n)\to\mathcal{S}'(\mathbb{R}^n)$.
The multi-index expansion from the previous step is the same as the exponential expansion. Indeed, for each $k\in\mathbb{N}_0$,
\begin{align*}
\left(\sum_{j=1}^n\partial_{x_j}\partial_{\xi_j}\right)^k a=\sum_{|\alpha|=k}\frac{k!}{\alpha!}\partial_x^\alpha\partial_\xi^\alpha a.
\end{align*}
Therefore, for every $N\in\mathbb{N}$,
\begin{align*}
b-\sum_{k=0}^{N-1}\frac{1}{k!}\left(-\frac{ih}{2}\right)^k\left(\sum_{j=1}^n\partial_{x_j}\partial_{\xi_j}\right)^k a\in h^N S^{m-N}_{1,0}(T^*\mathbb{R}^n).
\end{align*}
This is precisely the asserted asymptotic expansion of $b$.
[/step]
[step:Read off the principal symbol]
Taking the first two terms of the expansion gives
\begin{align*}
b=a-\frac{ih}{2}\sum_{j=1}^n\partial_{x_j}\partial_{\xi_j}a+O_{S^{m-2}_{1,0}}(h^2).
\end{align*}
The difference $b-a$ lies in $hS^{m-1}_{1,0}(T^*\mathbb{R}^n)$, so $a$ and $b$ define the same principal symbol in the quotient $S^m_{1,0}/hS^{m-1}_{1,0}$. Therefore the Weyl and Kohn--Nirenberg quantizations have the same principal symbol. This completes the proof.
[/step]