We study the Hermite operator H=−Δ+|x|2 in Rd and its fractional powers Hβ, β>0 in phase space. Namely, we represent functions f via the so-called short-time Fourier, alias Fourier-Wigner or Bargmann transform Vgf (g being a fixed window function), and we measure their regularity and decay by means of mixed Lebesgue norms in phase space of Vgf, that is in terms of membership to modulation spaces Mp,q, 0<p,q≤∞. We prove the complete range of fixed-time estimates for the semigroup e−tHβ when acting on Mp,q, for every 0<p,q≤∞, 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β 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−ct as the solution of the corresponding linear equation, where c=dβ is the bottom of the spectrum of Hβ. Global existence 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 (constant functions belong to M∞,1).