[proofplan]
We compute the composition kernel and rewrite it as the left quantization of an oscillatory integral symbol. Near the diagonal in the auxiliary variables, Taylor expansion of $b$ in the base variable converts powers of displacement into $h$ times derivatives of $a$ by integration by parts in the dual variable. Away from the diagonal, a nonstationary phase operator gives rapid decay in $h$, producing only a residual contribution. The Taylor remainder is estimated by the standard symbol seminorms and gives exactly an $h^N S^{m+m'-N}$ remainder.
[/proofplan]
[step:Write the composed kernel as an oscillatory integral symbol]
Fix $0<h\le h_0$ and, inside this proof only, write $a$ and $b$ for the fixed symbols $a_h$ and $b_h$. For $u \in \mathcal{S}(\mathbb{R}^n)$, the composed operator $\operatorname{Op}_h(a)\operatorname{Op}_h(b)u$ is the regularized fourfold oscillatory integral specified in the theorem statement. With the standard admissible Schwartz cutoffs inserted in all integration variables, the cutoff-regularized composition is
\begin{align*}\operatorname{Op}_h(a)\operatorname{Op}_h(b)u(x)=(2\pi h)^{-2n}\operatorname{Os}\!\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}e^{i(x-y)\cdot \eta/h}e^{i(y-z)\cdot \xi/h}a_h(x,\eta)b_h(y,\xi)u(z)\,d\mathcal{L}^n(z)\,d\mathcal{L}^n(\xi)\,d\mathcal{L}^n(y)\,d\mathcal{L}^n(\eta).\end{align*}
Define the formal composition symbol $c_h:T^*\mathbb{R}^n \to \mathbb{C}$ by the regularized oscillatory integral
\begin{align*}
c_h(x,\xi)=(2\pi h)^{-n}\operatorname{Os}\!\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}e^{-iY\cdot \theta/h}a_h(x,\xi+\theta)b_h(x+Y,\xi)\,d\mathcal{L}^n(Y)\,d\mathcal{L}^n(\theta).
\end{align*}
Here the regularization is the cutoff limit from the theorem statement, applied first with Schwartz cutoffs in $(y,\eta,z,\xi)$ and then passed to the standard oscillatory-integral limit. Because $u\in\mathcal{S}(\mathbb{R}^n)$ and the cutoff factors are Schwartz while $a$ and $b$ have polynomial growth in the covariables, each cutoff integral is absolutely convergent; hence Fubini's theorem applies at the cutoff level. The oscillatory-integral continuity assumption in the theorem statement permits passing the cutoff identity below to the regularized limit.
Define the auxiliary maps $Y:\mathbb{R}^n_y\times\mathbb{R}^n_x\to\mathbb{R}^n$ and $\theta:\mathbb{R}^n_\eta\times\mathbb{R}^n_\xi\to\mathbb{R}^n$ by $Y(y,x)=y-x$ and $\theta(\eta,\xi)=\eta-\xi$. For fixed $(x,\xi)$, the change of variables
the affine map $(y,\eta) \mapsto (Y,\theta)=(y-x,\eta-\xi)$
has Jacobian determinant $1$, so it preserves the product measure
$d\mathcal{L}^n(y)\,d\mathcal{L}^n(\eta)$.
Substituting $\eta=\xi+\theta$ and $y=x+Y$ gives
\begin{align*}
e^{i(x-y)\cdot \eta/h}e^{i(y-z)\cdot \xi/h}
=
e^{i(x-z)\cdot \xi/h}e^{-iY\cdot \theta/h}.
\end{align*}
Therefore
\begin{align*}\operatorname{Op}_h(a)\operatorname{Op}_h(b)u(x)=(2\pi h)^{-n}\operatorname{Os}\!\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}e^{i(x-z)\cdot \xi/h}c_h(x,\xi)u(z)\,d\mathcal{L}^n(z)\,d\mathcal{L}^n(\xi).\end{align*}
after passing from the compactly supported cutoff formula to the oscillatory limit. Thus the composition is $\operatorname{Op}_h(c_h)$ with the oscillatory-integral interpretation fixed in the theorem statement.
[/step]
[step:Extract the Taylor polynomial terms]
Fix $N \in \mathbb{N}$. For each fixed $(x,\xi)\in T^*\mathbb{R}^n$, apply Taylor's formula to the map
\begin{align*}
Y \mapsto b(x+Y,\xi)
\end{align*}
at $Y=0$ through order $N-1$. This gives
\begin{align*}
b(x+Y,\xi)
=
\sum_{|\alpha|<N}
\frac{Y^\alpha}{\alpha!}
(\partial_x^\alpha b)(x,\xi)
+
R_N(x,Y,\xi),
\end{align*}
where the integral remainder is
\begin{align*}
R_N(x,Y,\xi)
=
\sum_{|\alpha|=N}
\frac{N Y^\alpha}{\alpha!}
\int_0^1
(1-t)^{N-1}
(\partial_x^\alpha b)(x+tY,\xi)\,
d\mathcal{L}^1(t).
\end{align*}
For a multi-index $\alpha \in \mathbb{N}_0^n$, the identity
\begin{align*}
Y^\alpha e^{-iY\cdot \theta/h}
=
(ih)^{|\alpha|}\partial_\theta^\alpha e^{-iY\cdot \theta/h}
\end{align*}
holds. Integrating by parts in $\theta$ transfers $\partial_\theta^\alpha$ from the exponential to $a(x,\xi+\theta)$ and gives the sign factor $(-1)^{|\alpha|}$. Since
\begin{align*}
(-1)^{|\alpha|}i^{|\alpha|}
=
i^{-|\alpha|},
\end{align*}
the Taylor polynomial contributes
\begin{align*}
\frac{h^{|\alpha|}}{i^{|\alpha|}\alpha!}
(\partial_\xi^\alpha a)(x,\xi)
(\partial_x^\alpha b)(x,\xi).
\end{align*}
[guided]
We now explain why Taylor expansion in the $Y$ variable produces precisely the displayed coefficients. The oscillatory integral for $c_h$ has the phase $-Y\cdot \theta/h$. The variable $Y$ measures the displacement in the base point of $b$, while $\theta$ measures the displacement in the covariable of $a$. Thus Taylor expansion of $b(x+Y,\xi)$ in powers of $Y$ should become differentiation of $a(x,\xi+\theta)$ in $\theta$, hence in $\xi$.
Fix $N \in \mathbb{N}$. Taylor's formula with integral remainder for the smooth map $Y \mapsto b(x+Y,\xi)$ gives
\begin{align*}b(x+Y,\xi)=\sum_{|\alpha|<N}\frac{Y^\alpha}{\alpha!}(\partial_x^\alpha b)(x,\xi)+R_N(x,Y,\xi).\end{align*}
where
\begin{align*}R_N(x,Y,\xi)=\sum_{|\alpha|=N}\frac{N Y^\alpha}{\alpha!}\int_0^1(1-t)^{N-1}(\partial_x^\alpha b)(x+tY,\xi)\,d\mathcal{L}^1(t).\end{align*}
Consider one Taylor monomial. In the oscillatory factor, differentiation with respect to $\theta_j$ gives
\begin{align*}\partial_{\theta_j}e^{-iY\cdot \theta/h}=-\frac{iY_j}{h}e^{-iY\cdot \theta/h}.\end{align*}
Equivalently,
\begin{align*}Y_j e^{-iY\cdot \theta/h}=ih\,\partial_{\theta_j}e^{-iY\cdot \theta/h}.\end{align*}
Iterating this identity gives
\begin{align*}Y^\alpha e^{-iY\cdot \theta/h}=(ih)^{|\alpha|}\partial_\theta^\alpha e^{-iY\cdot \theta/h}.\end{align*}
We then integrate by parts in the $\theta$ variable in the regularized oscillatory-integral sense fixed in the theorem statement. The cutoff-limit continuity and integration-by-parts hypotheses ensure that no boundary term remains after the cutoff limit is taken. The derivative moves from the exponential to $a(x,\xi+\theta)$, producing $(-1)^{|\alpha|}$. Since $\partial_\theta^\alpha a(x,\xi+\theta)=\partial_\xi^\alpha a(x,\xi+\theta)$, the monomial contributes
\begin{align*}(2\pi h)^{-n}\frac{(ih)^{|\alpha|}(-1)^{|\alpha|}}{\alpha!}\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}e^{-iY\cdot \theta/h}(\partial_\xi^\alpha a)(x,\xi+\theta)(\partial_x^\alpha b)(x,\xi)\,d\mathcal{L}^n(Y)\,d\mathcal{L}^n(\theta).\end{align*}
The remaining oscillatory integral is the semiclassical Fourier inversion identity at the origin, so it evaluates the $\theta$-dependent factor at $\theta=0$. Therefore the contribution is
\begin{align*}\frac{h^{|\alpha|}}{i^{|\alpha|}\alpha!}(\partial_\xi^\alpha a)(x,\xi)(\partial_x^\alpha b)(x,\xi).\end{align*}
This is the desired coefficient.
[/guided]
[/step]
[step:Estimate the Taylor remainder in the expected symbol class]
Let $r_N:T^*\mathbb{R}^n \to \mathbb{C}$ denote the contribution of $R_N$ to $c_h$:
\begin{align*}
r_N(x,\xi)=(2\pi h)^{-n}\operatorname{Os}\!\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}e^{-iY\cdot \theta/h}a(x,\xi+\theta)R_N(x,Y,\xi)\,d\mathcal{L}^n(Y)\,d\mathcal{L}^n(\theta).
\end{align*}
For $|\alpha|=N$, write the corresponding summand of $R_N$ as $Y^\alpha q_\alpha(x,Y,\xi)$, where $q_\alpha:\mathbb{R}^n_x\times\mathbb{R}^n_Y\times\mathbb{R}^n_\xi\to\mathbb{C}$ is the map defined by
\begin{align*}
q_\alpha(x,Y,\xi)=\frac{N}{\alpha!}\int_0^1(1-t)^{N-1}(\partial_x^\alpha b)(x+tY,\xi)\,d\mathcal{L}^1(t).
\end{align*}
The map $q_\alpha:\mathbb{R}^n_x\times\mathbb{R}^n_Y\times\mathbb{R}^n_\xi\to\mathbb{C}$ is a symbol of order $m'$ in $\xi$, uniformly for $Y$ on compact sets and with polynomial control under $Y$-derivatives, because $b\in S^{m'}$ and differentiation in $Y$ differentiates $b$ in $x$ with a factor $t\in[0,1]$.
We estimate $r_N$ directly as a regularized oscillatory integral, without inserting a cutoff that would create lower-order Leibniz terms. For each $|\alpha|=N$, use
\begin{align*}
Y^\alpha e^{-iY\cdot\theta/h}=(ih)^N\partial_\theta^\alpha e^{-iY\cdot\theta/h}
\end{align*}
and integrate by parts in $\theta$ in the regularized oscillatory-integral sense. Since $q_\alpha(x,Y,\xi)$ is independent of $\theta$, every $\theta$-derivative falls on $a_h(x,\xi+\theta)$, and therefore
\begin{align*}
r_N(x,\xi)=h^N\sum_{|\alpha|=N}\operatorname{Os}\!\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}e^{-iY\cdot\theta/h}A_\alpha(x,Y,\theta,\xi)\,d\mathcal{L}^n(Y)\,d\mathcal{L}^n(\theta),
\end{align*}
where the amplitude $A_\alpha:\mathbb{R}^n_x\times\mathbb{R}^n_Y\times\mathbb{R}^n_\theta\times\mathbb{R}^n_\xi\to\mathbb{C}$ is defined by
\begin{align*}
A_\alpha(x,Y,\theta,\xi)=i^N(-1)^N(\partial_\xi^\alpha a_h)(x,\xi+\theta)q_\alpha(x,Y,\xi).
\end{align*}
The integration by parts has no boundary term because it is performed after inserting Schwartz cutoffs and then passing to the oscillatory-integral limit. Differentiating $A_\alpha$ in $x$, $Y$, $\theta$, and $\xi$ gives finite sums of products of derivatives of $\partial_\xi^\alpha a_h$ and derivatives of $q_\alpha$. The defining symbol estimates give
\begin{align*}
|\partial_x^{\beta_1}\partial_\theta^{\delta}\partial_\xi^{\gamma_1}(\partial_\xi^\alpha a_h)(x,\xi+\theta)|\le C\langle\xi+\theta\rangle^{m-N-|\delta|-|\gamma_1|}
\end{align*}
and
\begin{align*}
|\partial_x^{\beta_2}\partial_Y^{\mu}\partial_\xi^{\gamma_2}q_\alpha(x,Y,\xi)|\le C\langle\xi\rangle^{m'-|\gamma_2|},
\end{align*}
with constants depending only on finitely many symbol seminorms of $a_h$ and $b_h$. The bilinear oscillatory-product estimate assumed in the theorem statement applies to the amplitude $A_\alpha$. Its hypotheses are met because the displayed bounds give symbolic order $m-N$ in the shifted covariable factor $\xi+\theta$, symbolic order $m'$ in the unshifted covariable factor $\xi$, and polynomial control in the auxiliary variable $Y$. For fixed $\beta,\gamma\in\mathbb{N}_0^n$ and $N\in\mathbb{N}$, that assumed estimate supplies integers $J=J(\beta,\gamma,N,n,m,m')$ and $P=P(\beta,\gamma,N,n,m,m')$ such that the resulting seminorm is controlled by the finitely many seminorms
\begin{align*}
\sup_{(x,\zeta,h)}\langle \zeta\rangle^{-m+N+|\nu|}|\partial_x^\mu\partial_\zeta^\nu\partial_\xi^\alpha a_h(x,\zeta)|
\end{align*}
with $|\mu|+|\nu|\le J$, together with
\begin{align*}
\sup_{(x,Y,\xi,h)}\langle Y\rangle^{-P}\langle \xi\rangle^{-m'+|\nu|}|\partial_x^\mu\partial_Y^\lambda\partial_\xi^\nu q_{\alpha,h}(x,Y,\xi)|
\end{align*}
with $|\mu|+|\lambda|+|\nu|\le J$. These quantities are finite by the symbol estimates for $a_h$ and by the displayed estimates for $q_{\alpha,h}$. Hence there is a constant $C_{\beta\gamma N}>0$, depending only on these finitely many seminorms and on the fixed parameters of the assumed oscillatory-product estimate, such that
\begin{align*}
|\partial_x^\beta\partial_\xi^\gamma r_N(x,\xi)|\le C_{\beta\gamma N}h^N\langle\xi\rangle^{m+m'-N-|\gamma|}.
\end{align*}
Thus $r_N\in h^N S^{m+m'-N}(T^*\mathbb{R}^n)$.
[/step]
[step:Assemble the asymptotic product symbol]
Define $a\# b$ to be the symbol family $c=(c_h)_{0<h\le h_0}$ constructed by the oscillatory integral above. Its asymptotic expansion is
\begin{align*}a \# b\sim\sum_{\alpha \in \mathbb{N}_0^n}\frac{h^{|\alpha|}}{i^{|\alpha|}\alpha!}(\partial_\xi^\alpha a_h)(x,\xi)(\partial_x^\alpha b_h)(x,\xi).\end{align*}
The previous Taylor-remainder estimate holds for every $N\in\mathbb{N}$ and is an estimate for the actual oscillatory symbol $c_h$. Hence $a\# b=c$ satisfies each finite-order remainder estimate directly. If $c'$ is another symbol family with the same complete asymptotic expansion, then for every $N\in\mathbb{N}$ the difference $c-c'$ lies in $h^N S^{m+m'-N}(T^*\mathbb{R}^n)$; since $N$ is arbitrary, $c-c'\in h^\infty S^{-\infty}(T^*\mathbb{R}^n)$. Quantizing an element of $h^\infty S^{-\infty}$ gives a residual Schwartz kernel by repeated integration by parts in the defining kernel integral and the residual definition in the theorem statement. The previous steps show that, for every $N\in\mathbb{N}$,
\begin{align*}(a \# b)_h-\sum_{|\alpha|<N}\frac{h^{|\alpha|}}{i^{|\alpha|}\alpha!}(\partial_\xi^\alpha a_h)(\partial_x^\alpha b_h)\in h^N S^{m+m'-N}(T^*\mathbb{R}^n).\end{align*}
In particular the leading term belongs to $S^{m+m'}$, and all subsequent terms have lower symbolic order, so
\begin{align*}a\# b \in S^{m+m'}(T^*\mathbb{R}^n).\end{align*}
Combining this with the kernel computation gives
\begin{align*}\operatorname{Op}_h(a)\operatorname{Op}_h(b)=\operatorname{Op}_h(a\# b)\end{align*}
modulo a residual operator. This proves the stated composition formula and its full asymptotic expansion.
[/step]