[proofplan]
We take the pointwise positive square root $b_0=a^{1/2}$ and prove that ellipticity of $a$ relative to $m$ forces $b_0$ to lie in the symbol class $S(m^{1/2})$. The Weyl product of a real symbol with itself has no first-order $h$ correction because the Poisson bracket of a symbol with itself vanishes. The standard Weyl symbolic expansion therefore gives $b_0 \# b_0=a+h^2r_h$ with $r_h$ uniformly bounded in $S(m)$. Since $b_0$ is real-valued, Weyl quantization is formally self-adjoint, and the operator identity follows from the Weyl composition formula.
[/proofplan]
[step:Take the positive square root of the elliptic symbol]
Let $z=(x,\xi)$ denote a point of $T^*\mathbb{R}^n \cong \mathbb{R}^{2n}$. Define
\begin{align*}
b_0: T^*\mathbb{R}^n &\to \mathbb{R}
\end{align*}
by
\begin{align*}
b_0(z)=a(z)^{1/2}.
\end{align*}
This is well-defined and real-valued because $a(z)\geq c\,m(z)>0$ for every $z\in T^*\mathbb{R}^n$. For $0<h\leq h_0$, set
\begin{align*}
b_h:T^*\mathbb{R}^n &\to \mathbb{R}
\end{align*}
by $b_h(z)=b_0(z)$. Thus the symbol family is independent of $h$.
[/step]
[step:Prove that the square root has order $m^{1/2}$]
We prove $b_0\in S(m^{1/2})$. Let $\alpha\in\mathbb{N}_0^{2n}$ be a multi-index. Since $a\in S(m)$, for every multi-index $\beta\in\mathbb{N}_0^{2n}$ there is a constant $C_\beta>0$ such that
\begin{align*}
|\partial_z^\beta a(z)|\leq C_\beta m(z)
\end{align*}
for every $z\in T^*\mathbb{R}^n$.
For $\alpha=0$, the symbol bound follows from the order-function bound for $a$:
\begin{align*}
|b_0(z)|=a(z)^{1/2}\leq C_0^{1/2}m(z)^{1/2}.
\end{align*}
For $|\alpha|\geq 1$, the multivariable chain rule and Faà di Bruno formula express $\partial_z^\alpha b_0$ as a finite sum indexed by a finite set $\mathcal{J}_\alpha$. Each $J\in\mathcal{J}_\alpha$ determines an integer $k=k(J)$, a real combinatorial constant $C_{\alpha,J}$, and nonzero multi-indices $\beta_1,\dots,\beta_k\in\mathbb{N}_0^{2n}$, and the corresponding term has the form
\begin{align*}
C_{\alpha,J}\,a(z)^{1/2-k}\prod_{j=1}^{k}\partial_z^{\beta_j}a(z),
\end{align*}
where $1\leq k\leq |\alpha|$ and the multi-indices satisfy
\begin{align*}
\beta_1+\cdots+\beta_k=\alpha.
\end{align*}
Using $a(z)\geq c\,m(z)$ and $|\partial_z^{\beta_j}a(z)|\leq C_{\beta_j}m(z)$, each such term is bounded by a constant depending only on $\alpha$, $c$, and finitely many symbol seminorms of $a$, times
\begin{align*}
m(z)^{1/2-k}m(z)^k=m(z)^{1/2}.
\end{align*}
Hence for every $\alpha\in\mathbb{N}_0^{2n}$ there is a constant $C_\alpha'>0$ such that
\begin{align*}
|\partial_z^\alpha b_0(z)|\leq C_\alpha' m(z)^{1/2}
\end{align*}
for every $z\in T^*\mathbb{R}^n$. Therefore $b_0\in S(m^{1/2})$, and the family $b_h=b_0$ is uniformly bounded in $S(m^{1/2})$ for $0<h\leq h_0$.
[guided]
The only possible obstruction to taking the square root is differentiating $a^{1/2}$: derivatives of the square-root function produce negative powers of $a$. The ellipticity lower bound prevents these negative powers from becoming singular, because $a(z)\geq c\,m(z)$ and $m(z)>0$ everywhere.
Let $\alpha\in\mathbb{N}_0^{2n}$ be a multi-index. Since $a\in S(m)$, every derivative of $a$ is controlled by $m$: for each $\beta\in\mathbb{N}_0^{2n}$ there exists $C_\beta>0$ such that
\begin{align*}
|\partial_z^\beta a(z)|\leq C_\beta m(z)
\end{align*}
for all $z\in T^*\mathbb{R}^n$. For $\alpha=0$, this gives
\begin{align*}
|b_0(z)|=a(z)^{1/2}\leq C_0^{1/2}m(z)^{1/2}.
\end{align*}
Now suppose $|\alpha|\geq 1$. The Faà di Bruno formula for the composition $t\mapsto t^{1/2}$ with $a$ says that $\partial_z^\alpha(a^{1/2})$ is a finite sum of products in which $k$ derivatives of $a$ are multiplied by the $k$-th derivative of $t^{1/2}$ evaluated at $t=a(z)$. Thus every term has the form
\begin{align*}
C_{\alpha,J}\,a(z)^{1/2-k}\prod_{j=1}^{k}\partial_z^{\beta_j}a(z),
\end{align*}
where $J$ ranges over a finite set $\mathcal{J}_\alpha$, $1\leq k\leq |\alpha|$, the multi-indices $\beta_j$ are nonzero, and $\beta_1+\cdots+\beta_k=\alpha$.
We estimate this term factor by factor. The lower bound $a(z)\geq c\,m(z)$ gives
\begin{align*}
a(z)^{1/2-k}\leq c^{1/2-k}m(z)^{1/2-k}
\end{align*}
because $1/2-k<0$ when $k\geq 1$. The symbol estimates for $a$ give
\begin{align*}
\prod_{j=1}^{k}|\partial_z^{\beta_j}a(z)|
\leq
\left(\prod_{j=1}^{k}C_{\beta_j}\right)m(z)^k.
\end{align*}
Multiplying the two estimates yields
\begin{align*}
\left|C_{\alpha,J}\,a(z)^{1/2-k}\prod_{j=1}^{k}\partial_z^{\beta_j}a(z)\right|
\leq
C_{\alpha,J}'m(z)^{1/2}.
\end{align*}
There are only finitely many terms for the fixed multi-index $\alpha$, so their sum satisfies the same type of bound:
\begin{align*}
|\partial_z^\alpha b_0(z)|\leq C_\alpha' m(z)^{1/2}.
\end{align*}
This is precisely the defining estimate for $b_0\in S(m^{1/2})$.
[/guided]
[/step]
[step:Use the Weyl product expansion to remove the first-order term]
We use the standard semiclassical Weyl symbolic calculus for order functions: if $m_1:T^*\mathbb{R}^n\to(0,\infty)$ and $m_2:T^*\mathbb{R}^n\to(0,\infty)$ are order functions, if $p\in S(m_1)$ and $q\in S(m_2)$, and if $0<h\leq h_0$, then
\begin{align*}
p\# q = pq+\frac{ih}{2}\{p,q\}+h^2R_h(p,q)
\end{align*}
where $\{\cdot,\cdot\}$ is the canonical Poisson bracket on $T^*\mathbb{R}^n$ and where $R_h(p,q):T^*\mathbb{R}^n\to\mathbb{C}$ is a symbol family bounded uniformly in $S(m_1m_2)$ for $0<h\leq h_0$, with bounds controlled by finitely many seminorms of $p$ and $q$.
Apply this with the order functions $m_1=m^{1/2}$ and $m_2=m^{1/2}$ and with the symbols $p=q=b_0$. Since $b_0\in S(m^{1/2})$, the hypotheses of the expansion are satisfied and the remainder belongs uniformly to $S(m)$. Also
\begin{align*}
\{b_0,b_0\}=0
\end{align*}
by antisymmetry of the Poisson bracket. Therefore there exists a family $r_h$ bounded uniformly in $S(m)$ such that
\begin{align*}
b_0\# b_0=b_0^2+h^2r_h.
\end{align*}
Since $b_0^2=a$ pointwise, this becomes
\begin{align*}
b_h\# b_h-a=h^2r_h.
\end{align*}
Thus $b_h\# b_h-a\in h^2S(m)$.
[/step]
[step:Pass from the symbolic identity to the operator identity]
We interpret $*$ as the formal adjoint on Schwartz functions, as stated in the theorem. Because $b_h$ is real-valued, the Weyl quantization $\operatorname{Op}_h^w(b_h)$ is formally self-adjoint on Schwartz functions, hence
\begin{align*}
(\operatorname{Op}_h^w(b_h))^*=\operatorname{Op}_h^w(b_h).
\end{align*}
The Weyl composition formula for symbols in $S(m^{1/2})$ gives
\begin{align*}
\operatorname{Op}_h^w(b_h)\operatorname{Op}_h^w(b_h)=\operatorname{Op}_h^w(b_h\# b_h).
\end{align*}
Using the symbolic identity $b_h\# b_h=a+h^2r_h$, we obtain
\begin{align*}
(\operatorname{Op}_h^w(b_h))^*\operatorname{Op}_h^w(b_h)
=
\operatorname{Op}_h^w(a+h^2r_h).
\end{align*}
By linearity of Weyl quantization,
\begin{align*}
\operatorname{Op}_h^w(a+h^2r_h)=\operatorname{Op}_h^w(a)+h^2\operatorname{Op}_h^w(r_h).
\end{align*}
The family $r_h$ is bounded uniformly in $S(m)$ for $0<h\leq h_0$, so this is exactly the asserted operator identity on Schwartz functions.
[/step]
