  We study the Hermite operator $H=-\Delta+|x|^2$ in $\mathbb{R}^d$ and its
fractional powers $H^\beta$, $\beta>0$ in phase space. Namely, we represent
functions $f$ via the so-called short-time Fourier, alias Fourier-Wigner or
Bargmann transform $V_g f$ ($g$ being a fixed window function), and we measure
their regularity and decay by means of mixed Lebesgue norms in phase space of
$V_g f$, that is in terms of membership to modulation spaces $M^{p,q}$, $0<
p,q\leq \infty$. We prove the complete range of fixed-time estimates for the
semigroup $e^{-tH^\beta}$ when acting on $M^{p,q}$, for every $0< p,q\leq
\infty$, exhibiting the optimal global-in-time decay as well as phase-space
smoothing. As an application, we establish global well-posedness for the
nonlinear heat equation for $H^{\beta}$ with power-type nonlinearity (focusing
or defocusing), with small initial data in modulation spaces or in Wiener
amalgam spaces. We show that such a global solution exhibits the same optimal
decay $e^{-c t}$ as the solution of the corresponding linear equation, where
$c=d^\beta$ is the bottom of the spectrum of $H^\beta$. This is in sharp
contrast to what happens for the nonlinear focusing heat equation without
potential, where blow-up in finite time always occurs for (even small) constant
initial data - hence in $M^{\infty,1}$.
