In this work we deal with linear second order partial differential operators of the following type: H = partial derivative(t) - L = partial derivative(t) - Sigma(q)(i,j=1) a(ij) (t, x) X(i)X(j) - Sigma(q)(k=1) a(k) (t, X) X(k) - a(0) (t, x) where X(1), X(2), ..., X(q) is a system of real Hormander's vector fields in some bounded domain Omega subset of R(n), A = {a(ij) (t, x)}(i,j=1)(q) is a real symmetric uniformly positive definite matrix such that: lambda(-1) vertical bar xi vertical bar(2) <= Sigma(q)(i,j=1) a(ij) (t, x) xi(i)xi(j) <= lambda vertical bar xi vertical bar(2) for all xi is an element of R(q), x is an element of Omega, t is an element of (T(1), T(2)) for a suitable constant lambda > 0 a for some real numbers T(1) < T(2). The coefficients a(ij,) a(k), a(0) are Hblder continuous on (T(1), T(2)) x Omega with respect to the parabolic CG-metric d(p) ((t, x), (s, y)) = root d (x, y)(2) + vertical bar t - s vertical bar (where d is the Carnot-Caratheodory distance induced by the vector fields X(i)'s). We prove the existence of a fundamental solution h (t, x; s, y) for H, satisfying natural properties and sharp Gaussian bounds of the kind: [GRAPHICS] where vertical bar B (x, r)vertical bar denotes the Lebesgue measure of the d-ball B (x, r). We then use these properties of h as a starting point to prove a scaling invariant Harnack inequality for positive solutions to Hu = 0, when a(0) equivalent to 0. All the constants in our estimates and inequalities will depend on the coefficients a(ij), a(k), a(0) only through their Holder norms and the number lambda.