[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*}[/step]

[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: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]