[proofplan]
The squared modulus of the coherent state is a centered Gaussian density translated to $x_0$, so its first moment is obtained by the change of variables $y=(x-x_0)/h^{1/2}$ and oddness of the centered Gaussian. We compute the semiclassical Fourier transform explicitly by completing the square; its modulus squared is the same Gaussian density translated to $\xi_0$. The second moments then follow from the one-dimensional Gaussian identity $\int_{\mathbb{R}} t^2 e^{-t^2}\,d\mathcal{L}^1(t)=\frac{1}{2}\int_{\mathbb{R}} e^{-t^2}\,d\mathcal{L}^1(t)$ applied coordinate by coordinate.
[/proofplan]
[step:Compute the position density and reduce its moments to centered Gaussian moments]
For $x \in \mathbb{R}^n$, the phase factor $\exp(i x\cdot \xi_0/h)$ has modulus $1$, hence
\begin{align*}
|u_h(x)|^2 = (\pi h)^{-n/2}\exp\left(-\frac{|x-x_0|^2}{h}\right).
\end{align*}
Define the affine change of variables
\begin{align*}
T_h: \mathbb{R}^n &\to \mathbb{R}^n, \qquad y \mapsto x_0 + h^{1/2}y.
\end{align*}
Under $x=T_h(y)$, Lebesgue measure transforms as $d\mathcal{L}^n(x)=h^{n/2}\,d\mathcal{L}^n(y)$ and $|x-x_0|^2=h|y|^2$. Therefore
\begin{align*}
\int_{\mathbb{R}^n} x |u_h(x)|^2\,d\mathcal{L}^n(x) = \pi^{-n/2}\int_{\mathbb{R}^n} (x_0+h^{1/2}y)e^{-|y|^2}\,d\mathcal{L}^n(y).
\end{align*}
The function $y \mapsto y e^{-|y|^2}$ is odd componentwise and integrable on $\mathbb{R}^n$, so its integral over the symmetric domain $\mathbb{R}^n$ is $0$. Also
\begin{align*}
\pi^{-n/2}\int_{\mathbb{R}^n} e^{-|y|^2}\,d\mathcal{L}^n(y)=1.
\end{align*}
Thus
\begin{align*}
\int_{\mathbb{R}^n} x |u_h(x)|^2\,d\mathcal{L}^n(x)=x_0.
\end{align*}
[guided]
The spatial density is obtained by taking the squared modulus of $u_h$. The oscillatory factor $\exp(i x\cdot \xi_0/h)$ only records frequency information and has absolute value $1$, so it disappears from the probability density:
\begin{align*}
|u_h(x)|^2 = (\pi h)^{-n/2}\exp\left(-\frac{|x-x_0|^2}{h}\right).
\end{align*}
To see the center, we move the Gaussian back to the origin and scale it to unit width. Define
\begin{align*}
T_h: \mathbb{R}^n &\to \mathbb{R}^n, \qquad y \mapsto x_0+h^{1/2}y.
\end{align*}
With $x=T_h(y)$, the Jacobian determinant is $h^{n/2}$, so
\begin{align*}
d\mathcal{L}^n(x)=h^{n/2}\,d\mathcal{L}^n(y).
\end{align*}
Also $|x-x_0|^2=h|y|^2$. Substituting gives
\begin{align*}
\int_{\mathbb{R}^n} x |u_h(x)|^2\,d\mathcal{L}^n(x) = \pi^{-n/2}\int_{\mathbb{R}^n} (x_0+h^{1/2}y)e^{-|y|^2}\,d\mathcal{L}^n(y).
\end{align*}
Now separate the constant part from the odd part. The Gaussian integral satisfies
\begin{align*}
\pi^{-n/2}\int_{\mathbb{R}^n} e^{-|y|^2}\,d\mathcal{L}^n(y)=1,
\end{align*}
and each coordinate of $y e^{-|y|^2}$ is odd in one variable and integrable over the symmetric domain $\mathbb{R}^n$. Hence
\begin{align*}
\int_{\mathbb{R}^n} y e^{-|y|^2}\,d\mathcal{L}^n(y)=0.
\end{align*}
Therefore the only surviving contribution is the translated center:
\begin{align*}
\int_{\mathbb{R}^n} x |u_h(x)|^2\,d\mathcal{L}^n(x)=x_0.
\end{align*}
[/guided]
[/step]
[step:Compute the semiclassical Fourier transform by completing the square]
For $\xi \in \mathbb{R}^n$,
\begin{align*}
\mathcal{F}_h u_h(\xi) = (2\pi h)^{-n/2}(\pi h)^{-n/4}\int_{\mathbb{R}^n}\exp\left(-\frac{|x-x_0|^2}{2h}\right)\exp\left(\frac{i x\cdot(\xi_0-\xi)}{h}\right)\,d\mathcal{L}^n(x).
\end{align*}
Use the change of variables $z=x-x_0$, so $d\mathcal{L}^n(x)=d\mathcal{L}^n(z)$ and the domain remains $\mathbb{R}^n$. Then
\begin{align*}
\mathcal{F}_h u_h(\xi) = (2\pi h)^{-n/2}(\pi h)^{-n/4}\exp\left(\frac{i x_0\cdot(\xi_0-\xi)}{h}\right)\int_{\mathbb{R}^n}\exp\left(-\frac{|z|^2}{2h}\right)\exp\left(\frac{i z\cdot(\xi_0-\xi)}{h}\right)\,d\mathcal{L}^n(z).
\end{align*}
The standard Gaussian completion of the square gives
\begin{align*}
\int_{\mathbb{R}^n}\exp\left(-\frac{|z|^2}{2h}\right)\exp\left(\frac{i z\cdot(\xi_0-\xi)}{h}\right)\,d\mathcal{L}^n(z) = (2\pi h)^{n/2}\exp\left(-\frac{|\xi-\xi_0|^2}{2h}\right).
\end{align*}
Consequently
\begin{align*}
\mathcal{F}_h u_h(\xi) = (\pi h)^{-n/4}\exp\left(-\frac{|\xi-\xi_0|^2}{2h}\right)\exp\left(-\frac{i x_0\cdot(\xi-\xi_0)}{h}\right).
\end{align*}
Taking absolute values removes the phase factor and yields
\begin{align*}
|\mathcal{F}_h u_h(\xi)|^2 = (\pi h)^{-n/2}\exp\left(-\frac{|\xi-\xi_0|^2}{h}\right).
\end{align*}
[/step]
[step:Use the transformed density to identify the frequency center]
Define the affine change of variables
\begin{align*}
S_h: \mathbb{R}^n &\to \mathbb{R}^n, \qquad \eta \mapsto \xi_0+h^{1/2}\eta.
\end{align*}
Under $\xi=S_h(\eta)$, Lebesgue measure transforms as $d\mathcal{L}^n(\xi)=h^{n/2}\,d\mathcal{L}^n(\eta)$ and $|\xi-\xi_0|^2=h|\eta|^2$. Using the density computed above,
\begin{align*}
\int_{\mathbb{R}^n} \xi |\mathcal{F}_h u_h(\xi)|^2\,d\mathcal{L}^n(\xi) = \pi^{-n/2}\int_{\mathbb{R}^n}(\xi_0+h^{1/2}\eta)e^{-|\eta|^2}\,d\mathcal{L}^n(\eta).
\end{align*}
The same normalization and oddness argument as in the position computation gives
\begin{align*}
\int_{\mathbb{R}^n} \xi |\mathcal{F}_h u_h(\xi)|^2\,d\mathcal{L}^n(\xi)=\xi_0.
\end{align*}
[/step]
[step:Compute the position and frequency variances from the second centered Gaussian moment]
For the position variance, use $x=x_0+h^{1/2}y$ as above:
\begin{align*}
\int_{\mathbb{R}^n}|x-x_0|^2 |u_h(x)|^2\,d\mathcal{L}^n(x) = h\pi^{-n/2}\int_{\mathbb{R}^n}|y|^2 e^{-|y|^2}\,d\mathcal{L}^n(y).
\end{align*}
Since $|y|^2=\sum_{j=1}^n y_j^2$ and the Gaussian factorizes as $e^{-|y|^2}=\prod_{k=1}^n e^{-y_k^2}$, each coordinate contributes
\begin{align*}
\pi^{-n/2}\int_{\mathbb{R}^n} y_j^2 e^{-|y|^2}\,d\mathcal{L}^n(y)=\frac{1}{2}.
\end{align*}
Indeed, this follows from the one-dimensional identity
\begin{align*}
\int_{\mathbb{R}} t^2 e^{-t^2}\,d\mathcal{L}^1(t)=\frac{1}{2}\int_{\mathbb{R}} e^{-t^2}\,d\mathcal{L}^1(t),
\end{align*}
obtained by integration by parts using the differentiable function $t \mapsto e^{-t^2}$. Summing over $j=1,\dots,n$ gives
\begin{align*}
\int_{\mathbb{R}^n}|x-x_0|^2 |u_h(x)|^2\,d\mathcal{L}^n(x)=\frac{nh}{2}.
\end{align*}
The frequency density has the same form with $x_0$ replaced by $\xi_0$. Applying the same calculation with $\xi=\xi_0+h^{1/2}\eta$ gives
\begin{align*}
\int_{\mathbb{R}^n}|\xi-\xi_0|^2 |\mathcal{F}_h u_h(\xi)|^2\,d\mathcal{L}^n(\xi)=\frac{nh}{2}.
\end{align*}
Thus the total position and frequency variances are both proportional to $h$, and each coordinate variance is $h/2$.
[/step]
