On Li-Yau gradient estimate for sum of squares of vector fields up to higher step

In this paper, we generalize the Cao-Yau's gradient estimate for the sum of squares of vector fields up to higher step under assumption of the generalized curvature-dimension inequality. With its applications, by deriving a curvature-dimension inequality, we are able to obtain the Li-Yau gradient estimate for the CR heat equation in a closed pseudohermitian manifold of nonvanishing torsion tensors. As consequences, we obtain the Harnack inequality and upper bound estimate for the CR heat kernel.