We study the Hermite operator H = −∆ + |x| 2 in R d and its fractional powers H β , β > 0 in phase space. Namely, we represent functions f via the socalled 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 ≤ ∞. We prove the complete range of fixed-time estimates for the semigroup e −tH β when acting on M p,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 β . 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 ∞,1 .