[proofplan]
We first regularize the Weyl oscillatory integral by inserting a symmetric compact cutoff in the variables $(x,y,\xi)$, so every expression becomes an absolutely convergent Lebesgue integral. For the regularized operators, the adjoint identity follows from taking the complex conjugate in the $L^2$ inner product and interchanging the variables $x$ and $y$; the Weyl midpoint $(x+y)/2$ is exactly what makes the symbol transform into $\overline{a}((x+y)/2,\xi)$. Finally we pass to the oscillatory integral limit, using the defining continuity of Weyl quantization on $\mathcal{S}(\mathbb{R}^n)$ and $\mathcal{S}'(\mathbb{R}^n)$.
[/proofplan]
[step:Regularize the Weyl pairing by a symmetric cutoff]
Choose a function
\begin{align*}
\chi:\mathbb{R}^n \times \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R}
\end{align*}
with $\chi \in C_c^\infty(\mathbb{R}^{3n})$, $\chi(x,y,\xi)=\chi(y,x,\xi)$ for all $(x,y,\xi)$, and $\chi=1$ on a neighbourhood of $(0,0,0)$. For $\varepsilon \in (0,1]$, define
\begin{align*}
\chi_\varepsilon:\mathbb{R}^n \times \mathbb{R}^n \times \mathbb{R}^n &\to \mathbb{R}
\end{align*}
\begin{align*}
(x,y,\xi) &\mapsto \chi(\varepsilon x,\varepsilon y,\varepsilon \xi).
\end{align*}
For each $\varepsilon \in (0,1]$, define the regularized pairing as the map
\begin{align*}
I_\varepsilon(a;\cdot,\cdot):\mathcal{S}(\mathbb{R}^n) \times \mathcal{S}(\mathbb{R}^n) \to \mathbb{C}.
\end{align*}
For $u,v \in \mathcal{S}(\mathbb{R}^n)$, its value is
\begin{align*}
I_\varepsilon(a;u,v)
=
(2\pi h)^{-n}
\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}
\chi_\varepsilon(x,y,\xi)
e^{i(x-y)\cdot \xi/h}
a\!\left(\frac{x+y}{2},\xi\right)
u(y)\overline{v(x)}
\, d\mathcal{L}^n(\xi)\, d\mathcal{L}^n(y)\, d\mathcal{L}^n(x).
\end{align*}
Because $\chi_\varepsilon$ has compact support and the remaining factors are smooth on that compact set, this is an absolutely convergent Lebesgue integral. We use the standard cutoff-independence theorem for oscillatory integrals: if an amplitude is a symbol in the oscillatory variables and is paired with Schwartz functions in the remaining variables, then inserting any cutoff $\chi_\varepsilon(x,y,\xi)=\chi(\varepsilon x,\varepsilon y,\varepsilon \xi)$ with $\chi=1$ near the origin gives the same oscillatory integral in the limit $\varepsilon \to 0^+$. Its hypotheses apply here because $a \in S^m(T^*\mathbb{R}^n)$ has polynomial growth with symbol estimates, while $u,v \in \mathcal{S}(\mathbb{R}^n)$ provide rapid decay in the spatial variables; after the usual integration-by-parts operator in $\xi$ and in $x-y$, the regularized integrands are dominated by integrable seminorm bounds independent of $\varepsilon$ on compact parameter ranges. Therefore the cutoff limit represents the distributional pairing of the Weyl operator with $v$:
\begin{align*}
(\operatorname{Op}_h^w(a)u,v)_{L^2(\mathbb{R}^n)}
=
\lim_{\varepsilon \to 0^+} I_\varepsilon(a;u,v).
\end{align*}
[/step]
[step:Expand the regularized right-hand pairing and exchange the two spatial variables]
For each $\varepsilon \in (0,1]$, define the regularized right-hand pairing as the map
\begin{align*}
J_\varepsilon(\overline a;\cdot,\cdot):\mathcal{S}(\mathbb{R}^n) \times \mathcal{S}(\mathbb{R}^n) \to \mathbb{C}.
\end{align*}
For $u,v \in \mathcal{S}(\mathbb{R}^n)$, its value at $(v,u)$ is
\begin{align*}
J_\varepsilon(\overline a;v,u)
=
(2\pi h)^{-n}
\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}
\chi_\varepsilon(x,y,\xi)
u(x)
\overline{
e^{i(x-y)\cdot \xi/h}
\overline a\!\left(\frac{x+y}{2},\xi\right)
v(y)
}
\, d\mathcal{L}^n(\xi)\, d\mathcal{L}^n(y)\, d\mathcal{L}^n(x).
\end{align*}
This integral is absolutely convergent because $\chi_\varepsilon$ has compact support. Since $\chi_\varepsilon$ is real-valued, taking the complex conjugate inside the integrand gives
\begin{align*}
J_\varepsilon(\overline a;v,u)
=
(2\pi h)^{-n}
\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}
\chi_\varepsilon(x,y,\xi)
e^{i(y-x)\cdot \xi/h}
a\!\left(\frac{x+y}{2},\xi\right)
u(x)\overline{v(y)}
\, d\mathcal{L}^n(\xi)\, d\mathcal{L}^n(y)\, d\mathcal{L}^n(x).
\end{align*}
Now apply the linear change of variables $(x,y,\xi) \mapsto (y,x,\xi)$. Its absolute Jacobian determinant is $1$, so the product measure $\mathcal{L}^n \otimes \mathcal{L}^n \otimes \mathcal{L}^n$ is preserved. The cutoff is unchanged because $\chi_\varepsilon(y,x,\xi)=\chi_\varepsilon(x,y,\xi)$, and the Weyl midpoint is unchanged because
\begin{align*}
\frac{y+x}{2}=\frac{x+y}{2}.
\end{align*}
Therefore
\begin{align*}
J_\varepsilon(\overline a;v,u)
=
(2\pi h)^{-n}
\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}
\chi_\varepsilon(x,y,\xi)
e^{i(x-y)\cdot \xi/h}
a\!\left(\frac{x+y}{2},\xi\right)
u(y)\overline{v(x)}
\, d\mathcal{L}^n(\xi)\, d\mathcal{L}^n(y)\, d\mathcal{L}^n(x).
\end{align*}
The right-hand side is exactly $I_\varepsilon(a;u,v)$.
[guided]
The desired identity compares two different pairings, not $I_\varepsilon(a;u,v)$ with its own conjugate. Because the $L^2(\mathbb{R}^n)$ inner product is linear in the first entry and conjugate-linear in the second, the regularized version of $(u,\operatorname{Op}_h^w(\overline a)v)_{L^2(\mathbb{R}^n)}$ is
\begin{align*}
J_\varepsilon(\overline a;v,u)
=
(2\pi h)^{-n}
\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}
\chi_\varepsilon(x,y,\xi)
u(x)
\overline{
e^{i(x-y)\cdot \xi/h}
\overline a\!\left(\frac{x+y}{2},\xi\right)
v(y)
}
\, d\mathcal{L}^n(\xi)\, d\mathcal{L}^n(y)\, d\mathcal{L}^n(x).
\end{align*}
Now we take the conjugate of the operator output. The phase changes sign, the conjugated symbol $\overline a$ becomes $a$, and $v(y)$ becomes $\overline{v(y)}$:
\begin{align*}
\overline{
e^{i(x-y)\cdot \xi/h}
\overline a\!\left(\frac{x+y}{2},\xi\right)
v(y)
}
=
e^{i(y-x)\cdot \xi/h}
a\!\left(\frac{x+y}{2},\xi\right)
\overline{v(y)}.
\end{align*}
Thus
\begin{align*}
J_\varepsilon(\overline a;v,u)
=
(2\pi h)^{-n}
\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}
\chi_\varepsilon(x,y,\xi)
e^{i(y-x)\cdot \xi/h}
a\!\left(\frac{x+y}{2},\xi\right)
u(x)\overline{v(y)}
\, d\mathcal{L}^n(\xi)\, d\mathcal{L}^n(y)\, d\mathcal{L}^n(x).
\end{align*}
We now exchange the spatial variables. The map $T:\mathbb{R}^{3n}\to\mathbb{R}^{3n}$ defined by $T(x,y,\xi)=(y,x,\xi)$ is linear and has absolute Jacobian determinant $1$, so the measure $d\mathcal{L}^n(\xi)\, d\mathcal{L}^n(y)\, d\mathcal{L}^n(x)$ is preserved. The cutoff is symmetric, so $\chi_\varepsilon(y,x,\xi)=\chi_\varepsilon(x,y,\xi)$. The midpoint also survives the exchange:
\begin{align*}
\frac{y+x}{2}=\frac{x+y}{2}.
\end{align*}
After this change of variables, the phase $e^{i(y-x)\cdot \xi/h}$ becomes $e^{i(x-y)\cdot \xi/h}$, the factor $u(x)$ becomes $u(y)$, and $\overline{v(y)}$ becomes $\overline{v(x)}$. Hence
\begin{align*}
J_\varepsilon(\overline a;v,u)
=
(2\pi h)^{-n}
\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}\int_{\mathbb{R}^n}
\chi_\varepsilon(x,y,\xi)
e^{i(x-y)\cdot \xi/h}
a\!\left(\frac{x+y}{2},\xi\right)
u(y)\overline{v(x)}
\, d\mathcal{L}^n(\xi)\, d\mathcal{L}^n(y)\, d\mathcal{L}^n(x).
\end{align*}
This is exactly $I_\varepsilon(a;u,v)$. The midpoint in the Weyl quantization is the structural reason the argument works: exchanging $x$ and $y$ leaves $(x+y)/2$ fixed, so the same quantization convention reappears after the adjoint operation.
[/guided]
[/step]
[step:Pass from the regularized identity to the oscillatory integral identity]
From the preceding step, for every $\varepsilon \in (0,1]$,
\begin{align*}
I_\varepsilon(a;u,v)=J_\varepsilon(\overline a;v,u).
\end{align*}
Taking the limit as $\varepsilon \to 0^+$ on both sides gives
\begin{align*}
(\operatorname{Op}_h^w(a)u,v)_{L^2(\mathbb{R}^n)}
=
(u,\operatorname{Op}_h^w(\overline a)v)_{L^2(\mathbb{R}^n)}.
\end{align*}
The left-hand side is the cutoff-independent oscillatory limit for the pairing of $\operatorname{Op}_h^w(a)u$ with $v$. The right-hand side is the corresponding cutoff-independent oscillatory limit for the pairing of $u$ with $\operatorname{Op}_h^w(\overline a)v$; the same oscillatory-integral theorem applies because $\overline a \in S^m(T^*\mathbb{R}^n)$ whenever $a \in S^m(T^*\mathbb{R}^n)$, and $u,v \in \mathcal{S}(\mathbb{R}^n)$.
[/step]
[step:Identify the formal adjoint and the self-adjoint real-symbol case]
The formal adjoint $(\operatorname{Op}_h^w(a))^*$ is defined by the identity
\begin{align*}
((\operatorname{Op}_h^w(a))u,v)_{L^2(\mathbb{R}^n)}
=
(u,(\operatorname{Op}_h^w(a))^*v)_{L^2(\mathbb{R}^n)}
\end{align*}
for all $u,v \in \mathcal{S}(\mathbb{R}^n)$, interpreted as an identity of distributions. The identity proved above therefore gives
\begin{align*}
(\operatorname{Op}_h^w(a))^*=\operatorname{Op}_h^w(\overline a)
\end{align*}
as continuous maps $\mathcal{S}(\mathbb{R}^n)\to \mathcal{S}'(\mathbb{R}^n)$ after taking formal adjoints. If $a$ is real-valued, then $\overline a=a$, and hence
\begin{align*}
(\operatorname{Op}_h^w(a))^*=\operatorname{Op}_h^w(a).
\end{align*}
Thus $\operatorname{Op}_h^w(a)$ is formally self-adjoint.
[/step]
