[proofplan]
The proof is local, because both the smoothing symbol seminorms and residual kernel estimates are tested in coordinate charts on compact sets. For the forward implication, we differentiate the semiclassical oscillatory kernel and bound the resulting Fourier integral using sufficiently strong rapidly vanishing $S^{-\infty}$ seminorms; the factor $(2\pi h)^{-n}$ and the derivatives in $y$ only cost finitely many powers of $h^{-1}$. For the converse, we recover a local symbol from a residual kernel by the semiclassical Fourier transform in the difference variable $z=x-y$, and integration by parts in $z$ gives rapid decay in $\xi$. Proper support and compact localization ensure that all constants are uniform on the compact coordinate sets under consideration.
[/proofplan]
[step:Reduce the assertion to one coordinate patch with compact support]
Fix a coordinate patch $U \subset M$ identified with an open subset of $\mathbb{R}^n$, and let $K \subset U \times U$ be an arbitrary compact set. Choose compact sets $K_x,K_y \subset U$ with $K \subset K_x \times K_y$, for instance compact neighbourhoods of the coordinate projections of $K$. It is enough to prove the result for local kernels over $K_x \times K_y$, because a properly supported global quantization is obtained by summing finitely many such coordinate-local pieces on every compact subset of $M \times M$.
Let
\begin{align*}
a: U \times \mathbb{R}^n \times (0,h_0] \to \mathbb{C}
\end{align*}
denote the local representative of the symbol. The local operator
\begin{align*}
A_h:C_c^\infty(U) \to C^\infty(U)
\end{align*}
is defined by
\begin{align*}
(A_h u)(x)=(2\pi h)^{-n}\int_U\int_{\mathbb{R}^n}e^{i(x-y)\cdot \xi/h}a(x,\xi;h)u(y)\,d\mathcal{L}^n(\xi)\,d\mathcal{L}^n(y).
\end{align*}
Its local Schwartz kernel is the smooth function
\begin{align*}
K_a:U\times U\times (0,h_0]\to \mathbb{C}
\end{align*}
given by
\begin{align*}
K_a(x,y;h)=(2\pi h)^{-n}\int_{\mathbb{R}^n}e^{i(x-y)\cdot \xi/h}a(x,\xi;h)\,d\mathcal{L}^n(\xi).
\end{align*}
[/step]
[step:Differentiate the oscillatory kernel and absorb all powers of $h^{-1}$]
Let $\alpha,\beta \in \mathbb{N}_0^n$ be fixed. Differentiation under the integral sign is justified by dominated convergence applied on $K_x\times K_y$: after any derivative of total order at most $|\alpha|+|\beta|$, the integrand is dominated by a finite sum of functions of the form $C h^{-q}\langle \xi\rangle^{-m+q}$ with $m>n+q$, which are integrable in $\xi$ with respect to $\mathcal{L}^n$. Thus every derivative $\partial_x^\alpha \partial_y^\beta K_a$ is a finite sum of terms of the form
\begin{align*}
(2\pi h)^{-n}h^{-|\gamma|-|\delta|}\int_{\mathbb{R}^n}e^{i(x-y)\cdot \xi/h}\xi^\gamma \partial_x^\delta a(x,\xi;h)\,d\mathcal{L}^n(\xi),
\end{align*}
where $\gamma,\delta \in \mathbb{N}_0^n$ satisfy $|\gamma|+|\delta|\leq |\alpha|+|\beta|$. The multi-index $\gamma$ records derivatives landing on the phase, while $\delta$ records derivatives landing on $a$.
Fix $N\in\mathbb{N}$. Choose an integer $m>n+\max |\gamma|$, where the maximum is taken over the finitely many terms above. For each such term, the smoothing estimate for $a$ with exponent
\begin{align*}
N_1=N+n+|\alpha|+|\beta|
\end{align*}
gives
\begin{align*}
|\xi^\gamma \partial_x^\delta a(x,\xi;h)|\leq C_{\delta,0,m,N_1,K_x} h^{N_1}\langle \xi\rangle^{-m+|\gamma|}.
\end{align*}
Since $m-|\gamma|>n$, the function $\xi \mapsto \langle \xi\rangle^{-m+|\gamma|}$ is integrable over $\mathbb{R}^n$ with respect to $\mathcal{L}^n$. Therefore each differentiated kernel term is bounded on $K_x\times K_y$ by
\begin{align*}
(2\pi)^{-n}C_{\delta,0,m,N_1,K_x}h^{-n-|\gamma|-|\delta|}h^{N_1}\int_{\mathbb{R}^n}\langle \xi\rangle^{-m+|\gamma|}\,d\mathcal{L}^n(\xi).
\end{align*}
Because $|\gamma|+|\delta|\leq |\alpha|+|\beta|$, this is $O(h^N)$. Summing the finitely many differentiated terms and absorbing the finite Leibniz combinatorial coefficients into the constant gives
\begin{align*}
\sup_{(x,y)\in K_x\times K_y}|\partial_x^\alpha\partial_y^\beta K_a(x,y;h)|\leq C_{\alpha,\beta,N,K_x,K_y}h^N.
\end{align*}
[guided]
The only possible difficulty in this direction is the prefactor $(2\pi h)^{-n}$ and the powers of $h^{-1}$ produced when derivatives hit the phase. We handle them by asking for more decay in the symbol estimate than the final power $h^N$ we want.
Fix derivative orders $\alpha,\beta \in \mathbb{N}_0^n$. When $\partial_{x_j}$ hits the exponential factor $e^{i(x-y)\cdot \xi/h}$, it contributes $i\xi_j/h$. When $\partial_{y_j}$ hits the same exponential factor, it contributes $-i\xi_j/h$. When $\partial_{x_j}$ hits the symbol, it contributes $\partial_{x_j}a$. Thus $\partial_x^\alpha\partial_y^\beta K_a$ is a finite sum of integrals of the form
\begin{align*}
(2\pi h)^{-n}h^{-|\gamma|-|\delta|}\int_{\mathbb{R}^n}e^{i(x-y)\cdot \xi/h}\xi^\gamma\partial_x^\delta a(x,\xi;h)\,d\mathcal{L}^n(\xi),
\end{align*}
with $|\gamma|+|\delta|\leq |\alpha|+|\beta|$.
Now choose $m>n+\max|\gamma|$. This choice is made so that the majorant $\langle \xi\rangle^{-m+|\gamma|}$ is integrable over $\mathbb{R}^n$. Since $a$ is rapidly vanishing in every $S^{-\infty}$ seminorm, we may apply the symbol estimate with the larger exponent
\begin{align*}
N_1=N+n+|\alpha|+|\beta|.
\end{align*}
It gives, uniformly for $x\in K_x$ and $\xi\in\mathbb{R}^n$,
\begin{align*}
|\xi^\gamma\partial_x^\delta a(x,\xi;h)|\leq C_{\delta,0,m,N_1,K_x}h^{N_1}\langle \xi\rangle^{-m+|\gamma|}.
\end{align*}
Substituting this estimate into the differentiated kernel integral gives
\begin{align*}
|\partial_x^\alpha\partial_y^\beta K_a(x,y;h)|\leq C h^{-n-|\gamma|-|\delta|}h^{N_1}\int_{\mathbb{R}^n}\langle \xi\rangle^{-m+|\gamma|}\,d\mathcal{L}^n(\xi)
\end{align*}
for a constant $C$ independent of $h$. Since $|\gamma|+|\delta|\leq|\alpha|+|\beta|$, the exponent of $h$ is at least
\begin{align*}
N_1-n-|\alpha|-|\beta|=N.
\end{align*}
This proves the desired $O(h^N)$ bound for every differentiated kernel seminorm.
[/guided]
[/step]
[step:Conclude that rapidly vanishing smoothing symbols produce residual operators]
The previous step proves the defining residual estimate for every compact set in one coordinate patch. A properly supported global quantization is locally a finite sum of such kernels, multiplied by smooth cutoff functions whose derivatives are independent of $h$. Differentiating those cutoffs only changes the constants and does not change the power of $h$. Hence $\operatorname{Op}_h(a)$ is residual.
[/step]
[step:Recover a local symbol from a properly supported residual kernel]
Conversely, let $R_h$ be a properly supported residual operator. In a coordinate patch $U\subset\mathbb{R}^n$, let
\begin{align*}
K_R:U\times U\times(0,h_0]\to\mathbb{C}
\end{align*}
be its local smooth kernel. Choose cutoff functions $\chi,\psi \in C_c^\infty(U)$ with $\chi=1$ on $K_x$ and with $\psi=1$ on the $y$-projection of the properly supported kernel over $\operatorname{supp}\chi$. Replacing $K_R$ by $\chi(x)K_R(x,y;h)\psi(y)$, assume that $K_R(x,y;h)$ is supported in a compact subset of $K_x\times K_y\subset U\times U$.
Define the difference variable $z=x-y$. For each $x\in K_x$ and $h\in(0,h_0]$, define
\begin{align*}
F_x:\mathbb{R}^n\to\mathbb{C}
\end{align*}
by $F_x(z;h)=K_R(x,x-z;h)$ on the set where $x-z\in U$ and by $F_x(z;h)=0$ outside the compact $z$-projection of the localized support. Because the cutoff-localized kernel is smooth and vanishes in a neighbourhood of the boundary of this $z$-support, $F_x(\cdot;h)$ is a compactly supported smooth function of $z$, uniformly for $x\in K_x$.
Define the local recovered symbol
\begin{align*}
r:K_x\times\mathbb{R}^n\times(0,h_0]&\to\mathbb{C}
\end{align*}
by
\begin{align*}
r(x,\xi;h)=\int_{\mathbb{R}^n}e^{-iz\cdot \xi/h}F_x(z;h)\,d\mathcal{L}^n(z).
\end{align*}
The semiclassical Fourier inversion formula for compactly supported smooth functions, with the normalization used in the definition of $\operatorname{Op}_h$, gives
\begin{align*}
F_x(z;h)=(2\pi h)^{-n}\int_{\mathbb{R}^n}e^{iz\cdot\xi/h}r(x,\xi;h)\,d\mathcal{L}^n(\xi).
\end{align*}
Substituting $z=x-y$ yields
\begin{align*}
K_R(x,y;h)=(2\pi h)^{-n}\int_{\mathbb{R}^n}e^{i(x-y)\cdot\xi/h}r(x,\xi;h)\,d\mathcal{L}^n(\xi)
\end{align*}
as an identity of smooth kernels on the localized patch.
[/step]
[step:Prove rapid smoothing seminorm bounds for the recovered symbol]
Fix multi-indices $\alpha,\beta\in\mathbb{N}_0^n$, an integer $m\geq0$, and $N\in\mathbb{N}$. Differentiating the formula for $r$ gives a finite sum of terms of the form
\begin{align*}
h^{-|\beta_1|}\int_{\mathbb{R}^n}e^{-iz\cdot\xi/h}z^{\beta_1}\partial_x^{\alpha_1}\partial_y^{\alpha_2}K_R(x,x-z;h)\,d\mathcal{L}^n(z),
\end{align*}
where $\beta_1\leq \beta$ and $|\alpha_1|+|\alpha_2|\leq|\alpha|$. The residual hypothesis gives $O(h^{N_2})$ bounds for every differentiated kernel seminorm, with $N_2$ to be chosen.
To gain decay in $\xi$, use the identity
\begin{align*}
(1+|\xi|^2)^q e^{-iz\cdot\xi/h}=(1-h^2\Delta_z)^q e^{-iz\cdot\xi/h}
\end{align*}
for an integer $q\geq m/2$. Integrating by parts in $z$ moves $(1-h^2\Delta_z)^q$ onto the compactly supported smooth amplitude
\begin{align*}
z\mapsto z^{\beta_1}\partial_x^{\alpha_1}\partial_y^{\alpha_2}K_R(x,x-z;h).
\end{align*}
No boundary term appears because this amplitude is compactly supported in $z$. Each $z$-derivative is a finite linear combination of $y$-derivatives of $K_R$ multiplied by bounded powers of $z$, and each factor $h^2\Delta_z$ only improves the power of $h$. Therefore the residual kernel estimates imply
\begin{align*}
\langle\xi\rangle^m|\partial_x^\alpha\partial_\xi^\beta r(x,\xi;h)|\leq C_{\alpha,\beta,m,N,K_x}h^N
\end{align*}
after choosing $N_2$ larger than $N+|\beta|$. Thus $r$ is a rapidly vanishing $S^{-\infty}$ symbol.
[guided]
The recovered symbol is a semiclassical Fourier transform in the difference variable $z=x-y$. For each fixed $x\in K_x$ and $h\in(0,h_0]$, the localized amplitude $F_x:z\mapsto K_R(x,x-z;h)$ is a compactly supported smooth function of $z$, after extension by zero outside the compact $z$-projection of the cutoff-localized support. Therefore the symbol is
\begin{align*}
r(x,\xi;h)=\int_{\mathbb{R}^n}e^{-iz\cdot\xi/h}F_x(z;h)\,d\mathcal{L}^n(z),
\end{align*}
and semiclassical Fourier inversion gives
\begin{align*}
K_R(x,y;h)=(2\pi h)^{-n}\int_{\mathbb{R}^n}e^{i(x-y)\cdot\xi/h}r(x,\xi;h)\,d\mathcal{L}^n(\xi)
\end{align*}
on the localized patch.
Differentiating in $\xi$ is the only operation that creates negative powers of $h$: every $\partial_{\xi_j}$ brings down the factor $-iz_j/h$. Thus, after applying $\partial_x^\alpha\partial_\xi^\beta$, each term has the form
\begin{align*}
h^{-|\beta_1|}\int_{\mathbb{R}^n}e^{-iz\cdot\xi/h}z^{\beta_1}\partial_x^{\alpha_1}\partial_y^{\alpha_2}K_R(x,x-z;h)\,d\mathcal{L}^n(z).
\end{align*}
The support in $z$ is contained in a fixed compact set because the kernel has been localized and the operator is properly supported.
We still need decay in $\xi$. The useful identity is
\begin{align*}
(1+|\xi|^2)^q e^{-iz\cdot\xi/h}=(1-h^2\Delta_z)^q e^{-iz\cdot\xi/h}.
\end{align*}
This identity follows because $\Delta_z e^{-iz\cdot\xi/h}=-|\xi|^2h^{-2}e^{-iz\cdot\xi/h}$. Therefore applying $1-h^2\Delta_z$ to the exponential multiplies it by $1+|\xi|^2$.
Choose $q\geq m/2$. Multiplying by $\langle\xi\rangle^m$ is then controlled by multiplying by $(1+|\xi|^2)^q$. We integrate by parts in $z$ and move $(1-h^2\Delta_z)^q$ from the exponential onto the amplitude
\begin{align*}
z\mapsto z^{\beta_1}\partial_x^{\alpha_1}\partial_y^{\alpha_2}K_R(x,x-z;h).
\end{align*}
The compact support eliminates boundary terms. Derivatives in $z$ become derivatives in the second kernel variable $y$ together with derivatives of the polynomial factor $z^{\beta_1}$. Since $z$ stays in a fixed compact set, all powers of $z$ are uniformly bounded. The residual hypothesis says that every such differentiated kernel is $O(h^{N_2})$ for any prescribed $N_2$. Choosing $N_2>N+|\beta|$ absorbs the prefactor $h^{-|\beta_1|}$, and the additional factors $h^2$ from $(1-h^2\Delta_z)^q$ do not hurt the estimate. Hence
\begin{align*}
\langle\xi\rangle^m|\partial_x^\alpha\partial_\xi^\beta r(x,\xi;h)|\leq C_{\alpha,\beta,m,N,K_x}h^N.
\end{align*}
This is exactly the required rapidly vanishing smoothing symbol seminorm bound. Together with the Fourier inversion identity above, it proves the localized converse: the cutoff-localized residual kernel is represented by a rapidly vanishing smoothing symbol in the fixed left-quantization convention. The remaining global step is not an ordinary coordinate-free partition of a single formula; it is the locally finite sum of these fixed coordinate-local quantizations, and all terms not captured on a chosen coordinate patch are residual because their cutoff-localized kernels satisfy the same $O(h^N)$ smooth kernel bounds.
[/guided]
[/step]
[step:Patch the recovered local symbols and identify the operator modulo residual equality]
Choose the locally finite coordinate cover, partition of unity, and left-quantization cutoffs fixed in the statement. For each coordinate index $j$, let $\chi_j\in C_c^\infty(U_j)$ be the output cutoff and let $\psi_j\in C_c^\infty(U_j)$ be the corresponding input cutoff used in the fixed local quantization formula. Define
\begin{align*}
R_{j,h}:C_c^\infty(U_j)\to C^\infty(U_j)
\end{align*}
by the localized kernel
\begin{align*}
K_{R_j}(x,y;h)=\chi_j(x)K_R(x,y;h)\psi_j(y).
\end{align*}
The previous two steps construct a local symbol $r_j$ on the $j$-th coordinate cotangent patch such that the coordinate-local quantization of $r_j$ has exactly the same localized kernel as $R_{j,h}$ on the region where the cutoffs are identically one.
On every compact subset of $T^*M$, only finitely many coordinate pieces contribute. The family $(r_j)_j$ is therefore an element of the symbol class associated with this fixed quantization atlas: in each chart it satisfies the rapidly vanishing $S^{-\infty}$ seminorm estimates, and the locally finite sum of the coordinate-local quantizations is well-defined. This is the sense in which the local symbols assemble into a global representative
\begin{align*}
r\in S^{-\infty}(T^*M).
\end{align*}
No coordinate transformation compatibility beyond the fixed quantization convention is being asserted; the conclusion is representation modulo residual equality, not equality of a coordinate-invariant principal symbol.
The difference between $R_h$ and the sum of the quantizations of the $r_j$ is supported only in the cutoff transition and cross-chart pieces left over from the fixed partition construction. Each such leftover has a smooth kernel obtained by multiplying the residual kernel of $R_h$ by fixed smooth cutoffs, differentiating only those fixed cutoffs and the residual kernel. Since the cutoff derivatives are independent of $h$ and since $R_h$ satisfies $O(h^N)$ bounds for every differentiated kernel seminorm, every leftover kernel is residual. Therefore every properly supported residual operator is represented modulo residual equality by a rapidly vanishing smoothing symbol. This completes the proof.
[/step]
