Non-divergence operators structured on homogeneous H\"{o}rmander vector fields: heat kernels and global Gaussian bounds

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.