Uniform Gaussian bounds for subelliptic heat kernels and an application to the total variation flow of graphs over Carnot groups

In this paper we study heat kernels associated to a Carnot group $G$, endowed with a family of collapsing left-invariant Riemannian metrics $\sigma_\e$ which converge in the Gromov-Hausdorff sense to a sub-Riemannian structure on $G$ as $\e\to 0$. The main new contribution are Gaussian-type bounds on the heat kernel for the $\sigma_\e$ metrics which are stable as $\e\to 0$ and extend the previous time-independent estimates in \cite{CiMa-F}. As an application we study well posedness of the total variation flow of graph surfaces over a bounded domain in $(G,\s_\e)$. We establish interior and boundary gradient estimates, and develop a Schauder theory which are stable as $\e\to 0$. As a consequence we obtain long time existence of smooth solutions of the sub-Riemannian flow ($\e=0$), which in turn yield sub-Riemannian minimal surfaces as $t\to \infty$.


Introduction
Models of image processing based on total variation flow ∂ t u = div(∇u/|∇u|) have been first introduced in by Rudin, Osher, and Fatemi in [43], in order to perform edge preserving denoising. We refer the reader to the review paper [19] where many recent applications of total variation equation in the Euclidean setting are presented. The flow t → u(·, t) represents the gradient descent associated to the total variation energy |∇u| and as such has the property that both the total variation of the solution t → u(·, t) and the perimeter measure of fixed level sets {u(·, t) = const} are non-increasing in time. Aside from its usefulness in image processing, the flow also arises in connection with the limit of solutions of the parabolic p−Laplacian ∂ t u p = div(|∇u p | p−2 ∇u p ) as the parameter p → 1 + . In the case where the evolution of graphs S t = {(x, u(x, t))} is considered, i.e. ∂ t u = div(∇u/ 1 + |∇u| 2 ), then in both the total variation flow and the closely related mean curvature flow ∂ t u = 1 + |∇u| 2 div(∇u/ 1 + |∇u| 2 ), given appropriate boundary/initial conditions, global in time solutions asymptotically converge to minimal graphs. For further results, references and applications of the Euclidean total variation flow we refer the reader to [1], [2], [5] and for an anisotropic version of the flow to [36].
In this paper we study long time existence of graph solutions of the total variation flow in a special class of degenerate Riemannian ambient spaces: The so-called (sub-Riemannian) Carnot groups [26], [46]. Such spaces are nilpotent Lie groups endowed with a metric structure (G, σ 0 ) that arises as limit of collapsing left-invariant tame Riemannian structures (G, σ ǫ ). Our motivations are twofold: (a) On the one hand, minimal surfaces and mean curvature flow in the Carnot group setting (and in particular in the special case of the Heisenberg group) have been the subject of numerous papers in recent years, leading to partial solutions of long-standing problems such as the Pansu conjecture 1 . Despite such advances little is known about explicit construction of minimal surfaces. Since the asymptotic limit of graphical mean curvature flow provide minimal surfaces in the Riemannian setting it is only natural to follow the same approach in the sub-Riemannian setting. However this avenue runs into considerable (so far unsolved) technical difficulties due to the combined effect of the degeneracy of the metric and the nondivergence form aspect of the relevant PDE. In the present work we provide the basis for the construction of sub-Riemannian minimal graphs through the asymptotic behavior as t → ∞ of solutions of the sub-Riemannian total variation flow. For previous work on a similar theme see the work of Pauls [40] and Cheng, Hwang and Yang [14].
(b) A new class of image processing models based on the functionality of the visual cortex have been recently introduced in Lie groups. With a generalization of the classical Bargmann transform studied by Folland (see [27]) which lifts a L 2 function to a new one, defined on the phase space, an image can be lifted to a Lie group with a sub-Riemannian metric. The choice of the Lie group depends on the geometric property of the image to be studied: in [18], and [23], [24] 2D images are lifted to surfaces in the Lie group of Euclidean motions of the plane to study geometric properties of their level lines, in [25] cardiac images are lifted in the Heisenberg group to study their deformations, in [4] moving images are lifted in the Galilei group. The image processing is performed in these groups with algorithms expressed in terms of second order subelliptic differential equations. In particular sub-Riemannian mean curvature flows or total variation flows can be applied to perform image completion or impainting in this setting.
Our approach to the existence of global (in time) smooth solutions is based on a Riemannian approximation scheme: We study graph solutions of the total variation flow in the Riemannian spaces (G, σ ǫ ) where G is a Carnot group and σ ǫ is a family of tame Riemannian metrics that 'collapse' as ǫ → 0 to a sub-Riemannian metric σ 0 in G. The main technical novelties in the paper are a series of different a-priori estimates, which are stable as ǫ → 0: (1) Heat kernel estimates for sub-Laplacians and their elliptic regularizations (see Proposition 2.2).
These estimates are established in the setting of any Carnot group and provide a parabolic counterpart to the time-independent estimates proved by two of the authors in [16]. We remark that (3) and (4) can also be proved, with similar arguments, for solutions of the sub-elliptic mean curvature flow 2 . The limitation to Carnot groups of step two in (4) is due to the fact that for step three and higher there no known suitable barrier functions, or in other words, the class of known explicit minimal graphs is not sufficiently large. In order to state our results we need to introduce some notation.
1.1. Carnot group structure. Let G be an analytic and simply connected Lie group with topological dimension n and such that its Lie algebra G admits a stratification . . , r −1, and [V k , V r ] = 0, k = 1, . . . , r. Such groups are called in [26], [28], and [46] stratified nilpotent Lie groups. Set H(0) = V 1 , and for any x ∈ G we let H(x) = xH(0) = span(X 1 , . . . , X m )(x). The distribution x → H(x) is called the horizontal sub-bundle H. The metric structure is given by assuming that one has a left invariant positive definite form σ 0 defined in the sub-bundle H. We fix a orthonormal horizontal frame X 1 , . . . , X m which we complete it to a basis (X 1 , . . . , X n ) of G by choosing for every k = 2, . . . , r a basis of V k . If X i belongs to V k , then we will define its homogeneous degree as We will denote by xX = n i=1 x i X i a generic element of G. Since the exponential map exp : G → G is a global diffeomorphism we use exponential coordinates in G, and denote x = (x 1 , . . . , x n ) the point exp xX . We also set x H = (x 1 , . . . , x m ). Define non-isotropic dilations as δ s (x) = (s d(i) x i ), for s > 0. We We denote by (X 1 , . . . , X n ) (resp. (X r 1 , . . . , X r n )) the left invariant (resp. right invariant) translation of the frames (X 1 , . . . , X n ).
The vectors X 1 , . . . , X m and their commutators span all the Lie algebra G, and consequently verify Hörmander's finite rank condition ( [32]). This allows to use the results from [39], and define a control distance d 0 (x, y) associated to the distribution X 1 , . . . , X m , which is called the Carnot-Carathéodory metric (denote by d r,0 the corresponding right invariant distance). We call the couple (G, d 0 ) a Carnot Group.
We define a family of left invariant Riemannian metrics σ ǫ , ǫ > 0 in G by requesting that {X ǫ 1 , . . . , X ǫ n } := {X 1 , . . . , X m , ǫX m+1 , . . . , ǫX n } is an orthonormal frame. We will denote by d ǫ the corresponding distance functions. Correspondingly we use ∇ ǫ , (resp. ∇ r ǫ ) to denote the left (resp. right) invariant gradients. In particular, if φ ∈ C ∞ (G) we set We conclude by recalling the expression of the left invariant vector fields in exponential coordinates (see [42]) where p j ik (x) is an homogeneous polynomial of degree k − d(i) and depends only on 1.2. The total variation. The total variation flow is characterized by the fact that each point of the evolving surface graph moves in the direction of the upward unit normal with speed equal to the mean curvature times the volume element. In the setting of the approximating Riemannian metrics (G, σ ǫ ) and in terms of the functions t → u(·, t) : Ω ⊂ G → R describing the evolving graphs, the relevant equation reads: x ∈ Ω and t > 0, with u ǫ (x, 0) = ϕ(x), h ǫ is the mean curvature of the graph of u ǫ (·, t) and for all ξ ∈ R n . In the sub-Riemannian limit ǫ = 0 the equation reads for x ∈ Ω and t > 0, with u(x, 0) = ϕ(x).
We will be concerned with uniform (in the parameter ǫ as ǫ → 0) estimates and with the asymptotic behavior of solutions to the initial value problem for the mean curvature motion of graphs over bounded domains of a Carnot group G, Here ∂ p Q = (Ω × {t = 0}) ∪ (∂Ω × (0, T )) denotes the parabolic boundary of Q. Given appropriate hypothesis on the data, for instance convexity of Ω and ϕ twice continuously differentiable, the classical parabolic theory yields existence and uniqueness for smooth solutions u ǫ of (1.5). However classical parabolic theory will only provide estimates involving constants that degenerate as ǫ → 0 (in the transition from parabolic to degenerate parabolic regime). Our main goal consists in proving stable estimates.
Our first result consists in showing that if the initial/boundary data is sufficiently smooth then the solutions of (1.5) are Lipschitz up to the boundary uniformly in ǫ > 0. Theorem 1.1. (Global gradient bounds) Let G be a Carnot group of step two, Ω ⊂ G a bounded, open, convex 3 set and ϕ ∈ C 2 (Ω). For 1 ≥ ǫ > 0 denote by u ǫ ∈ C 2 (Ω × (0, T )) ∩ C 1 (Ω × (0, T )) the non-negative unique solution of the initial value problem (1.5). There exists C = C(G, ||ϕ|| C 2 (Ω) ) > 0 such that Remark 1.2. The hypothesis of the Theorem above can be slightly weakened by asking that G is a step two Carnot group, ϕ ∈ C 2 (Ω) but that ∂Ω be required to be convex in the Euclidean sense only in a neighborhood of its characteristic locus Σ(∂Ω) and that its mean curvature h ∂Ω ǫ ≤ −δ < 0 in ∂Ω \ Σ(∂Ω). Having established Lipschitz bounds, it is easy to show that the right derivatives X r i u ǫ of the solutions of (1.5) are themselves solutions of (3.1), a divergence form, degenerate parabolic PDE whose weak solutions satisfy a Harnack inequality (see [11] and Proposition 3.8 below). Consequently one obtains C 1,α interior estimates for the solution u ǫ which are uniform in ǫ > 0. At this point one rewrites the PDE in (1.5) in non-divergence form 4 of the total variation flow equation ∂ t u ǫ = a ǫ ij (∇ ǫ u ǫ )X ǫ i X ǫ j u ǫ and invokes the Schauder estimates in Proposition 4.4 to prove local higher regularity and long time existence. Since all the estimates are stable as ǫ → 0 one immediately obtains smoothness and the consequent global in time existence of the solution for the sub-Riemannian case ǫ = 0. Theorem 1.3. In the hypothesis of Theorem 1.1 (or Remark 1.2) one has that there exists a unique solution Since the estimates are uniform in ǫ and in time, and with respect to ǫ, we will deduce the following corollary: Under the assumptions of the Theorem 1.1, as ǫ → 0 the solutions u ǫ converge uniformly (with all its derivatives) on compact subsets of Q to the unique, smooth solution u 0 ∈ C ∞ (Ω × (0, ∞)) ∩ L ∞ ((0, ∞), C 1 (Ω)) of the sub-Riemannian total variation flow (1.4) in Ω × (0, ∞) with initial data ϕ. Corollary 1.5. Under the assumptions of Theorem 1.1, as T → ∞ the solutions u ǫ (·, t) converge uniformly on compact subsets of Ω to the unique solution of the minimal surface equation is the unique solution of the sub-Riemannian minimal surfaces equation h 0 = 0 in Ω, with boundary data ϕ.

Structure stability in the Riemannian limit
The Carnot-Carathéodory metric is equivalent to a more explicitly defined pseudo-distance function, that we will call (improperly) gauge distance, defined as The ball-box theorem in [39] states that there exists A = A(G, σ 0 ) such that for each x ∈ G, If x ∈ G and r > 0, we will denote by the balls in the gauge distance. For each ǫ > 0 we also define the distance function d ǫ corresponding to the Riemannian metric σ ǫ , Note that in the definition of d ǫ , if the curve for which the infimum is achieved happens to be horizontal then d e (x, y) = d 0 (x, y). In general we have sup ǫ>0 d ǫ (x, y) = d 0 (x, y) and it is well known 5 that (G, d ǫ ) converges in the Gromov-Hausdorff sense as ǫ → 0 to the sub-Riemannian space (G, d 0 ). The ball-box theorem in [39] and elementary considerations yield that there exists A = A(G, σ 0 ) > 0 independent of ǫ such that for all x, y ∈ G (see for example [11]).

2.1.
Stability of the homogenous structure as ǫ → 0. If G is a Carnot group, d ǫ is the distance function associated to σ ǫ , we will denote B ǫ (x, r) = {y ∈ G|d ǫ (x, y) < r}. If we denote by dx the Lebesgue measure and by |Ω| the corresponding measure of a subset Ω, then Rea and two of the authors have recently proved in [11] that Proposition 2.1. There is a constant C independent of ǫ such that for every x ∈ G and r > 0, Having this property the spaces (G, d ǫ , dx) are called homogenous with constant C > 0 independent of ǫ (see [20]).
Let τ > 0 and consider the spaceG = G × (0, τ ) with its product Lebesgue measure dxdt. InG define the pseudo-distance function In the paper [11] it is also shown that a Poincaré inequality holds with a choice of a constant which is stable as ǫ → 0. The stability of the homogenous space structure and of the constant in the Poincaré inequality are two of the key factors in the proof of Proposition 3.8.

2.2.
Stability in the estimates of the Heat Kernel. The results in this section are of independent interest and concern uniform Gaussian estimates for the heat kernel associated to elliptic regularizations of the Carnot group sub-Laplacians. They will be used in this paper in connection to the Schauder estimates and the higher regularity of the total variation flow. We will deal with non-divergence form operators similar to those in (1.3), but with constant coefficients. More precisely we will consider the operator ..n is a symmetric, real-valued n × n matrix, such that for some choice of constants Λ, C 1 , C 2 > 0 and for all ξ ∈ R n one has The point here is that the constants C 1 , C 2 , Λ are independent of ǫ and provide Λ−uniform coercivity of (2.5) in the horizontal directions. Formally, in the sub-Riemannian limit ǫ → 0 the equation becomes In order to ensure that the operator L ǫ,A (reps. L A ) is uniformly elliptic (reps. subelliptic), we will assume that the matrix A = (a ij ) belongs to a set of the form A is a symmetric n × n, real valued constant matrix, satisfying (2.6) for some choice of C 1 and C 2 } for some fixed Λ > 0.
Heat kernel estimates in Nilpotent groups are well known (see for instance [33] and references therein). We also refer to [6] where a self contained proof is provided. In our work we will need estimates which are uniform in the variable ǫ as ǫ → 0, in the same spirit as the results in [16]. We will denote Γ ǫ,A (x, t) the fundamental solution of (2.5), with matrix (a ij ) in M Λ , and Γ A (x, t) the fundamental solution of (2.7).

Proposition 2.2.
There exists constants C Λ > 0 depending on G, σ 0 , Λ but independent of ǫ such that for each ǫ > 0, x ∈ G and t > 0 one has For s ∈ N and k−tuple (i 1 , . . . , i k ) ∈ {1, . . . , n} k there exists C s,k > 0 depending only on k, s, G, σ 0 , Λ such that The proof of our result is directly inspired from the proof in [16], where the time independent case was studied by two of us, and where a general procedure was introduced for handling the dependence on ǫ and obtain independent estimates in the sub-Laplacian case.
Let us consider the groupḠ = G × G defined in terms of n new coordinates y ∈ G, so that points ofḠ will be denotedx = (x, y). Denote by Y 1 , . . . , Y n a copy of the vectors X 1 , . . . , X n , defined in terms of new variables y. The vector fields (X i ) and (Y i ) are defined in the product algebraḠ = G × G.
InḠ we will consider two different families of sub-Riemannian structures: • A sub-Riemannian structure determined by the choice of horizontal vector fields given by (2.12) (X 0 1 , · · · ,X 0 m+n ) = (X 1 , . . . , X m , Y 1 , . . . , Y n ). The sub-Laplacian/heat operator associated to this structure is We complete the horizontal frame to a basis of the whole Lie algebra by adding the commutators of the vectors (X i ): (X m+1 , . . . X n ). We denote the exponential coordinates of a pointx around a pointx 0 in terms of the full frame through the coefficients (v 0 i , w 0 i ) which are defined bȳ • For every ǫ ∈ [0, 1) consider a sub-Riemannian structure determined by the choice of horizontal vector fields given by (2.14) (X ǫ 1 , · · · ,X ǫ m+n ) = (X 1 , . . . , X m , Y 1 , . . . , Y m , X ǫ m+1 + Y m+1 , . . . , X ǫ n + Y n ).
The sub-Laplacian/heat operator associated to this structure is Analogously, the exponential coordinates associated to L ǫ,A will depend on the given horizontal frame (2.14) and the family (X m+1 , . . . X n ). The coordinates of a pointx are the coefficients v ǫ i and w ǫ i satisfying We denote byΓ 0,A andΓ ǫ,A the heat kernels of the corresponding heat operators. In both structures we define associated (pseudo)distance functionsd 0 (x,ȳ) andd ǫ (x,ȳ) that are equivalent to those defined in (2.1) and that assign unit weight to the corresponding horizontal vectors while (X m+1 , . . . X n ) will be weighted according to their degree d(i).
We will denote byB ǫ andB 0 the corresponding metric balls.
Proof. In order to estimate the fundamental solution of the operatorsL ǫ,A in terms ofL 0,A , we define a volume preserving change of variables on the Lie algebra and on G: DefineF ǫ :Ḡ →Ḡ asF ǫ (x) = exp(T ǫ (log(x)) with T ǫ (X 0 i ) =X ǫ i for i = 1, . . . , m + n. By definition the relation between the distancesd 0 andd ǫ is expressed by the formula Analogously we also have t), for i = m + 1, · · · , n, with similar identities holding for the iterated derivatives.
In view of the latter, assertions (2.16) and (2.17) immediately follow from the well known estimates ofΓ 0,A (see for instance there references cited above or [6]). The pointwise convergence (2.18) is also an immediate consequence of (2.19) and (2.20). In order to prove the dominated convergence result we need to relate the distancesd 0 andd ǫ : Expressing the exponential coordinates v ǫ i , w ǫ i in terms of v 0 i , w 0 i one easily obtains so that for all 6x ,x 0 ∈Ḡ where C 0 is independent of ǫ. The latter and (2.17) imply dominated convergence onḠ.
Although the vector fieldsX 0 i , i = 1, ..., m + n are not free, we can invoke the Rothschild and Stein lifting theorem [42] and lift them to a new family of free vector fields in a free Carnot group. For the sake of simplicity we will continue to use the same notationX 0 i , i = 1, ..., m + n to denote such family of free vector fields. To recover the desired estimate (2.22) for the original vector fields one needs to argue through projections down to the original space, exactly as done in [39] and [16]. A standard argument (see for instance [7]) yields that for every A ∈ M Λ , if we setĀ as in (2.13) then there exists a Lie group authomorphism F A such that where I denotes the identity matrix. The authomorphism F A is defined by and is extended to the whole Lie algebra as morphism (see the Appendix for more details). As in (2.15), we denote by (v i , w i ) the canonical coordinates of F A1 (x) around F A2 (x), then the Mean Value Theorem yields This estimate indicates the well known fact that at large scale the Riemannian approximating distances are equivalent to the sub-Riemannian distance By the result in the Appendix, the operators v i X i +w i Y i are zero order differential operators whose coefficients can be estimated by ||A 1 − A 2 ||. The conclusion follows by virtue of Proposition 5.1.
We conclude this section with the proof of the main result Proposition 2.2.
Proof of Proposition 2.2. From the definition of fundamental solution we have that for any x ∈ G and t > 0. In view of the (global) dominated convergence of the derivatives ofΓ ǫ,A to the corresponding derivatives ofΓ 0,A as ǫ → 0, we deduce that as ǫ → 0. The Gaussian estimates of Γ ǫ,A follow from the corresponding estimates onΓ ǫ,A and the fact that in view of (2.21),d Indeed the latter shows that there exists a constant C > 0 depending only on G, σ 0 such that for every The conclusion follows at once

Gradient estimates
In this section we prove Theorem 1.1. The proof is carried out in two steps: First we use the maximum principle to establish interior L ∞ bounds for the full gradient of the solution ∇ 1 u of (1.5) with respect to the Lipschitz norm of u on the parabolic boundary. Next, we construct appropriate barriers and invoke the comparison principle established in [8] to prove boundary gradient estimates. The combination of the two will yield the uniform global Lipschitz bounds.

Interior gradient estimates.
Recalling that the right invariant vector fields X r j commute with the left invariant frame X i , i = 1, . . . , n it is easy to show through a direct computation the following result.
. . , n. Then for every h = 0, . . . , n one has that v h is a solution of Note that in order to prove L ∞ bounds on the horizontal gradient of solutions of (1.5) one cannot invoke Lemma 3.1 with differentiation the horizontal frame, because these vector fields do not commute. The argument below bypasses this difficulty by using the algebraic relation between right and left invariant vector fields. Proposition 3.2. Let u ǫ ∈ C 3 (Q) be a solution to (1.5) with Ω bounded. There exists C = C(G, ||ϕ|| C 2 (Ω) ) > 0 such that for every compact subset K ⊂⊂ Ω one has where ∇ 1 is the full σ 1 −Riemannian gradient.
Proof. The proof is similar to the argument in [8,Proposition 5.1] and is based on the trivial observation that left-invariant vector fields at the layer k can be written as linear combination of right-invariant vector fields at level k and left (or right) invariant vector fields at levels k + 1, k + 2, . . . , r. The first step consists in differentiating the equation in (1.5) in the direction of the highest layer V r (which commutes with all the Lie algebra) and invoke the weak maximum principle applied to (3.1) to show (3.2) sup Next one proceeds in a similar fashion, iterating the argument by differentiating the equation in (1.5) along the right invariant vector fields in the layer V r−1 . The maximum principle and (3.1) yield L ∞ bounds on X r i u with d(X r i ) = r − 1. Coupling such bounds with (3.2) one obtains (3.3) sup The proof now follows by induction on the step of the derivatives.

Linear barrier functions.
In [8,Section 4.2] it is shown that in a step two Carnot group coordinate planes (i.e. images under the exponential of level sets of the form x k = 0) solve the minimal surface equation h 0 = 0. In the same paper it is also shown that this may fail for step three or higher. In the construction of the barrier function we will need the following slight refinement of this result, Lemma 3.3. Let G be a step two Carnot group. If f : G → R is linear (in exponential coordinates) then for every ǫ ≥ 0, the matrix with entries X ǫ i X ǫ j f is anti-symmetric, in particular every level set of f satisfies h ǫ = 0.
Proof. We need to show that for f ( We recall the expression (1.2) for the vector fields X i , d(i) = 1 in terms of exponential coordinates 3.3. Boundary gradient estimates. We say that a set Ω ⊂ G is convex in the Euclidean sense if exp −1 (Ω) ⊂ G is convex in the Euclidean space G. In a group of step two this condition is translation invariant.
Lemma 3.5. For each ǫ ≥ 0, if U ǫ is a bounded subsolution and V ǫ is a bounded supersolution of (1.5) then Let u ǫ ∈ C 2 (Q) be a solution of (1.5), and express the evolution PDE in non-divergence form Set v ǫ = u ǫ − ϕ so that v ǫ solves the homogenous 'boundary' value problem . We define our (weakly) parabolic operator for which the function v ǫ is a solution In the following we construct for each point p 0 = (x 0 , t 0 ) ∈ ∂Ω × (0, T ) a barrier function for Q, v ǫ : i.e., Lemma 3.6. Let G be a Carnot group of step two and Ω ⊂ G convex in the Euclidean sense. For each point p 0 = (x 0 , t 0 ) ∈ ∂Ω × (0, T ) one can construct a positive function w ∈ C 2 (Q) such that Proof. In the hypothesis that Ω is convex in the Euclidean sense we have that every x 0 ∈ ∂Ω supports a tangent hyperplane P defined by an equation of the form Π(x) = n i=1 a i x i = 0 with Π > 0 in Ω, Π(x 0 ) = 0, and normalized as d(i)=1,2 a 2 i = 1. Following the standard argument (see for instance [34, Chapter 10]) we select the barrier at (x 0 , t 0 ) ∈ ∂Ω × (0, T ) independent of time with with k and ν chosen appropriately so that conditions (3.11) will hold. We choose a neighborhood V = O × (0, T ) such that P ∩ O ∩ ∂Ω = {x 0 }. By an appropriate choice of k sufficiently large we can easily obtain w ǫ (p 0 ) = 0 and w ǫ ≥ v ǫ in ∂ p V ∩ Q.
To verify Q(w ǫ ) ≤ 0 we begin by observing that w satisfies with F = a ǫ ij (∇ e w + ∇ ǫ ϕ)X ǫ i wX ǫ j w. We will show: in a parabolic neighborhood of p 0 . The first claim holds with an equality as a ǫ ij is symmetric and X ǫ i X ǫ j Π is anti-symmetric in view of Lemma 3.3. To establish (Claim 2) we first note that Lemma 3.3 implies for some constant C(G) > 0. Consequently, for Φ ′ >> 1 sufficiently large one finds with C(G) > 0 a constant depending only on G (not always the same along the chain of inequalities). In view of the definition of b ǫ and (3.13) with an appropriate choice of ν = ν(G, ǫ, φ) > 0 and k = k(G, φ) >> 1 in (3.14), we conclude In view of Lemma 3.5, a comparison with the barrier constructed above yields that in V ∩ Q, with dist σ1 (x, x 0 ) being the distance between x and x 0 in the Riemannian metric σ 1 , concluding the proof of the boundary gradient estimates.
The proof of Theorem 1.1 now follows immediately from Proposition 3.2 and Proposition 3.4.
Having established uniform global Lipschitz bounds one now notes that equation (3.1) satisfies horizontal coercivity conditions uniformly in ǫ > 0. Such conditions are among the main hypothesis of the Harnack inequality in [11]. In this paper, G. Rea and two of the authors have proved that given a homogenous structure and a Poincaré inequality, then a sub-elliptic analogue of Aronson and Serrin's Harnack inequality [3] for quasilinear parabolic equations holds. As a consequence one obtains for some α ∈ (0, 1) that the solutions u ǫ to (1.5) satisfy C 1,α Hölder estimates, uniform in ǫ ∈ (0, 1), Definition 3.7. Let 0 < α < 1, Q ⊂ R n+1 and u be defined on Q. We say that u ∈ C α ǫ,X (Q) if there exists a positive constant M such that for every (x, t), , t), (x 0 , t 0 )). We put Iterating this definition, if k ≥ 1 we say that u ∈ C k,α ǫ,X (Q) if for all i = 1, . . . , m X i ∈ C k−1,α ǫ,X (Q). Where we have set C 0,α ǫ,X (Q) = C α ǫ,X (Q).

Regularity properties in the C k,α spaces
In this section we will prove uniform estimates for solution of (1.3) in the C k,α ǫ,X Hölder spaces. This is accomplished by using the uniform Gaussian bounds established in Section 2.2 to develop new uniform Schauder estimates for solutions of second order sub-elliptic differential equations in non divergence form in a cylinder Q = Ω × (0, T ) that are stable as ǫ → 0. As usual we will make use of the associated linear, constant coefficient frozen operator: where (x 0 , t 0 ) ∈ Q.
As a direct consequence of the definition of fundamental solution one has the following representation formula where we have denoted by Γ ǫ (x0,t0) the heat kernel for of L ǫ,(x0,t0) . We explicitly note that for ǫ > 0 fixed the operator L ǫ,(x0,t0) is uniformly parabolic. iIs heat kernel can be studied through standard singular integrals theory in the corresponding Riemannian balls. Hence, as noted in [29], one is allowed to differentiate twice the kernels defined in (4.1) with respect to any right or left invariant vector field. Proposition 4.2. Let 0 < α < 1 and w be a smooth solution of L ǫ w = f ∈ C α ǫ.X (Q) in the cylinder Q. Let K be a compact sets such that K ⊂⊂ Q, set 2δ = d 0 (K, ∂ p Q) and denote by K δ the δ−tubular neighborhood of K. Assume that there exists a constant C > 0 such that for every ǫ ∈ (0, 1) There exists a constant C 1 > 0 depending on δ, α, C and the constants in Proposition 2.2 such that The proof follows the outline of the standard case, as in [29,Theorem 4,Chapter 3], and rests crucially on the Gaussian estimates proved in Proposition 2.2, and the fact that the functions L ǫ,(x0,t0) − L ǫ (w φ), and f φ + wL ǫ φ + 2a ǫ ij X ǫ i wX ǫ j φ are Hölder continuous. Choose a parabolic sphere 7 B ǫ,δ ⊂⊂ K where δ > 0 will be fixed later and a cut-off function φ ∈ C ∞ 0 (R n+1 ) identically 1 on B ǫ,δ/2 and compactly supported in B ǫ,δ . This clearly implies that for some constant C > 0 depending only on G and σ 0 , |∇ ǫ φ| ≤ Cδ −1 , |L ǫ φ| ≤ Cδ −2 , in Q. Next, invoking (4.1) one has that for every multi-index I = (i 1 , i 2 ) ∈ {1, . . . , m} 2 and for every ( The uniform Hölder continuity of a ǫ ij , with Proposition 2.2 and Lemma 2.5 yield with C > 0 independent of ǫ. In view of the latter, using basic singular integral properties (see [26]) and proceeding as in [29, Theorem 2, Chapter 4], we obtain , where C 1 , andC 1 are stable as ǫ → 0. Similarly, if φ is fixed, the Hölder norm of the second term in the representation formula (4.2) is bounded by From (4.2), (4.3) and (4.4) we deduce that Choosing δ sufficiently small we prove the assertion on the fixed sphere B ǫ,δ The conclusion follows from a standard covering argument.
A direct calculation shows the following Lemma 4.3. Let k ∈ N and consider a k−tuple (i 1 , . . . , i k ) ∈ {1, . . . , m} k . There exists a differential operator B of order k + 1, depending on horizontal derivatives of a ǫ ij of order at most k, such that Proposition 4.4. Let w be a smooth solution of L ǫ w = f on Q. Let K be a compact sets such that K ⊂⊂ Q, set 2δ = d 0 (K, ∂ p Q) and denote by K δ the δ−tubular neighborhood of K.Assume that there exists a constant C > 0 such that ||a ǫ ij || C k,α ǫ,X (K δ ) ≤ C, for any ǫ ∈ (0, 1). There exists a constant C 1 > 0 depending on α, C, δ, and the constants in Proposition 2.2, but independent of ǫ, such that Proof. The proof is similar to the previous one for k = 1. We note that differentiating the representation formula (4.1) along any k−tuple (i 1 , . . . , i k ) ∈ {1, . . . , m} k , and using Lemma 4.3, yields where φ is as in the proof of Proposition 4.2 and B is a differential operator of order k + 1. The conclusion follows by further differentiating the previous representation formula along two horizontal directions X ǫ j1 X ǫ j2 and arguing as in the proof of Proposition 4.2.
Proof of Theorem 1.3. We will prove by induction that for every k ∈ N and for every compact set K ⊂⊂ Q there exists a positive constant C such that (4.5) ||u ǫ || C k,α ǫ,X (K) ≤ C, for every ǫ > 0. The assertion is true if k = 2, by Corollary 3.8 and Proposition 4.2. If (4.5) is true for k + 1, then the coefficients in (1.3) satisfy a ǫ ij ∈ C k,α ǫ,X uniformly as ǫ ∈ (0, 1) and (4.5) thus holds at the level k + 2 by virtue of Proposition 4.4.

Appendix
Let us assume that G be a free nilpotent Lie group, and let G = V 1 ⊕ ... ⊕ V r its associated Lie algebra. Denote by X 1 , . . . , X m a basis of the first layer of the Lie algebra. and by X m+1 , . . . , X n a list of vectors which complete a basis of the tangent space. If u(x) is an homogeneous polynomial of order p and X is a differential operator of degree q we will call u(x)X differential operator of degree q − p. If a ij = (A) ij is a real valued, constant coefficient, m × m positive definite matrix , one can define a Lie algebra automorphism a ij X j for every i = 1, . . . , m and extend it to the whole Lie algebra as a morphism Note that T A can be represented as a block matrixĀ = A 1 ⊕ A 2 · · · ⊕ A r , with A 1 = A and where the block A j act on vectors of degree j, and its coefficients are polynomial of order j of the coefficients of A 1 . In particular T A (X) and X have the same degree. Via the exponential map T A induces a group automorphism on the whole group F A : G → G, defined by F A (x) = exp(T A (log(x)).
In terms of exponential coordinates (around the origin) where d(i) denotes the homogenous degree defined as in (1.1). Since v i X i is a zero order homogeneous operator then also v i T A (X i ) is a differential operator of order zero.
Proposition 5.1. Let (w i ) i=1,...,n denote the canonical coordinates of F A (x) around F B (x), and (v j ) j=1,...,n the coordinates of x around 0, as in (5.1). There exists M ∈ N depending only on the group structure; constants c 1 , . . . , c M depending only on the group structure; zero-order differential operators Y 1 , . . . , Y M depending only on the group structure, on the coefficients (v j ) j=1,...,n and their derivatives along vector fields up to order r − 1 (but independent on A and B), and a constant C = C(||A||, ||B||, G) > 0 such that c l Y l and |c l | ≤ C||A − B|| for l = 1, . . . , M.
Proof. By definition n j=1 where * is the Baker-Campbell-Hausdorff operation (see for instance [42] and references therein). This formula allows to represent the left-hand side as a finite 8 sum of terms involving commutators (up to order r) of the two terms on the right-hand side. Let us consider separately the terms in this representation starting with those involving no commutators, i.e.
Clearly this expression is in the same form as (5.2). Next we consider the second term in the Baker-Campbell-Hausdorff sum, i.e. that involving commutators of the form It is evident that the coefficients of the commutators above are Lipschitz functions of A and B, vanishing for A = B, and consequently the whole sum satisfies (5.2). To conclude the proof we note that all further terms of the Baker-Campbell-Hausdorff sum involve operators as in (5.2) and zero order operators of the form . Such commutators give rise to differential operators of zero order whose coefficients are polynomials in A and B and vanish when A = B.