  Let $X_{1},...,X_{m}$ be a family of real smooth vector fields defined in
$\mathbb{R}^{n}$, $1$-homogeneous with respect to a nonisotropic family of
dilations and satisfying H\"{o}rmander's rank condition at $0$ (and therefore
at every point of $\mathbb{R}^{n}$). The vector fields are not assumed to be
translation invariant with respect to any Lie group structure. Let us consider
the nonvariational evolution operator $$
\mathcal{H}:=\sum_{i,j=1}^{m}a_{i,j}(t,x)X_{i}X_{j}-\partial_{t}% $$ where
$(a_{i,j}(t,x))_{i,j=1}^{m}$ is a symmetric uniformly positive $m\times m$
matrix and the entries $a_{ij}$ are bounded H\"{o}lder continuous functions on
$\mathbb{R}^{1+n}$, with respect to the "parabolic" distance induced by the
vector fields. We prove the existence of a global heat kernel
$\Gamma(\cdot;s,y)\in
C_{X,\mathrm{loc}}^{2,\alpha}(\mathbb{R}^{1+n}\setminus\{(s,y)\})$ for
$\mathcal{H}$, such that $\Gamma$ satisfies two-sided Gaussian bounds and
$\partial_{t}\Gamma, X_{i}\Gamma,X_{i}X_{j}\Gamma$ satisfy upper Gaussian
bounds on every strip $[0,T]\times\mathbb{R}^n$. We also prove a
scale-invariant parabolic Harnack inequality for $\mathcal{H}$, and a standard
Harnack inequality for the corresponding stationary operator $$
  \mathcal{L}:=\sum_{i,j=1}^{m}a_{i,j}(x)X_{i}X_{j}.
  $$
  with H\"{o}lder continuos coefficients.
