\section{Introduction}\label{NPD}
Let $M$ be a connected $C^\infty $ manifold of dimension $d$ and let $a$ be a $C^\infty $ non-negative quadratic form on the cotangent space $T^*M$.
We assume that $a$ has a sub-Riemannian structure, that is to say, there exist $m\in{\mathbb N}$ and $C^\infty $ vector fields $X_1,\dots,X_m$ on $M$ such that
\begin{equation}\label{LAC}
a(\xi,\xi)=\<\xi,a(x)\xi\>=\sum_{\ell=1}^m\<\xi,X_\ell(x)\>^2,\quad \xi\in T_x^*M
\end{equation}
and such that 
\begin{equation}\label{LAD}
\operatorname{span}\{Y(x):Y\in{\cal A}(X_1,\dots,X_m)\}=T_xM,\quad\text{for all $x\in M$}.
\end{equation}
Here ${\cal A}(X_1,\dots,X_m)$ denotes the smallest set of vector fields on $M$ containing $X_1,\dots,X_m$ and closed under the commutator product, given by
$$
[X,Y](x)=\sum_{i=1}^dX^i(x)\frac{\partial Y}{\partial x^i}(x)-Y^i(x)\frac{\partial X}{\partial x^i}(x).
$$
There is associated to the quadratic form $a$ an energy function $I$ on the set of continuous paths ${\Omega}=C([0,1],M)$.
An absolutely continuous path ${\omega}\in{\Omega}$ may have a driving path $\xi$, 
by which we mean a measurable path $\xi$ in $T^*M$ such that $\xi_t\in T^*_{{\omega}_t}M$ and $\dot{\omega}_t=a({\omega}_t)\xi_t$ for almost all $t$. 
Then ${\omega}$ has energy
$$
I({\omega})=\int_0^1\<\xi_t,a({\omega}_t)\xi_t\>dt.
$$
If ${\omega}$ is not absolutely continuous or has no driving path, then $I({\omega})=\infty$. 
Write $H$ for the subset of ${\Omega}$ consisting of paths of finite energy.
For $x,y\in M$, set 
$$
H^{x,y}=\{{\omega}\in H:{\omega}_0=x,{\omega}_1=y\}.
$$
It is well known that $H^{x,y}$ is non-empty under the bracket condition \eqref{LAD} and that the sub-Riemannian distance
\begin{equation}\label{DIST}
d(x,y)=\inf_{{\omega}\inH^{x,y}}\sqrt{I({\omega})}
\end{equation}
defines a metric compatible with the topology of $M$.

In this paper we investigate the restriction of the energy function $I$ to $H^{x,y}$ and its relation to heat flow from $x$ to $y$.
We make no assumption of completeness for $M$ so in general the infimum in \eqref{DIST} may not be achieved.
Nevertheless, we will show in Theorem \ref{VAR} that the small-time logarithmic asymptotics for heat kernels are determined by the distance function 
and indeed that the associated diffusion measures concentrate exponentially fast in small time around paths of nearly minimal energy.

The main work in the paper addresses the case where $x$ and $y$ are chosen so that $I$ achieves a minimum on $H^{x,y}$ uniquely, say at ${\gamma}$.
We will then construct, under a regularity condition on ${\gamma}$, a vector space ${T_\g H^{x,y}}$ of absolutely continuous paths $v$ in $TM$, with $v_t\inT_{\g_t}M$ for all $t$ and $v_0=v_1=0$,
along with an equivalence class of norms on ${T_\g H^{x,y}}$, each making ${T_\g H^{x,y}}$ into a Hilbert space.
The paths in ${T_\g H^{x,y}}$ can be thought of as the infinitesimal variations of ${\gamma}$ in $H^{x,y}$.
We will further construct a continuous non-negative quadratic form $Q$ on ${T_\g H^{x,y}}$ such that $Q(v)$ is the minimal second variation of $I$ in the direction $v$, in a sense to be made precise.
These constructions are the content of Section \ref{SUB}.

The sub-Riemannian cut locus was defined by Bismut \cite{B} in terms of the bicharacteristic flow associated to $a$.
An alternative characterization was then given by Ben Arous \cite{BA} in terms of a quadratic form on control paths associated to a compatible sub-Riemannian structure.
We will show in Theorem \ref{2.7} that $(x,y)$ lies outside the sub-Riemannian cut locus if and only if our regularity condition on ${\gamma}$ holds and $Q$ is positive-definite.
Since the vector space ${T_\g H^{x,y}}$ and the quadratic form $Q$ are intrinsic to $a$, 
this provides an intrinsic characterization of the cut locus in terms of the energy function in this general setting.
To clarify the distinction between objects intrinsic to $a$ and those constructed from a compatible sub-Riemannian structure, 
we show by example in Section \ref{SUB} that, when $a$ has non-constant rank, it may admit two inequivalent sub-Riemannian structures.
For control theory, the sub-Riemannian structure is the basic given quantity, but for heat flow it is natural to seek objects intrinsic to $a$.

The sub-Riemannian cut locus is less well understood than its Riemannian counterpart. 
It is known to be a closed and symmetric subset of $M\times M$. 
Rifford \& Tr\'elat \cite{MR2139255} and Agrachev \cite{MR2513150} have proved results which limit its size.
The reader may find in the lecture notes of Agrachev, Barilari \& Boscain \cite[Theorems 10.4 and 10.11]{ABB} a proof that, for all $x\in M$ and any $r\in(0,\infty)$ such that $B=\{y\in M:d(x,y)\leqslant r\}$ is compact, 
the set of points $y\in B$ such that $(x,y)$ lies outside the cut locus is dense in $B$.

Let ${\cal L}$ be a second order differential operator on $M$ with $C^\infty $ coefficients, such that $\cL1=0$ and such that ${\cal L}$ has principal symbol $a/2$.
In each coordinate chart, ${\cal L}$ thus takes the form
\begin{equation}\label{LAB}
{\cal L}=\frac12\sum_{i,j=1}^da^{ij}(x)\frac{\partial^2}{\partial x^i\partial x^j}+\sum_{i=1}^db^i(x)\frac{\partial}{\partial x^i}
\end{equation}
for some $C^\infty $ functions $b^i$ on $M$.
In this context we refer to $a$ as the diffusivity.
Given that $a$ has sub-Riemannian structure $(X_1,\dots,X_m)$, 
we obtain the same family of operators by choosing another $C^\infty $ vector field $X_0$ on $M$ and setting\footnote{%
Here, $X_\ell$ is used to denote the differential operator $\sum_{i=1}^dX^i_\ell(x)\partial/\partial x^i$.}
\begin{equation}\label{HSS}
{\cal L}=\frac12\sum_{\ell=1}^mX_\ell^2+X_0.
\end{equation}
We will sometimes impose the following global condition
\begin{equation}\label{URD}
M={\mathbb R}^d\quad\text{and}\quad X_0,X_1,\dots,X_m\in C^\infty _b({\mathbb R}^d,{\mathbb R}^d).
\end{equation}
At other times, we will impose no further condition on $M$ or $a$, in particular no condition of completeness, but then we will insist that there is a $C^\infty $ $1$-form ${\beta}$ satisfying the sector condition
\begin{equation}\label{SECT}
\|a({\beta},{\beta})\|_\infty<\infty
\end{equation}
and there is a locally invariant\footnote{%
When $M$ is compact, the local invariance property does not restrict the class of operators. 
To see this, suppose that ${\cal L}$ has the form \eqref{CFO} for some positive $C^\infty $ measure $\nu$, not necessarily invariant.
Then ${\cal L}$ has a positive $C^\infty $ invariant density function $\rho$ with respect to $\nu$ on $M$.
Thus, for all $C^\infty $ functions $f$, we have.
$$
\int_Ma({\beta}\rho-\tfrac12\nabla\rho,\nabla f)d\nu=\int_M{\cal L} f\rho d\nu=0.
$$
Set $\tilde{\beta}={\beta}-\frac12\nabla\log\rho$ and $\tilde\nu=\rho\nu$. 
Then ${\cal L} f=\frac12\widetilde{\operatorname{div}}(a\nabla f)+a(\tilde{\beta},\nabla f)$, $\tilde{\beta}$ satisfies \eqref{SECT}, and $\tilde\nu$ is invariant.
In particular, when $M$ is compact and $a$ is positive-definite, we can write all operators ${\cal L}$ of the form \eqref{LAB} in the form \eqref{CFO}.}
positive $C^\infty $ measure $\nu$ on $M$ such that
\begin{equation}\label{CFO}
{\cal L} f=\tfrac12\operatorname{div}(a\nabla f)+a({\beta},\nabla f)
\end{equation}
where the divergence is understood with respect to $\nu$.
Here, by a positive $C^\infty $ measure $\nu$ we mean any measure having a positive $C^\infty $ density with respect to Lebesgue measure in each coordinate chart.
We say that $\nu$ is locally invariant for ${\cal L}$ if, for all $C^\infty $ functions $f$ of compact support in $M$, we have
\begin{equation}\label{INV}
\int_Ma({\beta},\nabla f)d\nu=0.
\end{equation}
The divergence of a vector field $X$ with respect to $\nu$ is characterized by 
$$
\int_Mf\operatorname{div} Xd\nu=-\int_MXfd\nu
$$
for all such functions $f$.
From \eqref{CFO}, we find that $X_0=\sum_{\ell=1}^m{\alpha}_\ell X_\ell$, where ${\alpha}_\ell=\frac12\operatorname{div} X_\ell+\<{\beta},X_\ell\>$.
In particular, $X_0(x)\in\operatorname{span}\{X_1(x),\dots,X_m(x)\}$ for all $x\in M$, which is not required in \eqref{URD}.
The sector condition and the local invariance property are satisfied trivially in the case ${\beta}=0$.
Thus our results apply to all sub-Riemannian Laplacians.

There is a family of probability measures on the set $\Omega^{x,y}=\{{\omega}\in{\Omega}:{\omega}_0=x,{\omega}_1=y\}$ which is naturally associated to the operator ${\cal L}$.
Fix ${\varepsilon}>0$ and $x\in M$. 
There exists a diffusion process starting from $x$ and having generator ${\varepsilon}{\cal L}$.
Since, in general, the coefficients of ${\cal L}$ may be unbounded and we make no assumption of completeness for $M$, this diffusion may explode with positive probability, that is to say it may leave all compact sets in finite time.
We will write $\mu^x_\ve$ for the unique sub-probability measure on ${\Omega}$ which is the law of this diffusion restricted to paths which do not explode by time $1$.
Under our assumptions, there is a unique family of probability measures $(\mu^{x,y}_\ve:y\in M)$ on ${\Omega}$ which is weakly continuous in $y$, 
with $\mu^{x,y}_\ve$ supported on $\Omega^{x,y}$ for all $y$, and such that
$$
\mu^x_\ve(d{\omega})=\int_M\mu^{x,y}_\ve(d{\omega})p({\varepsilon},x,dy)
$$
where $p({\varepsilon},x,.)$ is the (sub-)law of ${\omega}_1$ under $\mu^x_\ve$.
More explicitly, the finite-dimensional distributions of each measure $\mu^{x,y}_\ve$ may be written as follows.
Choose a positive $C^\infty $ measure $\nu$ on $M$.
It does not have to be the same measure as in the preceding paragraph.
There is a positive $C^\infty $ function $p$ on $(0,\infty)\times M\times M$ such that
$$
p({\varepsilon},x,dy)=p({\varepsilon},x,y)\nu(dy).
$$
This function $p$ is the Dirichlet heat kernel for ${\cal L}$ with respect to $\nu$.
Then, for all $k\in {\mathbb N}$, all $t_1,\dots,t_k\in(0,1)$ with $t_1<t_2<\dots<t_k$ and all $x_1,\ldots,x_k\in M$, we have
\begin{align*}
&\mu^{x,y}_\ve(\{{\omega}:{\omega}_{t_1}\in dx_1,\ldots,{\omega}_{t_k}\in dx_k\})\\
&\quad\q=\frac{p({\varepsilon} t_1,x,x_1)p({\varepsilon}(t_2-t_1),x_1,x_2)\ldots p({\varepsilon}(1-t_k),x_k,y)}{p({\varepsilon},x,y)}
\nu(dx_1)\dots\nu(dx_k).
\end{align*}
It is straightforward to see that these finite-dimensional distributions are consistent and do not depend on the choice of $\nu$.

Our main aim is to understand the behaviour of the diffusion bridge measures $\mu^{x,y}_\ve$ in the limit ${\varepsilon}\to0$.
We will say that ${\gamma}\inH^{x,y}$ is strongly minimal if there exist ${\delta}>0$ and a relatively compact open set $U$ in $M$ such that 
$$
I({\omega})\geqslant I({\gamma})\text{ for all ${\omega}\inH^{x,y}$}\quad\text{and}\quad I({\omega})\geqslant I({\gamma})+{\delta}\text{ for all ${\omega}\inH^{x,y}$ which leave $U$}.
$$
When $M$ is complete for the sub-Riemannian distance, the second condition follows from the first by the Hopf--Rinow theorem.
Moreover, if ${\gamma}$ is the only strongly minimal path in $H^{x,y}$, 
then the second condition holds for all relatively compact domains $U$ containing ${\gamma}$ by a weak compactness argument.
Here is our first main result.
\begin{theorem}\label{VAR}
Let $M$ be a connected $C^\infty $ manifold.
Let ${\cal L}$ be a second order differential operator on $M$ of the form \eqref{CFO}, 
$$
{\cal L} f=\tfrac12\operatorname{div}(a\nabla f)+a({\beta},\nabla f)
$$
where the diffusivity $a$ has a sub-Riemannian structure, 
where the divergence is taken with respect to a positive $C^\infty $ locally invariant measure, and where ${\beta}$ is a $C^\infty $ $1$-form satisfying the sector condition \eqref{SECT}.
Then, for all $x,y\in M$, as $t\to0$, we have
\begin{equation}\label{LOGP}
t\log p(t,x,y)\to-d(x,y)^2/2.
\end{equation}
Moreover, for all ${\delta}>0$, for $r={\delta}^{1/4}(d(x,y)^2+{\delta})^{1/2}$, we have
\begin{equation}\label{LDEV}
\limsup_{{\varepsilon}\to0}{\varepsilon}\log\mu^{x,y}_\ve(\{{\omega}\in\Omega^{x,y}:d({\omega}_t,{\Gamma}_t({\delta}))\geqslant r\text{ for some }t\in[0,1]\})\leqslant-{\delta}/2
\end{equation}
where
$$
{\Gamma}_t({\delta})=\{{\gamma}_t:{\gamma}\inH^{x,y},\, I({\gamma})\leqslant d(x,y)^2+{\delta}\}.
$$
In particular, in the case where there is a unique strongly minimal path ${\gamma}\inH^{x,y}$, we have $\mu^{x,y}_\ve\to\d_\g$ as ${\varepsilon}\to0$ weakly on $\Omega^{x,y}$, 
where $\d_\g$ is the unit mass at ${\gamma}$.
\end{theorem}
The small-time logarithmic asymptotics \eqref{LOGP} were proved by Varadhan \cite{MR0208191} in the case when $M={\mathbb R}^d$ and $a$ is everywhere positive-definite.
They were generalized by L\'eandre \cite{MR871256,MR904825} to the sub-Riemannian case, under the global assumption \eqref{URD}, using probabilistic techniques which make this assumption unavoidable.
While the lower bound extends easily, for the upper bound it appears that analytic techniques are needed, and we arrive at \eqref{LOGP} by checking the applicability of a general result of Sturm \cite{MR1355744} in this context.
In the case of Brownian motion in a complete Riemannian manifold, the exponential estimate \eqref{LDEV} follows from a result of Hsu \cite{MR1027823}. 
Under the assumption \eqref{URD} and subject to the condition that $a(x)$ is positive-definite, \eqref{LDEV} follows from work of Inahama \cite{INA}. 

We now go beyond the law of large numbers and exponential asymptotics contained in Theorem \ref{VAR} to identify a small-time equivalent for the heat kernel
and to obtain a central limit theorem for the diffusion bridge measures.
For this we will assume that $x$ and $y$ are chosen so that
\begin{equation}\label{UPM}
\text{there is a unique strongly minimal path ${\gamma}\inH^{x,y}$}
\end{equation}
and further that $(x,y)$ lies outside the sub-Riemannian cut locus. 
We first state a small-time equivalent for the heat kernel, which was proved by Ben Arous \cite{BA} under the condition \eqref{URD}.
In fact, Ben Arous obtained a full asymptotic expansion, with local regularity statements, but we will be content with the basic form.

\begin{theorem}\label{BALOC}
Under the conditions of Theorem \ref{VAR}, in the case where there is a unique strongly minimal path ${\gamma}\inH^{x,y}$ and $(x,y)$ lies outside the cut locus, there is a constant $c(x,y)>0$ such that, in the limit $t\to0$,
$$
p(t,x,y)=c(x,y)t^{-d/2}\operatorname{exp}\{-d(x,y)^2/(2t)\}(1+o(1)).
$$
\end{theorem}

We turn to the second order result which is the true focus of this work.
The quadratic form $Q$ mentioned above allows us to give ${T_\g H^{x,y}}$ an intrinsic Hilbert norm.
Write ${T_\g\oxy}$ for the set of continuous paths $v$ in $TM$ such that $v_t\inT_{\g_t}M$ for all $t$ and $v_0=v_1=0$.
We make ${T_\g\oxy}$ into a Banach space using the uniform norm $\|v\|_\infty=\sup_{t\in[0,1]}|v_t|$ corresponding to a choice of Riemannian metric on $M$.
The associated topology on ${T_\g\oxy}$, which is all that matters for us, does not depend on the choice of metric.
We will show in Theorem \ref{SRMR} that there is a unique zero-mean Gaussian measure $\mu_\g$ on ${T_\g\oxy}$ such that
$$
\int_{T_\g\oxy}\phi(v)^2\mu_\g(dv)=Q(\tilde\phi)
$$
for all continuous linear functionals $\phi$ on ${T_\g\oxy}$, where $\tilde\phi\in{T_\g H^{x,y}}$ is determined by $\phi(v)=Q(\tilde\phi,v)$ for all $v\in{T_\g H^{x,y}}$.

We rescale the fluctuations of the diffusion bridge around the minimal path ${\gamma}$ to obtain a non-degenerate limit.
To do this, we choose a $C^\infty $ map ${\theta}:M\to{\mathbb R}^d$ such that the derivative ${\theta}^*({\gamma}_t):T_{\g_t}M\to{\mathbb R}^d$ is invertible for all $t\in[0,1]$.
We can always choose a chart $U$ along ${\gamma}$ and obtain a suitable function ${\theta}$ by extending the coordinate map from a neighbourhood of ${\gamma}$.
Define ${\sigma}_{\varepsilon}:\Omega^{x,y}\to{T_\g\oxy}$ by
$$
{\sigma}_{\varepsilon}({\omega})_t={\theta}^*({\gamma}_t)^{-1}({\theta}({\omega}_t)-{\theta}({\gamma}_t))/\sqrt{\varepsilon}.
$$
Then we obtain a probability measure $\tilde\mu^{x,y}_\ve$ on ${T_\g\oxy}$ by setting
$$
\tilde\mu^{x,y}_\ve=\mu^{x,y}_\ve\circ{\sigma}_{\varepsilon}^{-1}.
$$
It is straightforward to check that if $\tilde\mu^{x,y}_\ve$ converges weakly on ${T_\g\oxy}$ as ${\varepsilon}\to0$ for one choice of the function ${\theta}$, 
then it does so for all such choices and with the same limit.

\begin{theorem}\label{MR}
Let $M$ be a connected $C^\infty $ manifold.
Let ${\cal L}$ be a second order differential operator on $M$, which has the form \eqref{LAB} in each coordinate chart, and whose diffusivity has a sub-Riemannian structure.
Assume that either \eqref{HSS} and \eqref{URD} hold or \eqref{SECT} and \eqref{CFO} hold.
Let $x,y\in M$.
Suppose that there is a unique strongly minimal path ${\gamma}\inH^{x,y}$ and that $(x,y)$ lies outside the cut locus.
Then the rescaled diffusion bridge measures satisfy $\tilde\mu^{x,y}_\ve\to\mu_\g$ weakly on ${T_\g\oxy}$ as ${\varepsilon}\to0$.
\end{theorem}

In the case of Brownian motion on a compact Riemannian manifold, the convergence $\tilde\mu^{x,y}_\ve\to\mu_\g$, with $\mu_\g$ characterized by its covariance, 
was derived by Molchanov \cite{Mo}, but a full proof was not given.
Many of the techniques we use follow ideas pioneered by Bismut \cite{B} and Ben Arous \cite{BA} in their studies of the heat kernel, indeed, it turns out that the diffusion bridge asymptotics are another side of the same story.
Two recent works have taken the fundamental ideas of Bismut and Ben Arous in different directions to ours. 
First, Deuschel, Friz, Jacquier \& Violante \cite{MR3139426,MR3149845} have obtained small time expansions for marginal distributions of the heat flow.
Second, Barilari, Boscain \& Neel \cite{MR3005058} have found estimates of the heat kernel actually on the cut locus.

In Section \ref{SRN}, we recall the notions of the bicharacteristic flow associated to $a$ and the sub-Riemannian cut locus.
We also give a characterization of the limit fluctuation measure $\mu_\g$ in terms of the bicharacteristic flow.
Then, in Section \ref{RIEM}, we discuss the case where $a(x)$ is positive-definite for all $x$, and thus defines a Riemannian metric. 
In this case the limit measure $\mu_\g$ has several different characterizations using classical objects of Riemannian geometry. 
The proof of Theorem \ref{MR} is given under the global assumption \eqref{URD} in Section \ref{LAP}, 
adapting a method of Azencott, Bismut and Ben Arous, and relying on ideas of Malliavin calculus.
Then in Section \ref{LL}, we show how to localize to a neighbourhood of the minimal path, 
and we prove Theorem \ref{BALOC} and complete the proof of Theorem \ref{MR} using analytic estimates in Section \ref{NFB}.
Section \ref{NFB} also contains the proof of Theorem \ref{VAR} and does not depend on material in other sections.
As outlined above, Sections \ref{SUB} and \ref{QG} are devoted to the geometry of some spaces of paths, leading to the intrinsic construction of the minimal second
variation $Q$ of the energy function and its associated Gaussian measure. 
These two sections may be read independently of the preceding probabilistic and analytic parts.

We are grateful to Emmanuel Tr\'elat and Fabrice Baudoin for helpful advice in the course of this work.

\section{Cut locus and fluctuation measure in terms of the bicharacteristic flow}\label{SRN}
The bicharacteristic flow is the maximal flow $(\psi_t({\lambda}):{\lambda}\in T^*M,t\in({\zeta}^-({\lambda}),{\zeta}^+({\lambda})))$ of the vector field $V$ on $T^*M$
given by ${\beta}(V,.)=d{\cal H}$, where ${\beta}$ is the canonical symplectic $2$-form on $T^*M$ and ${\cal H}:T^*M\to[0,\infty)$ is the Hamiltonian
$$
{\cal H}({\lambda})=\tfrac12\<{\lambda},a(x){\lambda}\>,\quad {\lambda}\in T_x^*M.
$$
Thus, for all ${\lambda}\in T^*M$, we have $\psi_0({\lambda})={\lambda}$ and ${\zeta}^-({\lambda})<0<{\zeta}^+({\lambda})$, and
$$
\dot\psi_t({\lambda})=V(\psi_t({\lambda})),\quad t\in({\zeta}^-({\lambda}),{\zeta}^+({\lambda}))
$$
and $\psi_t({\lambda})$ leaves all compact sets in $T^*M$ as $t\to{\zeta}^+({\lambda})$ if ${\zeta}^+({\lambda})<\infty$ and as $t\to{\zeta}^-({\lambda})$ if ${\zeta}^-({\lambda})>-\infty$.
If $M$ is complete, then ${\zeta}^\pm({\lambda})=\pm\infty$ for all ${\lambda}\inT^* M$.
In a coordinate chart $U$ for $M$, for ${\lambda}\in U$, write ${\lambda}_t=\psi_t({\lambda})=(x_t,p_t)$ while ${\lambda}_t$ remains in $T^*U$, with $x_t\in U$ and $p_t\in{\mathbb R}^d$.
Then, for any sub-Riemannian structure $(X_1,\dots,X_m)$ for $a$, we have
\begin{align}\notag
\dot x_t&=\sum_{\ell=1}^m\<p_t,X_\ell(x_t)\>X_\ell(x_t),\quad x_0=x\\
\label{BCF}
\dot p_t&=-\sum_{\ell=1}^m\<p_t,X_\ell(x_t)\>\<p_t,\nabla X_\ell(x_t)\>,\quad p_0=p
\end{align}
where $\<p,\nabla X(x)\>_i=(\partial/\partial x^i)\<p,X(x)\>$.

Fix $x,y\in M$ and suppose that condition \eqref{UPM} holds, that is to say, there is a unique strongly minimal path ${\gamma}\inH^{x,y}$.
Write $\pi$ for the projection $T^*M\to M$.
We assume that ${\gamma}$ is a normal minimizer, meaning that
\begin{equation}\label{PBC}
\text{there exists a bicharacteristic $({\lambda}_t)_{t\in[0,1]}$ such that ${\gamma}_t=\pi{\lambda}_t$ for all $t\in[0,1]$.}
\end{equation}
For $t\in[0,1]$, write ${\lambda}_t=\psi_t({\lambda}_0)$ and define linear maps $J_t:T^*_xM\toT_{\g_t}M$ and $K_t:T^*_yM\toT_{\g_t}M$ by
\begin{equation}\label{JFE}
J_t\xi_0
=\left.\frac{\partial}{\partial{\varepsilon}}\right|_{{\varepsilon}=0}\pi\psi_t({\lambda}_0+{\varepsilon}\xi_0),\quad
K_t\xi_1
=\left.\frac{\partial}{\partial{\varepsilon}}\right|_{{\varepsilon}=0}\pi\psi_{-(1-t)}({\lambda}_1-{\varepsilon}\xi_1).
\end{equation}
We assume that $(x,y)$ is non-conjugate for $({\lambda}_t)_{t\in[0,1]}$, meaning that 
\begin{equation}\label{NC}
\text{$J_1$ is invertible.}
\end{equation}
Following Bismut \cite{B}, when conditions \eqref{UPM}, \eqref{PBC} and \eqref{NC} hold, the pair $(x,y)$ is said to lie outside the cut locus of $a$.
By a simple and well known argument, \eqref{NC} implies that the bicharacteristic projecting to ${\gamma}$ is unique.

The following statement is part of Theorem \ref{SRMR}.
\begin{theorem}\label{MRM}
Assume the hypotheses of Theorem \ref{MR}. Then the limit measure $\mu_\g$ is the unique zero-mean Gaussian measure on ${T_\g\oxy}$ with covariance given for $s\leqslant t$ by
\begin{equation}\label{ZMG}
\int_{T_\g\oxy} v_s\otimes v_t\,\mu_{\gamma}(dv)=J_sJ_1^{-1}K_t^*.
\end{equation}
\end{theorem}
The characterization of the cut locus and limit measure $\mu_\g$ in terms of the bicharacteristic flow is computationally effective and does not require the
construction of the quadratic form $Q$.
On the other hand, the formulation of Theorem \ref{MR} in terms of $Q$ confirms in the setting of an infinite-dimensional path space 
an intuition derived from analogous considerations for functions and measures in finite dimensions.
Moreover, the proof of Theorem \ref{MR} is by analysis in path space and leads naturally to the formulation given.

We conclude this section with some remarks on symmetry under time-reversal.
The following calculation shows that $J_1=K_0^*$ and hence that the cut locus is symmetric.
Since $V$ is Hamiltonian, its flow preserves the symplectic form ${\beta}$.
See for example \cite{M}.
For $\xi\inT^* M$, write $\tilde\xi$ for the corresponding vertical vector in $TT^* M$ and write $\psi_t^*$ for the action of $\psi_t$ on $TT^* M$. 
Then
$$
\<J_1\xi_0,\xi_1\>
=\<\pi^*\psi_1^*\tilde\xi_0,\xi_1\>
={\beta}(\psi_1^*\tilde\xi_0,\tilde\xi_1)
={\beta}(\tilde\xi_0,\psi_{-1}^*\tilde\xi_1)
=-\<\xi_0,\pi^*\psi_{-1}^*\tilde\xi_1\>
=\<\xi_0,K_0\xi_1\>.
$$
When ${\cal L}$ has the form \eqref{CFO}, for compactly supported $C^\infty $ functions $f,g$ on $M$, we have
$$
\int_Mf{\cal L} g\,d\nu=\int_Mg\hat{\cal L} f\,d\nu.
$$
where 
$$
\hat{\cal L} f=\tfrac12\operatorname{div}(a\nabla f)-a({\beta},\nabla f).
$$
and the associated bridge measure satisfies $\hat\mu_{\varepsilon}^{y,x}=\mu^{x,y}_\ve\circ{\wedge}^{-1}$,
where $\wedge:\Omega^{x,y}\to{\Omega}^{y,x}$ is the time-reversal map, given by $\hat{\omega}_t={\omega}_{1-t}$.
Hence, in this case, the whole set-up is invariant under time-reversal.

\section{Brownian bridge on a Riemannian manifold}\label{RIEM}
Before proving Theorem \ref{MR}, we discuss its consequences in the case where $a$ is everywhere positive-definite.
Let $M$ be a connected $C^\infty $ Riemannian manifold.
Write $\Delta$ for the Laplace--Beltrami operator on $M$ and set 
$$
{\cal L}=\tfrac12\Delta+\bar X_0
$$
where $\bar X_0$ is a $C^\infty $ vector field on $M$. 
Then ${\cal L}$ has the form \eqref{LAB} and its diffusivity $a$ is the inverse of the metric tensor.
We obtain in this way all operators ${\cal L}$ of the form \eqref{LAB} with $a$ everywhere positive-definite.
In this context, it is well known that there exist $m\in{\mathbb N}$ and $C^\infty$ vector fields $X_1,\dots,X_m$ on $M$ such that \eqref{LAC} and \eqref{LAD} hold, with
$$
\operatorname{span}\{Y(x):Y\in{\cal A}(X_1,\dots,X_m)\}=\operatorname{span}\{X_1(x),\dots,X_m(x)\}=T_xM\quad\text{for all $x\in M$.}
$$
We now review a number of equivalent characterizations of the cut locus for Riemannian manifolds, 
and we give some further descriptions of the limit Gaussian measure $\mu_{\gamma}$ using objects of Riemannian geometry.
In the case $\bar X_0=0$, the measure $\mu^{x,y}_\ve$ is the law of the Riemannian Brownian bridge from $x$ to $y$ of speed $\sqrt{\varepsilon}$.

Fix $x,y\in M$ and assume, as above, that there is a unique strongly minimal path ${\gamma}\inH^{x,y}$.
It is well known that in this context ${\gamma}$ is always the projection of a bicharacteristic.
Write $\nabla$ for the Levi--Civita covariant derivative and $R$ for the Riemann curvature tensor. 
Define a linear map $R_t:T_{\g_t}M\toT_{\g_t}M$, symmetric with respect to the metric, by $R_t=R(.,\dot{\gamma}_t)\dot{\gamma}_t$. 
It is straightforward to see that the set ${T_\g H^{x,y}}$, defined more generally in Section \ref{SUB},
is here given by the set of all absolutely continuous paths $v$ in ${T_\g\oxy}$ such that
$$
\int_0^1|\nabla v_t|^2dt<\infty.
$$
Moreover, we can define, for $\eta\in{\mathbb R}$ sufficiently small, a path ${\gamma}^\eta\inH^{x,y}$ by ${\gamma}^\eta_t=\operatorname{exp}_{{\gamma}_t}(\eta v_t)$. Then the map $\eta\mapsto I({\gamma}^\eta)$ is twice differentiable near $0$ with
\begin{equation}\label{QR}
Q(v)=\left.\frac{\partial^2}{\partial\eta^2}\right|_{\eta=0}I({\gamma}^\eta)=\int_0^1|\nabla v_t|^2\, dt-\int_0^1\<v_t,R_tv_t\>\, dt.
\end{equation}
See, for example \cite{K}.
By a standard calculation, the processes $(J_t)_{t\in[0,1]}$ and $(K_t)_{t\in[0,1]}$, defined above, are Jacobi fields along ${\gamma}$ and satisfy the differential equations
$$
\nabla^2J_t+R_tJ_t=0,\quad J_0=0,\quad\nabla J_0=a(x)
$$
and
$$
\nabla^2K_t+R_tK_t=0,\quad K_1=0,\quad\nabla K_1=-a(y).
$$
Let $(b_t)_{t\in[0,1]}$ be a Brownian motion in $T_x M$, starting from $0$.
Set $z_t=b_t-tb_1$. 
Then $(z_t)_{t\in[0,1]}$ is a Brownian bridge in $T_x M$ from $0$ to $0$ in time $1$. 
Let $({\tau}_t)_{t\in[0,1]}$ denote parallel translation along ${\gamma}$, thus ${\tau}_t\inT_{\g_t}M\otimesT^*_x M$ and ${\tau}_0=\operatorname{id}$, $\nabla{\tau}_t=0$.
Let $\bar\mu$ denote the law of $({\tau}_tz_t)_{t\in[0,1]}$ on ${T_\g\oxy}$.
The following result is well known.
See \cite{K} and, for (vi), \cite[Theorem 4.17]{B}.
Although it is framed geometrically, it is in fact only a result about quadratic equations in matrices and Gaussian processes,
as one can see by choosing a chart along ${\gamma}$ such that $a({\gamma}_t)={\tau}_t=\operatorname{id}$ for all $t$. 
\begin{proposition}\label{1.5}
The following are equivalent
\begin{itemize}
\item[{\rm (i)}] $J_t$ is invertible for all $t\in(0,1]$,
\item[{\rm (ii)}] $K_t$ is invertible for all $t\in[0,1)$,
\item[{\rm (iii)}] there exists a $C^1$ path $(A_t)_{t\in[0,1)}$ along ${\gamma}$, with $A_t\in T_{{\gamma}_t}M\otimes T^*_{{\gamma}_t}M$,
solving the Riccati equation
$$
\nabla A_t+A_t^2+R_t=0,\quad (1-t)A_t\rightarrow-\operatorname{id}\quad {\text as}\quad t\to 1,
$$
\item[{\rm (iv)}] there exists a $C^1$ path $(B_t)_{t\in(0,1]}$ along ${\gamma}$, with $B_t\in T_{{\gamma}_t}M\otimes T^*_{{\gamma}_t}M$,
solving the Riccati equation
$$
\nabla B_t+B_t^2+R_t=0,\quad tB_t\rightarrow\operatorname{id}\quad {\text as}\quad t\to 0,
$$
\item[{\rm (v)}] $Q$ is positive-definite on ${T_\g H^{x,y}}$,
\item[{\rm (vi)}] we have
$$
\int_{T_\g\oxy}\operatorname{exp}\left\{\frac12\int_0^1\<v_t,R_tv_t\>\, dt\right\}\bar\mu(dv)<\infty.
$$
\end{itemize}
Moreover, under these conditions, we have 
$$
J_1=K_0^*,\quad \nabla K_t=A_tK_t,\quad \nabla J_t=B_tJ_t
$$
and for $t<1$
$$
J_t=K_t\int_0^tK_s^{-1}a({\gamma}_s)(K_s^{-1})^*\,ds\,K_0^*.
$$
\end{proposition}
It is standard that the condition that ${\gamma}$ is minimal implies already that $Q$ is non-negative and that $J_t$ is invertible for all $t\in(0,1)$.
In this context then, condition (i) is equivalent to the condition that $x$ and $y$ are non-conjugate along ${\gamma}$, that is,
there is no non-trivial vector field $(v_t)_{t\in[0,1]}$ along ${\gamma}$ vanishing at the endpoints and such that $\nabla^2 v_t+R_tv_t=0$ for all $t$.
In Section \ref{SUB}, we will generalize the equivalence of (i), (ii) and (v) to the sub-Riemannian case.

The following result gives three further characterizations for the Gaussian measure $\mu_{\gamma}$,
which by Theorem \ref{MR} describes the small-time fluctuations of the Brownian bridge in $M$.
As in the preceding result, by choice of a suitable chart along ${\gamma}$, we can reduce to the case where $M={\mathbb R}^d$ with $a({\gamma}_t)={\tau}_t=\operatorname{id}$ for all $t$. 
The result is then of a standard type for Gaussian processes.
See for example \cite{NV}.
\begin{theorem}\label{RMR}
Let $M$ be a connected Riemannian manifold and let $x,y\in M$.
Suppose that there is a unique minimal path ${\gamma}\inH^{x,y}$ and that $(x,y)$ is non-conjugate along ${\gamma}$.
Then there exists a zero-mean Gaussian probability measure $\mu_\g$ on ${T_\g\oxy}$
with the following properties:
\begin{itemize}
\item[{\rm (i)}] for $s\leqslant t$, we have
$$
\int_{{T_\g\oxy}} v_s\otimes v_t\,\mu_\g(dv)=J_sJ_1^{-1}K_t^*,
$$
\item[{\rm (ii)}] for all continuous linear functionals $\phi$ on ${T_\g\oxy}$, we have
$$
\int_{T_\g\oxy}\phi(v)^2\mu_\g(dv)=Q(\tilde\phi)
$$
where $\tilde\phi\in{T_\g H^{x,y}}$ is determined by $\phi(v)=Q(\tilde\phi,v)$ for all $v\in{T_\g H^{x,y}}$,
\item[{\rm (iii)}] under $\mu_\g$, the coordinate process $v$ on ${T_\g\oxy}$ satisfies a covariant linear stochastic differential equation over ${\gamma}$ of the form
$$
\nabla v_t={\tau}_tdb_t+A_tv_tdt,\quad v_0=0,
$$
\item[{\rm (iv)}] $\mu_\g$ is absolutely continuous with respect to $\bar\mu$, with Radon--Nikodym derivative 
\begin{equation}\label{exp}
\frac{d\mu_\g}{d\bar\mu}(v)\propto\operatorname{exp}\left\{\frac12\int_0^1\<v_t,R_tv_t\>\, dt\right\}.
\end{equation}
\end{itemize}
Moreover, any one of these properties characterizes $\mu_\g$ uniquely.
\end{theorem}
The equivalence of (i) and (ii) will be established in a more general context in Section \ref{QG}.
Here is an illustrative calculation, taking advantage of the reduction to ${\tau}_t=\operatorname{id}$ mentioned above. 
Suppose that $\mu$ satisfies (iii). 
Set $w_t=K_t^{-1}v_t$, then $dw_t=K_t^{-1}db_t$, so
$$
w_t=\int_0^tK_s^{-1}db_s.
$$
Hence, for $0\leqslant s\leqslant t\leqslant 1$,
$$
\int_{{T_\g\oxy}} v_s\otimes v_t\,\,\mu(dv)=K_s\left(\int_0^sK_r^{-1}(K_r^{-1})^*dr\right)K_t^*=J_sJ_1^{-1}K_t^*.
$$
Hence $\mu$ satisfies (i).

We examine now how our result specializes in some simple cases. 
When $M={\mathbb R}^d$ with Euclidean metric, the analysis is trivial, because $\tilde\mu^{x,y}_\ve=\mu_\g$ for all ${\varepsilon}>0$. 
Then
$$
J_s=sI,\quad K_t=(1-t)I,\quad R_t=0,\quad A_t=-I/(1-t)
$$
so the alternatives in Theorem \ref{RMR} are simply some of the standard
descriptions of the Brownian bridge in ${\mathbb R}^d$.

In the case where $M$ is a sphere or hyperbolic space, we can rewrite \eqref{exp} in the form
$$
\frac{d\mu_\g}{d\bar\mu}(y)
\propto\operatorname{exp}\left\{\frac{ K d(x,y)^2}{2}\int_0^1|y_t|^2\, dt\right\},
$$
where $K$, the scalar curvature, is $1$ for the sphere and $-1$ for hyperbolic space. 
Thus, on a sphere, the variance of the fluctuations is larger than in ${\mathbb R}^d$, whereas, in hyperbolic space it is less.
This does not contradict the tendency of Brownian paths to separate quickly in hyperbolic space because we are conditioning on the endpoint.
Thus we tend to see those paths which have never deviated far from the geodesic.

\section{Laplace's method on Wiener space}\label{LAP}
In this section we prove Theorem \ref{MR} in the case where conditions \eqref{HSS} and \eqref{URD} hold, following closely the method used by Ben Arous \cite{BA} to study the heat kernel in this case.
It will suffice to consider the case where ${\theta}(z)=z$ for all $z\in{\mathbb R}^d$.
The vector fields $X_0,X_1,\dots,X_m$ provide both a means to construct a diffusion process $(x^{\varepsilon}_t)_{t\ge0}$ with generator ${\varepsilon}{\cal L}$ starting from $x$ 
and a parametrization of the set $H^x$ of finite energy paths starting from $x$. 
Let $(B_t)_{t\in[0,1]}$ be a Brownian motion in ${\mathbb R}^m$, which we assume is realized as the coordinate process on ${\Omega}^0({\mathbb R}^m)=\{w\in C([0,1],{\mathbb R}^m):w_0=0\}$ under Wiener measure $\mathbb{P}$.
Define a vector field $\tilde X_0$ on ${\mathbb R}^d$ by
$$
\tilde X_0^i(x)=X_0^i(x)+\frac12\sum_{\ell=1}^m\<\nabla X^i_\ell(x),X_\ell(x)\>.
$$
Consider the It\^o stochastic differential equation in ${\mathbb R}^d$
$$
dx_t^{\varepsilon}=\sqrt{\varepsilon}\sum_{\ell=1}^mX_\ell(x_t^{\varepsilon})dB_t^\ell+{\varepsilon}\tilde X_0(x_t^{\varepsilon})dt,\quad x_0^{\varepsilon}=x.
$$
This has a unique strong solution $(x_t^{\varepsilon})_{t\in[0,1]}$, whose law on ${\Omega}=C([0,1],{\mathbb R}^d)$ is $\mu_{\varepsilon}^x$.
There is no explosion.
Write $H^0({\mathbb R}^m)$ for the set of Cameron--Martin paths $h$ in ${\mathbb R}^m$, that is to say, the set of absolutely continuous functions $h:[0,1]\to{\mathbb R}^m$, starting from $0$, such that
$$
\int_0^1|\dot h_t|^2dt<\infty.
$$
Given $h\in H^0({\mathbb R}^m)$, consider the differential equation in ${\mathbb R}^d$
\begin{equation}\label{FEP}
d\phi_t=\sum_{\ell=1}^mX_\ell(\phi_t)dh^\ell_t,\quad \phi_0=x.
\end{equation}
There is a unique solution $\phi(x,h)=(\phi_t(x,h))_{t\in[0,1]}$ in ${\Omega}$.
As we will show in Proposition \ref{2.1}, in fact $\phi(x,h)\in H^x$ and $\phi(x,.)$ maps $H^0({\mathbb R}^m)$ onto $H^x$.
Recall that ${\gamma}$ is the minimizing path in $H^{x,y}$ and ${\gamma}$ is the projection of a bicharacteristic ${\lambda}$.
Define $h\in H^0({\mathbb R}^m)$ by $\dot h_t^\ell=\<{\lambda}_t,X_\ell({\gamma}_t)\>$.
Then ${\gamma}=\phi(x,h)$ and $I({\gamma})=\int_0^1|\dot h_t|^2dt$.
We reserve the notation $h$ for this minimizing control path from now on.
For $s\in[0,1]$, we will write $(\phi_{ts}(x,h))_{t\in[s,1]}$ for the solution to \eqref{FEP} starting from $x$ at time $s$.
We denote the derivative in $x$ by $\phi^*_{ts}(x,h)$ and set $u_t=\phi_t^*(x,h)$.
Then $u_t$ is invertible for all $t$ and 
$$
du_t=\sum_{\ell=1}^m\nabla X_\ell({\gamma}_t)u_tdh^\ell_t,\quad u_0=I.
$$
Moreover, we have $\phi^*_{ts}({\gamma}_s,h)=u_tu_s^{-1}$.
We can define a continuous linear map $v:{\Omega}^0({\mathbb R}^m)\to{T_\g\oxy}$ by\footnote{%
Since $M={\mathbb R}^d$ in the present discussion, ${T_\g\oxy}$ can be naturally identified with a subset of ${\Omega}$, but we will keep the distinction anyway.}
\begin{equation}\label{VEQ}
v_t(w)=\sum_{\ell=1}^m\int_0^t\phi^*_{ts}({\gamma}_s,h)X_\ell({\gamma}_s)dw_s^\ell=\sum_{\ell=1}^mu_t\int_0^tu_s^{-1}X_\ell({\gamma}_s)dw_s^\ell
\end{equation}
where the integral is understood by a formal integration by parts.
Set $Y_t=v_t(B)$ and note that $Y_1$ is a zero-mean Gaussian random variable in ${\mathbb R}^d$ having covariance matrix $\bar C_1=u_1C_1u_1^*$, where
$$
C_1=\sum_{\ell=1}^m\int_0^1(u_t^{-1}X_\ell({\gamma}_t))\otimes(u_t^{-1}X_\ell({\gamma}_t))dt.
$$
In \cite{B}, Bismut called $C_1$ the deterministic Malliavin covariance matrix.
Condition \eqref{NC} implies that $C_1$ is invertible. 
This follows in particular from \eqref{al}.
Hence $Y_1$ has a density function with respect to Lebesgue measure on ${\mathbb R}^d$, given by
$$
\bar p(z)=(2\pi)^{-d/2}(\det\bar C_1)^{-1/2}\operatorname{exp}\{-\<z,\bar C_1^{-1}z\>/2\}.
$$
Set\footnote{%
Formally, we have
$$
Y_t=\frac\partial{\partial h}\phi_t(x,h)(B),\quad Y^{(2)}_1=\frac{\partial^2}{\partial h^2}\phi_1(x,h)(B,B)
$$
which may be seen by differentiating \eqref{FEP} twice in $h$ in the direction $B$ and solving the resulting equations by variation of constants.
}
$$
Y^{(2)}_1=\sum_{\ell=1}^m\int_0^1\phi^*_{1t}({\gamma}_t,h)\left\{\nabla^2X_\ell({\gamma}_t)(Y_t,Y_t)dh_t^\ell+2\nabla X_\ell({\gamma}_t)Y_tdB_t^\ell\right\},\quad S=\<{\lambda}_1,Y_1^{(2)}\>.
$$
where the integral $Y_tdB_t$ is understood in the sense of It\^o.
Define a linear map ${\tau}:{\mathbb R}^d\toH^0(\R^m)$ by
$$
\dot{\tau}_t^\ell(z)=\left\<C_1^{-1}u_1^{-1}z,u_t^{-1}X_\ell({\gamma}_t)\right\>.
$$
Set $K=\{k\inH^0(\R^m):v_1(k)=0\}$.
From \eqref{VEQ} we see that ${\tau}$ maps ${\mathbb R}^d$ onto the orthogonal complement $K^\perp$ of $K$ in $H^0(\R^m)$. 
Moreover, $v_1({\tau}(z))=z$ and the restriction of the map ${\tau}\circ v_1$ to $H^0(\R^m)$ is the orthogonal projection $H^0(\R^m)\to K^\perp$.
Set
\begin{equation}\label{BWW}
W_t=B_t-W'_t,\quad W'_t={\tau}_t(v_1(B)).
\end{equation}
Then $(W_t)_{t\in[0,1]}$ and $(W'_t)_{t\in[0,1]}$ are independent continuous Gaussian processes, and $v_1(W)=0$ and $W'\in K^\perp$ almost surely.
For $z\in{\mathbb R}^d$, set $W_t(z)=W_t+{\tau}_t(z)$ and $Y_t(z)=v_t(W(z))$ and $S(z)=\<{\lambda}_1,Y_1^{(2)}(z)\>$, where
\begin{equation}\label{SDEF}
Y^{(2)}_1(z)=\sum_{\ell=1}^m\int_0^1\phi^*_{1t}({\gamma}_t,h)\left\{\nabla^2X_\ell({\gamma}_t)(Y_t(z),Y_t(z))dh_t^\ell+2\nabla X_\ell({\gamma}_t)Y_t(z)dW_t^\ell(z)\right\}.
\end{equation}
We interpret the integral with respect to $W(z)$ by writing 
$$
Y_t(z)dW_t(z)=v_t(B-W'+{\tau}(z))d(B_t-W_t'+{\tau}_t(z))
$$ 
and expanding. 
The term $v_t(B)dB_t$ is considered as an It\^o integral, 
while the remaining terms can be considered as integrals with respect to Lebesgue measure, sometimes after a formal integration by parts.
It is straightforward to see, using the same version of the It\^o integral for all $z$, 
that the family of random variables $(Y^{(2)}_1(z):z\in{\mathbb R}^d)$ is continuous in $z$.
Also, $Y(z)$ and $Y_1^{(2)}(z)$ are independent of $v_1(B)$ for all $z$.
Note that $W(v_1(B))=B$ and $Y(v_1(B))=Y$, and that $S(v_1(B))=S$ almost surely. 
Note also that $Y^{(2)}_1(z)$ belongs to the sum of the zeroth, first and second Wiener chaoses in $L^2({\Omega}^0({\mathbb R}^m),\mathbb{P})$ for all $z$.
By Theorem \ref{SRMR}, we have ${\mathbb E}(e^{pS(z)/2})<\infty$ for all $z\in{\mathbb R}^d$ for some $p>1$.

Consider the function 
$$
E(z)=d(x,z)^2/2=\inf_{{\omega}\in H^{x,z}}I({\omega})/2.
$$ 
As Bismut \cite[Theorem 1.26]{B} showed, $E$ is $C^\infty $ in a neighbourhood $N$ of $y$, with $E'(y)={\lambda}_1$.
Following Ben Arous \cite[Lemma 3.8]{BA}, there then exists a function $F\in C^\infty_b({\mathbb R}^d)$ such that the map $F+E$ has a unique and non-degenerate 
minimum at $y$, with minimum value $0$.
To see this, choose a neighbourhood $N_0$ of $y$ and a $C^\infty$ function $\chi$ of compact support such that $1_{N_0}\leqslant\chi\le1_N$.
Fix a constant ${\alpha}>0$ and consider the function
$$
F(z)=\chi(z)\left({\alpha}|y-z|^2-E(y)-E'(y)(y-z)\right)+(1-\chi(z)).
$$
Then $F\in C^\infty_b({\mathbb R}^d)$ and $F(y)+E(y)=0$ and $F'(y)+E'(y)=0$.
Moreover, by choosing ${\alpha}$ sufficiently large, we can ensure that $F''(y)+E''(y)$ is positive-definite and $F(z)+E(z)>0$ for all $z\not=y$,
so $F+E$ has a non-degenerate minimum at $y$, which is also its global minimum.
We fix a choice of a function $F$ with the given properties for the rest of the analysis.

Set ${\gamma}^0={\gamma}$ and define $({\gamma}_t^{\varepsilon})_{t\in[0,1]}$ for ${\varepsilon}>0$ as the strong solution of the stochastic differential equation
\begin{equation}\label{GSDE}
d{\gamma}_t^{\varepsilon}=\sum_{\ell=1}^mX_\ell({\gamma}_t^{\varepsilon})dh_t^\ell+\sqrt{\varepsilon}\sum_{\ell=1}^mX_\ell({\gamma}_t^{\varepsilon})dB_t^\ell+{\varepsilon}\tilde X_0({\gamma}_t^{\varepsilon})dt,\quad {\gamma}_0^{\varepsilon}=x.
\end{equation}
By standard results on stochastic differential equations, for all $t\in[0,1]$ and all $p\in[1,\infty)$, the map ${\varepsilon}\mapsto{\gamma}_t^{\varepsilon}:[0,\infty)\to L^p(\mathbb{P})$ is continuous.
Furthermore, we can and do choose versions so that, almost surely, the map ${\sigma}\mapsto{\gamma}_t^{{\sigma}^2}:[0,\infty)\to{\mathbb R}^d$ is $C^\infty $.
Moreover, the first and second derivatives at ${\sigma}=0$ then satisfy
$$
\left.\frac\partial{\partial{\sigma}}\right|_{{\sigma}=0}{\gamma}_t^{{\sigma}^2}=Y_t,\quad 
\left.\left(\frac\partial{\partial{\sigma}}\right)^2\right|_{{\sigma}=0}{\gamma}_1^{{\sigma}^2}=Y^{(2)}_1+Z_1,\quad
Z_1=\int_0^1\phi_{1t}^*({\gamma}_t,h)\tilde X_0({\gamma}_t)dt.
$$
Now the map $f({\sigma})=F({\gamma}_1^{{\sigma}^2})$ is $C^\infty $ on $[0,\infty)$ and $F(y)=-d(x,y)^2/2$ and $F'(y)=-{\lambda}_1$, so
$$
f(0)=F(y)=-\frac12\int_0^1|\dot h_t|^2dt,\quad
f'(0)=F'(y)Y_1=-\sum_{\ell=1}^m\int_0^1\dot h_t^\ell dB_t^\ell
$$
and
$$
f''(0)=F'(y)(Y^{(2)}_1+Z_1)+F''(y)(Y_1,Y_1).
$$
Set
$$
R({\varepsilon})=\int_0^1(1-{\theta})f''({\theta}\sqrt{\varepsilon})d{\theta}.
$$
Then, by Taylor's theorem,
$$
F({\gamma}_1^{\varepsilon})=f(\sqrt{\varepsilon})=f(0)+\sqrt{\varepsilon} f'(0)+{\varepsilon} R({\varepsilon})=-\frac12\int_0^1|\dot h_t|^2dt-\sqrt{\varepsilon}\sum_{\ell=1}^m\int_0^1\dot h_t^\ell dB_t^\ell+{\varepsilon} R({\varepsilon}).
$$
Set 
$$
\tilde x^{\varepsilon}_t=\frac{x_t^{\varepsilon}-{\gamma}_t}{\sqrt{\varepsilon}},\quad\tilde{\gamma}_t^{\varepsilon}=\frac{{\gamma}^{\varepsilon}_t-{\gamma}_t}{\sqrt{\varepsilon}}
$$
and note that, for all $t\in[0,1]$, we have $\tilde{\gamma}^{\varepsilon}_t\to Y_t$ as ${\varepsilon}\to0$ almost surely. 
For ${\varepsilon}>0$, consider the new probability measure $\mathbb{P}^{\varepsilon}$ on ${\O^0(\R^m)}$ given by $d\mathbb{P}^{\varepsilon}/d\mathbb{P}=\rho_1^{\varepsilon}$, where
$$
\rho_1^{\varepsilon}=\operatorname{exp}\left\{-\frac1{\sqrt{\varepsilon}}\sum_{\ell=1}^m\int_0^1\dot h_t^\ell dB_t^\ell-\frac1{2{\varepsilon}}\int_0^1|\dot h_t|^2dt\right\}=\operatorname{exp}\left\{\frac{F({\gamma}_1^{\varepsilon})}{\varepsilon}-R({\varepsilon})\right\}.
$$
Note that \eqref{GSDE} can be rewritten in the form
$$
d{\gamma}_t^{\varepsilon}=\sqrt{\varepsilon}\sum_{\ell=1}^mX_\ell({\gamma}_t^{\varepsilon})dB_t^{{\varepsilon},\ell}+{\varepsilon}\tilde X_0({\gamma}_t^{\varepsilon})dt,\quad {\gamma}_0^{\varepsilon}=x
$$
where $B^{{\varepsilon},\ell}_t=B^\ell_t+h_t^\ell/\sqrt{\varepsilon}$.
By the Cameron-Martin formula, under $\mathbb{P}^{\varepsilon}$, the process $(B^{\varepsilon}_t)_{t\in[0,1]}$ is a Brownian motion, so ${\gamma}^{\varepsilon}$ has law $\mu_{\varepsilon}^x$.

At this point we modify the argument of Ben Arous by introducing a smooth cylindrical function $G$ on ${\Omega}$ which serves to keep track of the paths of the diffusion bridge.
Fix $t_1,\dots,t_k\in(0,1)$ with $t_1<\dots<t_k$ and a $C^\infty $ function $g$ on $({\mathbb R}^d)^k$ of polynomial growth.
Set $G({\omega})=g({\omega}_{t_1},\dots,{\omega}_{t_k})$.
Define for $z\in{\mathbb R}^d$ and ${\varepsilon}>0$
\begin{align*}
G_{\varepsilon}(z)
&={\varepsilon}^{d/2}p({\varepsilon},x,y+\sqrt{\varepsilon} z)e^{-F(y+\sqrt{\varepsilon} z)/{\varepsilon}}\int_{T_\g\oxy} G({\omega})\tilde\mu_{\varepsilon}^{x,y+\sqrt{\varepsilon} z}(d{\omega}),\\
G_0(z)
&=\bar p(z)e^{\<{\lambda}_1,Z_1\>/2-F''(y)(z,z)/2}{\mathbb E}(G(Y(z))e^{S(z)/2}).
\end{align*}
Then $G_{\varepsilon}$ and $G_0$ are continuous integrable functions on ${\mathbb R}^d$.
Consider the Fourier transform
$$
\hat G_{\varepsilon}(\xi)=\int_{{\mathbb R}^d}G_{\varepsilon}(z)e^{i\<\xi,z\>}dz.
$$
Then
\begin{align}\notag
\hat G_{\varepsilon}(\xi)
&=\int_{{\mathbb R}^d}\int_{\Omega} p({\varepsilon},x,y')e^{-F(y')/{\varepsilon}}G\left(\frac{{\omega}-{\gamma}}{\sqrt{\varepsilon}}\right)\mu_{\varepsilon}^{x,y'}(d{\omega})e^{i\<\xi,(y'-y)/\sqrt{\varepsilon}\>}dy'\\\notag
&=\int_{\Omega} e^{-F({\omega}_1)/{\varepsilon}}G\left(\frac{{\omega}-{\gamma}}{\sqrt{\varepsilon}}\right)e^{i\<\xi,({\omega}_1-{\gamma}_1)/\sqrt{\varepsilon}\>}\mu_{\varepsilon}^x(d{\omega})\\\notag
&={\mathbb E}\left(G(\tilde x^{\varepsilon})\operatorname{exp}\{i\<\xi,\tilde x_1^{\varepsilon}\>-F(x_1^{\varepsilon})/{\varepsilon}\}\right)\\\notag
&={\mathbb E}\left(G(\tilde{\gamma}^{\varepsilon})\operatorname{exp}\{i\<\xi,\tilde {\gamma}_1^{\varepsilon}\>-F({\gamma}_1^{\varepsilon})/{\varepsilon}\}\rho_1^{\varepsilon}\right)\\\label{GEP}
&={\mathbb E}\left(G(\tilde{\gamma}^{\varepsilon})\operatorname{exp}\{i\<\xi,\tilde {\gamma}_1^{\varepsilon}\>-R({\varepsilon})\}\right)
\end{align}
and
\begin{align*}
\hat G_0(\xi)
&=\int_{{\mathbb R}^d}\bar p(z)e^{\<{\lambda}_1,Z_1\>/2-F''(y)(z,z)/2}{\mathbb E}(G(Y(z))e^{S(z)/2})e^{i\<\xi,z\>}dz\\
&={\mathbb E}\left(G(Y)\operatorname{exp}\{i\<\xi,Y_1\>+(S+\<{\lambda}_1,Z_1\>-F''(y)(Y_1,Y_1))/2\}\right)\\
&={\mathbb E}\left(G(Y)\operatorname{exp}\{i\<\xi,Y_1\>-f''(0)/2\}\right).
\end{align*}
Consider the limit ${\varepsilon}\to0$.
We have $\tilde{\gamma}^{\varepsilon}_t\to Y_t$ almost surely and in $L^p$ for all $p<\infty$, for all $t\in[0,1]$.
Now $R({\varepsilon})\to f''(0)/2$ almost surely as ${\varepsilon}\to0$ and we have the following key estimate proved by Ben Arous \cite[Lemma 3.25]{BA}: there exists $p>1$ such that
\begin{equation}\label{BKE}
\limsup_{{\varepsilon}\to0}{\mathbb E}(e^{-pR({\varepsilon})})<\infty.
\end{equation}
Hence $\hat G_{\varepsilon}(\xi)\to\hat G_0(\xi)$ for all $\xi\in{\mathbb R}^d$.
We will prove the following key lemma at the end of this section.
\begin{lemma}\label{KL}
There exists a constant $C(G)<\infty$ such that, for all ${\varepsilon}\in(0,1]$ and all $\xi\in{\mathbb R}^d$, we have
\begin{equation}\label{HRR}
|\hat G_{\varepsilon}(\xi)|\leqslant C(G)/(1+|\xi|^{d+1}).
\end{equation}
Moreover, there exists $C<\infty$ such that, uniformly in $s,t\in[0,1]$, in the case where $G({\omega})=|{\omega}_s-{\omega}_t|^4$, 
for all ${\varepsilon}\in(0,1]$ and all $\xi\in{\mathbb R}^d$, we have
\begin{equation}\label{GST}
|\hat G_{\varepsilon}(\xi)|\leqslant C|s-t|^2/(1+|\xi|^{d+1}).
\end{equation}
\end{lemma}
The lemma allows us to use dominated convergence in the Fourier inversion formula to deduce that
$$
G_{\varepsilon}(0)=(2\pi)^{-d}\int_{{\mathbb R}^d}\hat G_{\varepsilon}(\xi)d\xi\to(2\pi)^{-d}\int_{{\mathbb R}^d}\hat G_0(\xi)d\xi=G_0(0).
$$
that is to say
\begin{equation}\label{KF}
{\varepsilon}^{d/2}p({\varepsilon},x,y)e^{d(x,y)^2/(2{\varepsilon})}\int_{T_\g\oxy} G({\omega})\tilde\mu_{\varepsilon}^{x,y}(d{\omega})\to(2\pi)^{-d/2}(\det\bar C_1)^{-1/2}e^{\<{\lambda}_1,Z_1\>/2}{\mathbb E}(G(Y(0))e^{S(0)/2})
\end{equation}
where we used the fact that $\bar p(0)=(2\pi)^{-d/2}(\det\bar C_1)^{-1/2}$.
On taking $g=1$, we recover the heat kernel asymptotics shown by Ben Arous
\begin{equation}\label{BBF}
{\varepsilon}^{d/2}p({\varepsilon},x,y)e^{d(x,y)^2/(2{\varepsilon})}\to(2\pi)^{-d/2}(\det\bar C_1)^{-1/2}e^{\<{\lambda}_1,Z_1\>/2}{\mathbb E}(e^{S(0)/2}).
\end{equation}
By Theorem \ref{SRMR}, $\mu_\g$ is the law of $Y(0)$ on ${T_\g\oxy}$ under the probability measure $\tilde\mathbb{P}$, where $d\tilde\mathbb{P}/d\mathbb{P}\propto e^{S(0)/2}$. 
Hence, on dividing \eqref{KF} by \eqref{BBF}, we obtain
$$
\int_{T_\g\oxy} G({\omega})\tilde\mu^{x,y}_\ve(d{\omega})\to\int_{T_\g\oxy} G({\omega})\mu_\g(d{\omega})
$$
from which we deduce that the finite-dimensional distributions of $\tilde\mu^{x,y}_\ve$ converge weakly to those of $\mu_{\gamma}$.
Finally, on taking $G({\omega})=|{\omega}_s-{\omega}_t|^4$ and using the estimate \eqref{GST}, we find by Fourier inversion that 
$$
{\varepsilon}^{d/2}p({\varepsilon},x,y)e^{d(x,y)^2/(2{\varepsilon})}\int_{T_\g\oxy}|{\omega}_s-{\omega}_t|^4\tilde\mu^{x,y}_\ve(d{\omega})=G_{\varepsilon}(0)\leqslant C|s-t|^2.
$$
Since the right-hand side in \eqref{BBF} is positive, we deduce that for some $C<\infty$ we have, for all $s,t\in[0,1]$,
$$
\sup_{{\varepsilon}\in(0,1]}\int_{T_\g\oxy}|{\omega}_s-{\omega}_t|^4\tilde\mu^{x,y}_\ve(d{\omega})\leqslant C|s-t|^2.
$$
Hence, by standard arguments, the family of laws $(\tilde\mu^{x,y}_\ve:{\varepsilon}\in(0,1])$ is tight on ${T_\g\oxy}$ and so $\tilde\mu^{x,y}_\ve\to\mu_\g$ weakly as ${\varepsilon}\to0$.
In particular, we also have $\mu^{x,y}_\ve\to\d_\g$, so the proof is complete.

\begin{proof}[Proof of Lemma \ref{KL}]
The idea is to integrate by parts in \eqref{GEP}, $d+1$ times, using Malliavin calculus. 
We will vary the argument of \cite{BA} in three respects.
First, we will present the argument using Bismut's integration by parts formula \cite{MR621660,MR942019} for solutions of stochastic differential equations. 
This is more direct and may be followed in detail without knowledge of the general apparatus of Malliavin calculus. 
Second, we will use the corrected argument \cite{MES} for the uniform non-degeneracy of the Malliavin covariance matrix.
Third, we will include the simple modifications needed to go beyond the case $G=1$ which is covered in \cite{BA}.
Within the proof, we will use a few notations which conflict with usage elsewhere.

Fix ${\sigma}=\sqrt{\varepsilon}>0$.
Define processes $(x_t)_{t\in[0,1]}$ in ${\mathbb R}^d$ and $(u_t)_{t\in[0,1]}$, $(v_t)_{t\in[0,1]}$ in ${\mathbb R}^d\otimes({\mathbb R}^d)^*$ as the strong
solutions of the following system of stochastic differential equations
\begin{align}\label{DGE}
dx_t
&=\sum_{\ell=1}^mX_\ell(x_t)dh^\ell_t+\sum_{\ell=1}^m{\sigma} X_\ell(x_t)dB^\ell_t+{\sigma}^2\tilde X_0(x_t)dt, \quad x_0=x\\\notag
du_t
&=\sum_{\ell=1}^m\nabla X_\ell(x_t)u_tdh^\ell_t+\sum_{\ell=1}^m{\sigma}\nabla X_\ell(x_t)u_tdB^\ell_t
+{\sigma}^2\nabla\tilde X_0(x_t)u_tdt, \quad u_0=I\\
\label{VGE}
dv_t
&=-\sum_{\ell=1}^mv_t\nabla X_\ell(x_t)dh^\ell_t-\sum_{\ell=1}^m{\sigma} v_t\nabla X_\ell(x_t)dB^\ell_t
-{\sigma}^2v_t\left\{\nabla\tilde X_0-\sum_{\ell=1}^m(\nabla X_\ell)^2\right\}(x_t)dt, \quad v_0=I.
\end{align}
Then $x_t={\gamma}^{{\sigma}^2}_t$ and, by It\^o's formula, $v_t=u_t^{-1}$ for all $t$.
It is well known that the random variables $\sup_{t\in[0,1]}|x_t|$, $\sup_{t\in[0,1]}|u_t|$ and $\sup_{t\in[0,1]}|v_t|$ have moments of all orders, 
which are bounded uniformly in ${\sigma}\in[0,1]$.

Consider the following stochastic differential equation in ${\mathbb R}^N$ with $C^\infty $ coefficients
\begin{equation}\label{ZT}
dz_t=\sum_{\ell=1}^mW_\ell(z_t)dh_t^\ell+\sum_{\ell=1}^mZ_\ell(z_t)dB^\ell_t+Z_0(z_t)dt,\quad z_0=z.
\end{equation}
We assume that the coefficients have a graded Lipschitz structure.
By this we mean that the coefficients and all their derivatives satisfy polynomial growth bounds, 
and that there exist $k\in{\mathbb N}$ and $N_1,\dots,N_k\in{\mathbb N}$ and a decomposition $z_t=(z_t^1,\dots,z_t^k)$ such that $N_1+\dots+N_k=N$,
and, for $j=1,\dots,k$, the component $z_t^j$ takes values in ${\mathbb R}^{N_j}$ and, for all $\ell$, 
the corresponding components $W_\ell^j$ and $Z_\ell^j$ of the coefficients depend only on $(z^1,\dots,z^j)$,
with $\partial W_\ell^j/\partial z^j$ and $\partial Z_\ell^j/\partial z^j$ uniformly bounded.
We will allow the coefficients in \eqref{ZT} to depend measurably on ${\sigma}$ and $t$, without making this explicit in the notation, 
while assuming that the bounds in their graded Lipschitz structure hold uniformly in ${\sigma}\in[0,1]$ and $t\in[0,1]$.
This is sufficient to ensure the existence and uniqueness of a strong solution $(z_t)_{t\in[0,1]}$ to \eqref{ZT} 
such that $\sup_{t\in[0,1]}|z_t|$ has moments of all orders, bounded uniformly in ${\sigma}\in[0,1]$.
See \cite{MR942019} for more details.

Fix $i\in\{1,\dots,d\}$ and consider for $\eta\in{\mathbb R}$ a perturbed process $(B_t^\eta)_{t\in[0,1]}$ in ${\mathbb R}^m$ given by
$$
dB_t^{\eta,\ell}=dB_t^\ell+\eta(v_tX_\ell(x_t))^idt,\quad B_0^\eta=0.
$$
Write $(z_t^\eta)_{t\in[0,1]}$ for the strong solution of the stochastic differential equation which is obtained when we replace $(B_t)_{t\in[0,1]}$ 
in \eqref{ZT} by this perturbed process.
We can and do choose a version of the family of processes $(z_t^\eta)_{t\in[0,1]}$ which is almost surely $C^\infty $ in $\eta$.
We associate to $(z_t)_{t\in[0,1]}$ the derived process $(z_t')_{t\in[0,1]}=((z_t')^1,\dots,(z_t')^d)_{t\in[0,1]}$ in ${\mathbb R}^N\otimes{\mathbb R}^d$ given by 
\begin{equation}\label{DDE}
(z_t')^i=(\partial/\partial\eta)|_{\eta=0}z_t^\eta.
\end{equation}
Then $(z'_t)_{t\in[0,1]}$ satisfies the following stochastic differential equation in ${\mathbb R}^N\otimes{\mathbb R}^d$
\begin{equation}\label{ZPT}
dz'_t=\sum_{\ell=1}^m\nabla W_\ell(z_t)z_t'dh_t^\ell+\sum_{\ell=1}^m\nabla Z_\ell(z_t)z_t'dB^\ell_t+\nabla Z_0(z_t)z'_tdt
+\sum_{\ell=1}^mZ_\ell(z_t)\otimes(v_tX_\ell(x_t))dt.
\end{equation}
Define processes $(m_t)_{t\in[0,1]}$ in ${\mathbb R}$ and $(r_t)_{t\in[0,1]}$ in ${\mathbb R}^d$ by
\begin{align*}
dm_t
&=\sum_{\ell=1}^m\dot h_t^\ell dB^\ell_t,\quad m_0=0\\
dr_t
&=\sum_{\ell=1}^mv_tX_\ell(x_t)dB^\ell_t,\quad r_0=0.
\end{align*}
Write $(x_t')_{t\in[0,1]}$ for the derived process associated to the stochastic differential equation \eqref{DGE} and set
$$
y_t^{(0)}=(x_{t\wedge t_1},\dots,x_{t\wedge t_k},x_t,v_t,m_t,r_t,x_t').
$$ 
Then $(y_t^{(0)})_{t\in[0,1]}$ satisfies a stochastic differential equation of the form \eqref{ZT}.
The stopped components are obtained by multiplying the coefficients of \eqref{DGE} by the indicator function $1_{[0,t_i]}(t)$.
For $n\ge0$, recursively, set $z^{(n)}_t=(y^{(0)}_t,\dots,y^{(n)}_t)$, 
note that $(z_t^{(n)})_{t\in[0,1]}$ satisfies a stochastic differential equation of the form \eqref{ZT} with graded Lipschitz coefficients,
define $((z^{(n)})'_t)_{t\in[0,1]}$ by solving the associated derived equation, write $(z^{(n)})_t'=((y^{(0)})_t',\dots,(y^{(n)})_t')$, and set
$$
y^{(n+1)}_t=(y^{(n)})_t'.
$$
Where a process is presented in components, for example $y^{(0)}_t$, we will write the corresponding decomposition of its derived process in the obvious way.
Thus 
$$
(y^{(0)})_t'=(x_{t\wedge t_1}',\dots,x_{t\wedge t_k}',x_t',v_t',m_t',r_t',x_t''). 
$$
It is straightforward to check that this notation is consistent when a given process is a component in two different autonomous processes.
By a standard calculation, we have $x_t'={\sigma} u_tc_t$, where $c_t$ is the Malliavin covariance matrix
\begin{equation}\label{CTM}
c_t=\sum_{\ell=1}^m\int_0^t(v_sX_\ell(x_s))\otimes(v_sX_\ell(x_s))ds.
\end{equation}
It is well known that, under the bracket condition \eqref{LAD}, the Malliavin covariance matrix $c_1$ is invertible and its inverse has moments of all orders.
The basic form of Bismut's integration by parts formula is the identity in ${\mathbb R}^d$
$$
{\mathbb E}(\nabla\phi(z_1)z_1')={\mathbb E}(\phi(z_1)r_1)
$$
valid for $C_b^1$ functions $\phi$ on ${\mathbb R}^N$ and for $(z_t)_{t\in[0,1]}$ and $(z_t')_{t\in[0,1]}$ satisfying \eqref{ZT} and \eqref{ZPT} respectively.
Set 
$$
\tilde x_1=(x_1-{\gamma}_1)/{\sigma},\quad \tilde x_1'=x_1'/{\sigma}=u_1c_1,\quad \tilde x_1''=x_1''/{\sigma}.
$$
Define random variables $y$ in $({\mathbb R}^d)^*\otimes({\mathbb R}^d)^*$ and $R$ in ${\mathbb R}$ by
$$
y=(\tilde x'_1)^{-1},\quad R=(F(x_1)+\tfrac12d(x,y)^2)/{\sigma}^2+m_1/{\sigma}
$$
and define $y'$ in $({\mathbb R}^d)^*\otimes({\mathbb R}^d)^*\otimes{\mathbb R}^d$ and $R'$ in ${\mathbb R}^d$ by
$$
y'=-y\tilde x_1''y,\quad R'=\nabla F(x_1)x_1'/{\sigma}^2+m_1'/{\sigma}.
$$
Fix $n\ge0$ and apply the integration by parts formula with $(z_t)_{t\in[0,1]}$ replaced by $(z^{(n)}_t)_{t\in[0,1]}$ 
and with $\phi(z_1)$ replaced by $f(\tilde x_1)y\phi(y,z_1^{(n)})e^{-R}$ to obtain
\begin{align*}
&{\mathbb E}(y\otimes(\nabla f(\tilde x_1)\tilde x_1')\phi(y,z_1^{(n)})e^{-R})
+{\mathbb E}(f(\tilde x_1)y'\phi(y,z_1^{(n)})e^{-R})
-{\mathbb E}(f(\tilde x_1)\phi(y,z_1^{(n)})(y\otimes R')e^{-R})\\
&\quad+{\mathbb E}(f(\tilde x_1)y\otimes(\nabla_y\phi(y,z_1^{(n)})y')e^{-R})
+{\mathbb E}(f(\tilde x_1)y\otimes(\nabla_z\phi(y,z_1^{(n)})(z^{(n)})_1')e^{-R})\\
&\quad\q={\mathbb E}(f(\tilde x_1)\phi(y,z_1^{(n)})(y\otimes r_1)e^{-R}).
\end{align*}
We assume here that $f$ is $C_b^1$ and that $\phi$ is polynomial in $y$, with coefficients $C^1$ in $z_1^{(n)}$ and of polynomial growth, along with their derivatives.
Then all factors in the integrands of the preceding formula have moments of all orders.
A straightforward approximation by $C_b^1$ functions then establishes the formula.
Define a linear map ${\tau}_i:({\mathbb R}^d)^*\otimes({\mathbb R}^d)^*\otimes{\mathbb R}^d\to{\mathbb R}$ by 
$$
{\tau}_i(e_j^*\otimes e^*_{i'}\otimes e_{j'})={\delta}_{ii'}{\delta}_{jj'}.
$$
Then
$$
{\tau}_i(y\otimes(\nabla f(\tilde x_1)\tilde x_1'))=\nabla_if(\tilde x_1)
$$
so 
\begin{equation}\label{IBP}
{\mathbb E}(\nabla_if(\tilde x_1)\phi(y,z_1^{(n)})e^{-R})={\mathbb E}(f(\tilde x_1)\nabla_i^*\phi(y,z_1^{(n+1)})e^{-R})
\end{equation}
where
\begin{align*}
\nabla_i^*\phi(y,z_1^{(n+1)})
&={\tau}_i(y\otimes r_1+y\tilde x''_1y+y\otimes R')\phi(y,z_1^{(n)})\\
&+{\tau}_i(y\otimes(\nabla_y\phi(y,z_1^{(n)})y\tilde x''_1y))-{\tau}_i(y\otimes(\nabla_z\phi(y,z_1^{(n)})(z^{(n)})_1').
\end{align*}
Take 
$$
\phi_0(y,z^{(0)}_1)=G(y,z^{(0)}_1)=G(\tilde x)=g(\tilde x_{t_1},\dots,\tilde x_{t_k}).
$$
We see inductively that \eqref{IBP} is valid for $\phi_n=\nabla^*_{i_n}\dots\nabla^*_{i_1}\phi_0$ for all $n\ge0$.
So we can iterate \eqref{IBP} to obtain, for any multi-index ${\alpha}=(i_1,\dots,i_n)$
$$
{\mathbb E}(\nabla^{\alpha} f(\tilde x_1)G(\tilde x)e^{-R})={\mathbb E}(f(\tilde x_1)(\nabla^*)^{\alpha} G(y,z^{(n)}_1)e^{-R}).
$$
We take $f(x)=e^{i\<\xi,x\>}$ to deduce that $|\xi^{\alpha}||\hat G_{\varepsilon}(\xi)|\leqslant C_{\varepsilon}({\alpha},G)$ where 
$$
C_{\varepsilon}({\alpha},G)={\mathbb E}(|(\nabla^*)^{\alpha} G(y,z^{(n)}_1)|e^{-R}).
$$
Now $y=c_1^{-1}v_1$ and $(\nabla^*)^{\alpha} G$ is of polynomial growth in $(y,z^{(n)}_1)$.
Given the estimate \eqref{BKE}, to prove \eqref{HRR}, it will suffice to show that, for $n=d+1$ and for all $p<\infty$, we have
\begin{equation}\label{ZPE}
\sup_{{\sigma}\in(0,1]}{\mathbb E}(|z^{(n)}_1|^p)<\infty
\end{equation}
and
\begin{equation}\label{CPE}
\sup_{{\sigma}\in(0,1]}{\mathbb E}(|c_1^{-1}|^p)<\infty.
\end{equation}

We already noted the availability of moment estimates of all orders for $z^{(n)}_1$ and that these are uniform in ${\sigma}\in[0,1]$ for the components
derived from $(v_t)_{t\in[0,1]}$ and $(r_t)_{t\in[0,1]}$.
Recall that $\tilde x_t=(x_t-{\gamma}_t)/{\sigma}$ and $R=(F(x_1)+d(x,y)^2/2)/{\sigma}^2+m_1/{\sigma}$. 
We will use first and second order mean value theorems to see that there is in fact no singularity as ${\sigma}\to0$ in these processes or any processes derived from them.
Consider, for ${\theta}\in[0,1]$, the stochastic differential equation in ${\mathbb R}^d$
$$
dx_t({\theta})=\sum_{\ell=1}^mX_\ell(x_t({\theta}))dh^\ell_t+\sum_{\ell=1}^m{\theta}{\sigma} X_\ell(x_t({\theta}))dB^\ell_t+{\theta}^2{\sigma}^2\tilde X_0(x_t({\theta}))dt,\quad x_0({\theta})=x.
$$
There exists a unique family of strong solutions $(x_t({\theta}))_{t\in[0,1]}$ which are almost surely jointly continuous in ${\theta}$ and $t$.
Moreover, the following derivatives then exist for all ${\theta}$ and $t$, almost surely
$$
\bar x_t({\theta})=\frac1{\sigma}\frac\partial{\partial{\theta}}x_t({\theta}),\quad\hat x_t({\theta})=\frac1{{\sigma}^2}\left(\frac\partial{\partial{\theta}}\right)^2x_t({\theta}).
$$
Moreover, the processes $(\bar x_t({\theta}))_{t\in[0,1]}$ and $(\hat x_t({\theta}))_{t\in[0,1]}$ satisfy
\begin{align*}
d\bar x_t({\theta})
&=\sum_{\ell=1}^m\nabla X_\ell(x_t({\theta}))\bar x_t({\theta})dh^\ell_t+\sum_{\ell=1}^m{\theta}{\sigma}\nabla X_\ell(x_t({\theta}))\bar x_t({\theta})dB^\ell_t
+{\theta}^2{\sigma}^2\nabla\tilde X_0(x_t({\theta}))\bar x_t({\theta})dt\\
&\quad\q+\sum_{\ell=1}^mX_\ell(x_t({\theta}))dB^\ell_t+2{\theta}{\sigma}\tilde X_0(x_t({\theta}))dt,\quad \bar x_0({\theta})=0\\
d\hat x_t({\theta})
&=\sum_{\ell=1}^m\nabla X_\ell(x_t({\theta}))\hat x_t({\theta})dh^\ell_t+\sum_{\ell=1}^m{\theta}{\sigma}\nabla X_\ell(x_t({\theta}))\hat x_t({\theta})dB^\ell_t
+{\theta}^2{\sigma}^2\nabla\tilde X_0(x_t({\theta}))\hat x_t({\theta})dt\\
&\quad+\sum_{\ell=1}^m\nabla^2X_\ell(x_t({\theta}))(\bar x_t({\theta}),\bar x_t({\theta}))dh^\ell_t
+\sum_{\ell=1}^m{\theta}{\sigma}\nabla^2X_\ell(x_t({\theta}))(\bar x_t({\theta}),\bar x_t({\theta}))dB^\ell_t\\
&\quad\q\quad\q\quad\q+{\theta}^2{\sigma}^2\nabla^2\tilde X_0(x_t({\theta}))(\bar x_t({\theta}),\bar x_t({\theta}))dt\\
&\quad+\sum_{\ell=1}^m2\nabla X_\ell(x_t({\theta}))\bar x_t({\theta})dB^\ell_t+4{\theta}{\sigma}\nabla\tilde X_0(x_t({\theta}))\bar x_t({\theta})dt\\
&\quad+2\tilde X_0(x_t({\theta}))dt,\quad\hat x_0({\theta})=0.
\end{align*}
Set $z_t({\theta})=(x_t({\theta}),\bar x_t({\theta}),\hat x_t({\theta}))$.
Note that the process $(z_t({\theta}))_{t\in[0,1]}$ is the solution of a stochastic differential equation with graded Lipschitz coefficients, 
and that the coefficient bounds of the graded structure are uniform in ${\theta}\in[0,1]$ and ${\sigma}\in[0,1]$.
Hence $(z_t({\theta}))_{t\in[0,1]}$ and its derived process $(z_t'({\theta}))_{t\in[0,1]}$ have moments of all orders, and these are bounded uniformly in ${\theta}$ and ${\sigma}$.
Now 
$$
\tilde x_t=(x_t-{\gamma}_t)/{\sigma}=\int_0^1\bar x_t({\theta})d{\theta}
$$
and
\begin{equation}\label{EFR}
R=\int_0^1(1-{\theta})\{F'(x_1({\theta}))\hat x_1({\theta})+F''(x_1({\theta}))(\bar x_1({\theta}),\bar x_1({\theta}))\}d{\theta}.
\end{equation}
We can take the derivative in \eqref{DDE} under the integral sign in \eqref{EFR} to express $R'$ also by such an integral over ${\theta}$, 
with integrand expressed in terms of $z_1({\theta})$ and $z_1'({\theta})$.
Hence $\sup_{t\in[0,1]}|\tilde x_t|$, $R$ and $R'$ have moments of all orders which are bounded uniformly in ${\sigma}\in(0,1]$.
The same reasoning can be extended to conclude that all processes derived from $(\tilde x_t)_{t\in[0,1]}$ and $R$ 
have moments of all orders bounded uniformly in ${\sigma}\in(0,1]$.
Alternatively, it can be observed that $(x_t,\tilde x_t',v_t/{\sigma},\tilde x_t''/{\sigma},m_t''/{\sigma})$ satisfies a stochastic differential equation 
with graded Lipschitz coefficients, with bounds uniform in ${\sigma}\in(0,1]$, and from this observation we can draw the same conclusion.
Hence we have shown \eqref{ZPE}.

We turn to the proof of \eqref{CPE}.
Write $(v_t(0))_{t\in[0,1]}$ for the solution to \eqref{VGE} when ${\sigma}=0$.
Recall that
$$
c_1=\sum_{\ell=1}^m\int_0^1(v_tX_\ell(x_t))\otimes(v_tX_\ell(x_t))dt
$$
and note that $c_1$ depends continuously on $(x_t,v_t)_{t\in[0,1]}$ in uniform norm.
We have assumed that the deterministic Malliavin covariance matrix $c_1(0)$ is invertible.
Hence there are constants $C<\infty$ and ${\delta}>0$ such that $|c_1^{-1}|\leqslant C$ on the event
$$
{\Omega}({\delta})=\{|x_t-{\gamma}_t|\leqslant{\delta}\quad\text{and}\quad|v_t-v_t(0)|\leqslant{\delta}\quad\text{for all}\quad t\in[0,1]\}.
$$
By standard $L^p$ estimates for stochastic differential equations, for all $p<\infty$, there is a constant $C({\delta},p)<\infty$ such that, for all ${\sigma}\in(0,1]$,
$$
\mathbb{P}({\Omega}({\delta})^c)\leqslant C({\delta},p){\sigma}^p.
$$
On the other hand, it is shown in \cite{AKS} that, for all $p<\infty$, there are constants $C(p)<\infty$ and $D\in{\mathbb N}$ such that, for all ${\sigma}\in(0,1]$,
$$
\|c_1^{-1}\|_p\leqslant C(p){\sigma}^{-D}.
$$
Now $|c_1^{-1}|\leqslant C+|c_1^{-1}|1_{{\Omega}({\delta})^c}$ so, by H\"older's inequality, for all $p<\infty$ and all ${\sigma}\in(0,1]$,
$$
\|c_1^{-1}\|_p\leqslant C+C(2p)C({\delta},2Dp)^{1/2p}.
$$

Finally, consider the case where $G({\omega})=|{\omega}_s-{\omega}_t|^4$ for some $s,t\in[0,1]$.
Set $\tilde x_t^{(0)}=\tilde x_t$ and, recursively for $n\ge0$, set $\tilde x_t^{(n+1)}=(\tilde x_t,(\tilde x^{(n)})'_t)$.
Then, by standard estimates, for all $p\in[1,\infty)$, there is a constant $C<\infty$ such that, uniformly in $s,t\in[0,1]$ and in ${\sigma}\in(0,1]$, we have
$$
{\mathbb E}\left(|\tilde x_s^{(n)}-\tilde x_t^{(n)}|^{4p}\right)\leqslant C|s-t|^{2p}.
$$
The adjoint operators $\nabla_i^*$ have an explicit form written above, from which we deduce that, for all $n$ and for ${\alpha}=(i_1,\dots,i_n)$,
there is a random variable $K_{\alpha}$, having moments of all orders bounded uniformly in ${\sigma}\in(0,1]$, such that
$$
(\nabla^*)^{\alpha} G(y,z_1^{(n)})=K_{\alpha}|\tilde x_s^{(n)}-\tilde x_t^{(n)}|^4.
$$
Hence, by H\"older's inequality, there is a constant $C_{\alpha}<\infty$ such that, uniformly in $s,t\in[0,1]$ and in ${\varepsilon}\in(0,1]$, we have
$$
C_{\varepsilon}({\alpha},G)\leqslant C_{\alpha}|s-t|^2.
$$
This implies \eqref{GST}.
\end{proof}

\section{Localization}\label{LL}
Having established the conclusion of Theorem \ref{MR} under the global condition \eqref{URD}, 
we now show that it extends to general manifolds under a condition \eqref{LOC} which limits the probability that the bridge strays too far from the minimal path. 
At the same time we obtain that this condition allows the extension of the small-time equivalent for the heat kernel.
In the next section, we will find conditions on the operator which imply \eqref{LOC}.

\begin{theorem}\label{NFBT}
Let $M$ be a connected $C^\infty $ manifold.
Let ${\cal L}$ be a second order differential operator on $M$ of the form \eqref{LAB}
$$
{\cal L}=\frac12\sum_{i,j=1}^da^{ij}(x)\frac{\partial^2}{\partial x^i\partial x^j}+\sum_{i=1}^db^i(x)\frac{\partial}{\partial x^i}
$$
where the diffusivity $a$ has a sub-Riemannian structure $(X_1,\dots,X_m)$.
Suppose that $x,y\in M$ are joined by a unique minimal path ${\gamma}$ and that $(x,y)$ lies outside the cut locus.
Suppose further that there is a chart $\phi:U_0\to{\mathbb R}^d$ of $M$ and a domain $U$ in $M$, containing ${\gamma}$ and compactly contained in $U_0$, such that
\begin{equation}\label{LOC}
\limsup_{{\varepsilon}\to0}{\varepsilon}\log p({\varepsilon},x,M\setminus U,y)<-d(x,y)^2/2
\end{equation}
where $p({\varepsilon},x,M\setminus U,y)=p({\varepsilon},x,y)-p_U({\varepsilon},x,y)$ and $p_U({\varepsilon},x,y)$ is the Dirichlet heat kernel in $U$ of the restriction of ${\cal L}$ to $U$.
Then, as ${\varepsilon}\to0$, we have
$$
p({\varepsilon},x,y)=c(x,y){\varepsilon}^{-d/2}\operatorname{exp}\{-d(x,y)^2/(2{\varepsilon})\}(1+o(1))
$$
and $\tilde\mu^{x,y}_\ve\to\mu_\g$ weakly on ${T_\g\oxy}$.
\end{theorem}
\begin{proof}
It will suffice, and it will lighten the notation, to consider the case where $U_0$ is a domain both in $M$ and in ${\mathbb R}^d$, 
where $\phi$ is the identity map, and where ${\theta}$ restricts to the identity map on $U$.
We will show that there exist vector fields $\bar X_0,\bar X_1,\dots,\bar X_{m+d}$ on ${\mathbb R}^d$, which are bounded with bounded derivatives of all orders and such that 
\begin{itemize} 
\item[{\rm (a)}] $(\bar X_1,\dots,\bar X_{m+d})$ is a sub-Riemannian structure on ${\mathbb R}^d$, 
\item[{\rm (b)}] we have $a=\bar a$ and ${\cal L}=\bar{\cal L}$ on $U$, where
$$
\bar{\cal L}=\frac12\sum_{\ell=1}^{m+d}\bar X_\ell^2+\bar X_0
$$
and where $a$ and $\bar a$ are the diffusivities of ${\cal L}$ and $\bar{\cal L}$ respectively,
\item[{\rm (c)}] ${\gamma}$ is the unique $\bar a$-minimal path in $H^{x,y}({\mathbb R}^d)$.
\end{itemize}
From these properties, given that $(x,y)$ lies outside the cut locus of ${\cal L}$ in $M$, we deduce that $(x,y)$ also lies outside the cut locus of $\bar{\cal L}$ in ${\mathbb R}^d$.
Write $\bar p$ for the heat kernel of $\bar{\cal L}$ on ${\mathbb R}^d$ and write $\mu^{x,y,\R^d}_\ve$ for the bridge measure on $\Omega^{x,y}({\mathbb R}^d)$ associated to the operator ${\varepsilon}\bar{\cal L}$.
Similarly, write $\tilde\mu^{x,y,\R^d}_\ve$ for the rescaled bridge measure on ${T_\g\oxy}$, given by $\tilde\mu^{x,y,\R^d}_\ve=\mu^{x,y,\R^d}_\ve\circ\bar{\sigma}_{\varepsilon}^{-1}$, where $\bar{\sigma}_{\varepsilon}({\omega})_t=({\omega}_t-{\gamma}_t)/\sqrt{\varepsilon}$.
Then, as we showed in Section \ref{LAP}, as ${\varepsilon}\to0$, we have 
$$
\bar p({\varepsilon},x,y)=c(x,y){\varepsilon}^{-d/2}\operatorname{exp}\{-d(x,y)^2/(2{\varepsilon})\}(1+o(1))
$$
and $\mu^{x,y,\R^d}_\ve\to\d_\g$ weakly on $\Omega^{x,y}({\mathbb R}^d)$ and $\tilde\mu^{x,y,\R^d}_\ve\to\mu_\g$ weakly on ${T_\g\oxy}$.

Write $\mu^{x,y,U}_\ve$ for the diffusion bridge measure on $\Omega^{x,y}(U)$ associated to the restriction of the operator ${\varepsilon}{\cal L}$ to $U$.
Then, for all measurable sets $A$ in $\Omega^{x,y}(M)$, we have
$$
p({\varepsilon},x,y)\mu^{x,y}_\ve(A)=p_U({\varepsilon},x,y)\mu^{x,y,U}_\ve(A\cap\Omega^{x,y}(U))+p({\varepsilon},x,y)\mu^{x,y}_\ve(A\setminus\Omega^{x,y}(U))
$$
and
$$
\bar p({\varepsilon},x,y)\mu^{x,y,\R^d}_\ve(A\cap\Omega^{x,y}(U))=p_U({\varepsilon},x,y)\mu^{x,y,U}_\ve(A\cap\Omega^{x,y}(U)).
$$
For any bounded measurable set $B$ in ${T_\g\oxy}$ and for ${\varepsilon}>0$ sufficiently small, we have ${\sigma}_{\varepsilon}^{-1}(B)=\bar{\sigma}_{\varepsilon}^{-1}(B)\subseteq\Omega^{x,y}(U)$, so
$$
p({\varepsilon},x,y)\tilde\mu^{x,y}_\ve(B)=\bar p({\varepsilon},x,y)\tilde\mu^{x,y,\R^d}_\ve(B).
$$
Now $\mu^{x,y,\R^d}_\ve(\Omega^{x,y}(U))\to1$ so, taking $A=\Omega^{x,y}(M)$, we find
$$
p_U({\varepsilon},x,y)=c(x,y){\varepsilon}^{-d/2}\operatorname{exp}\{-d(x,y)^2/(2{\varepsilon})\}(1+o(1))
$$
and, using \eqref{LOC}, 
$$
p({\varepsilon},x,y)=p_U({\varepsilon},x,y)+p({\varepsilon},x,M\setminus U,y)=c(x,y){\varepsilon}^{-d/2}\operatorname{exp}\{-d(x,y)^2/(2{\varepsilon})\}(1+o(1)).
$$
Moreover, for all $B\in{T_\g\oxy}$, we then have $\tilde\mu^{x,y}_\ve(B)=\tilde\mu^{x,y,\R^d}_\ve(B)(1+o(1))$, so the claimed weak limit for $\tilde\mu^{x,y}_\ve$ follows from that for $\tilde\mu^{x,y,\R^d}_\ve$.

It remains to show the existence of vector fields $\bar X_0, \bar X_1,\dots,\bar X_{m+d}$ with the claimed properties.
Fix an open set $U_1$ such that $U$ is compactly contained in $U_1$ and $U_1$ is compactly contained in $U_0$.
There exists a $C^\infty $ function $\chi$ on ${\mathbb R}^d$ such that $1_U\leqslant\chi\le1$ and which is uniformly positive on $U_1$ and such that $\{\chi>0\}$ is compactly contained in $U_0$.
Write $X_0$ for the vector field on $M$ given by \eqref{HSS}.
For $\ell=0,1,\dots,m$, set $\bar X_\ell(z)=\chi(z)X_\ell(z)$ for $z\in U_0$ and set $\bar X_\ell(z)=0$ for $z\in{\mathbb R}^d\setminus U_0$.
Then $\bar X_\ell$ is a bounded vector field on ${\mathbb R}^d$ with bounded derivatives of all orders.
Define $\tilde a=\sum_{\ell=1}^m\bar X_\ell\otimes\bar X_\ell$ and write $\tilde I$ for the associated energy function.
Then $\tilde I({\omega})\geqslant I({\omega})$ for all ${\omega}\inH^{x,y}(U_0)$ with equality if ${\omega}$ is contained in $U$.
Moreover, by choosing $\chi$ sufficiently small near $\partial U_1$, we can and do ensure that, for all $z\in\partial U_1$ and all ${\omega}\in H^{x,z}(U_0)$, 
we have $\tilde I({\omega})>I({\gamma})$.
Choose now another $C^\infty $ function $\bar\chi$ on ${\mathbb R}^d$ such that $1_{U_1}\le1-\bar\chi\le1_{U_0}$ and such that $\chi+\bar\chi$ is everywhere positive.
For $i=1,\dots,d$, set $\bar X_{m+i}(z)=\bar\chi(z)e_i$, where $e_1,\dots,e_d$ is the standard basis in ${\mathbb R}^d$.
Then the brackets of $\bar X_1,\dots,\bar X_m$ span ${\mathbb R}^d$ on $\{\chi>0\}$, while $\bar X_{m+1},\dots,\bar X_{m+d}$ themselves span ${\mathbb R}^d$ on $\{\bar\chi>0\}$.
Hence (a) holds.
Also, (b) holds because $\bar X_\ell=X_\ell$ on $U$ for $\ell=0,1,\dots,m$ and $\bar X_\ell=0$ on $U$ for $\ell=m+1,\dots,m+d$.
Now, if ${\omega}\inH^{x,y}({\mathbb R}^d)$ is contained in $U_1$, then $\bar I({\omega})=\tilde I({\omega})\geqslant I({\omega})$, so $\bar I({\omega})\geqslant I({\gamma})$ with equality only if ${\omega}={\gamma}$.
On the other hand, if ${\omega}$ is not contained in $U_1$, set ${\tau}=\inf\{t\in[0,1]:{\omega}_{\tau}\in\partial U_1\}$ and set $\tilde{\omega}_t={\omega}_{{\tau} t}$.
Then $\tilde{\omega}\in H^{x,z}(U_0)$, where $z={\omega}_{\tau}\in\partial U_1$.
So $\bar I({\omega})\geqslant \tilde I(\tilde{\omega})>I({\gamma})$.
Hence (c) holds.
\end{proof}

\section{Heat kernel upper bounds}\label{NFB}
Let $M$ be a connected $C^\infty $ manifold.
Consider a second order differential operator ${\cal L}$ on $M$ of the form \eqref{CFO}
$$
{\cal L} f=\tfrac12\operatorname{div}(a\nabla f)+a({\beta},\nabla f)
$$
where the diffusivity $a$ has a sub-Riemannian structure $(X_1,\dots,X_m)$,
where the divergence is understood with respect to a positive $C^\infty $ locally invariant measure $\nu$ and where ${\beta}$ is a $C^\infty $ $1$-form satisfying the sector condition \eqref{SECT}.
Recall the definition \eqref{DIST} of the sub-Riemannian distance.
For a closed subset $K$ of $M$ and $x,y\in M\setminus K$, consider the distance from $x$ to $y$ through $K$, given by
$$
d(x,K,y)=\inf\{\sqrt{I({\omega})}:{\omega}\inH^{x,y},\,{\omega}_t\in K\text{ for some $t\in[0,1]$}\}.
$$
The heat kernel for diffusion from $x$ to $y$ through $K$ is defined by
$$
p(t,x,K,y)=p(t,x,y)-p_{M\setminus K}(t,x,y)
$$
where $p_{M\setminus K}$ is the Dirichlet heat kernel of ${\cal L}$ in $M\setminus K$.
We will prove in this section the following upper bound.\footnote{%
In the case where $a$ is everywhere positive-definite, there is an alternative argument of Hsu \cite{MR1089046} for heat kernel upper bounds in incomplete manifolds, 
relying  on estimates of Azencott \cite{MR634964}, which does not require the condition \eqref{SECT}. 
Hsu instead requires that $d(x,y)\leqslant d(x,\infty)+d(\infty,y)$. 
The same paper \cite{MR1089046} also contains an example where heat flow does not localize near the shortest path.}
\begin{theorem}\label{HTK}
For all $x,y\in M$ and all closed subsets $K$ of $M$, we have
$$
\limsup_{t\to0}t\log p(t,x,K,y)\leqslant-d(x,K,y)^2/2.
$$
\end{theorem}

Before giving the proof, we show how this estimate allows us to prove Theorems \ref{VAR}, \ref{BALOC} and \ref{MR}.

\begin{proof}[Proof of Theorem \ref{VAR}]
The upper bound in \eqref{LOGP} is the simple case $K=M$ of Theorem \ref{HTK}, 
while the lower bound follows by a standard argument from L\'eandre's lower bound in ${\mathbb R}^d$.

We turn to the proof of \eqref{LDEV}. 
Consider the operator 
$$
\tilde{\cal L}={\cal L}+\tfrac12(\partial/\partial{\tau})^2
$$ 
on $\tilde M=M\times{\mathbb R}$, where ${\tau}$ denotes the coordinate in ${\mathbb R}$.
Then the diffusivity of $\tilde{\cal L}$ also has a sub-Riemannian structure.
Set $\tilde x=(x,0)$ and $\tilde y=(y,1)$, and define
$$
\tilde K=\tilde M\setminus\tilde U,\quad \tilde U=\{({\gamma}_t,{\sigma}_t):({\gamma},{\sigma})\in\tilde{\Gamma}({\delta}),\,t\in[0,1]\}
$$
where
$$
\tilde{\Gamma}({\delta})=\left\{({\gamma},{\sigma})\inH^{x,y}\times H^{0,1}({\mathbb R}):I({\gamma})+\int_0^1\dot{\sigma}_t^2dt<d(x,y)^2+1+{\delta}\right\}.
$$
Then $\tilde K$ is closed in $\tilde M$.
Write ${\beta}^{0,1}_{\varepsilon}$ for the law on ${\Omega}^{0,1}({\mathbb R})$ of a Brownian bridge from $0$ to $1$ of speed ${\varepsilon}$.
Then, with obvious notation,
$$
\tilde p(t,\tilde x,\tilde y)=p(t,x,y)\frac1{\sqrt{2\pi}}e^{-1/(2t)},\quad\tilde\mu^{\tilde x,\tilde y}_{\varepsilon}(d{\omega},d{\tau})=\mu^{x,y}_\ve(d{\omega}){\beta}_{\varepsilon}^{0,1}(d{\tau}).
$$
By Theorem \ref{HTK}, we have
$$
\limsup_{t\to0}t\log\tilde p(t,\tilde x,\tilde K,\tilde y)\leqslant-\tilde d(\tilde x,\tilde K,\tilde y)^2/2=-(d(x,y)^2+1+{\delta})/2
$$
so
\begin{align}\notag
&\limsup_{{\varepsilon}\to0}{\varepsilon}\log\tilde\mu_{\varepsilon}^{\tilde x,\tilde y}(\{({\omega},{\tau}):({\omega}_t,{\tau}_t)\in\tilde K\text{ for some }t\in[0,1]\})\\
\label{LD1}
&\quad\q\quad\q\leqslant\limsup_{{\varepsilon}\to0}{\varepsilon}\log\tilde p({\varepsilon},\tilde x,\tilde K,\tilde y)-\operatorname{liminf}_{{\varepsilon}\to0}{\varepsilon}\log\tilde p({\varepsilon},\tilde x,\tilde y)\leqslant-{\delta}/2
\end{align}
where we have used the lower bound from \eqref{LOGP}.
By standard estimates, we also have
\begin{equation}\label{LD2}
\lim_{{\varepsilon}\to0}{\varepsilon}\log{\beta}_{\varepsilon}^{0,1}(\{{\tau}:|{\tau}_t-t|\geqslant\sqrt{\delta}/2\text{ for some }t\in[0,1]\})=-{\delta}/2.
\end{equation}
Suppose then that ${\omega}\in\Omega^{x,y}$ and ${\tau}\in{\Omega}^{0,1}({\mathbb R})$ satisfy $({\omega}_t,{\tau}_t)\in\tilde U$ and $|{\tau}_t-t|<\sqrt{\delta}/2$ for all $t\in[0,1]$.
Then, for each $t\in[0,1]$, there exist $s\in[0,1]$ and ${\gamma}\inH^{x,y}$ and ${\sigma}\in H^{0,1}({\mathbb R})$ such that 
$$
{\omega}_t={\gamma}_s,\quad {\tau}_t={\sigma}_s,\quad I({\gamma})<d(x,y)^2+{\delta},\quad \int_0^1\dot{\sigma}_r^2dr<1+{\delta}.
$$
Then $|{\sigma}_s-s|\leqslant\sqrt{\delta}/2$ so $|t-s|\leqslant\sqrt{\delta}$ and so 
$$
d({\omega}_t,{\Gamma}_t({\delta}))^2\leqslant d({\omega}_t,{\gamma}_t)^2=d({\gamma}_s,{\gamma}_t)^2\leqslant|t-s|I({\gamma})\leqslant{\delta}^{1/2}(d(x,y)^2+{\delta}).
$$
The estimates \eqref{LD1} and \eqref{LD2} thus imply \eqref{LDEV}.
\end{proof}  

\begin{proof}[Proof of Theorems \ref{BALOC} and \ref{MR}]
The hypothesis that ${\gamma}\inH^{x,y}$ is strongly minimal implies that, for any domain $U$ containing ${\gamma}$, we have $d(x,M\setminus U,y)>d(x,y)$.
We can choose a chart $U_0$ in $M$ and a domain $U$ containing ${\gamma}$ such that $U$ is compactly contained in $U_0$.
We suppose that ${\cal L}$ has the form \eqref{CFO} and satisfies the sector condition \eqref{SECT}.
Then Theorem \ref{HTK} applies to show that
$$
\limsup_{t\to0}t\log p(t,x,M\setminus U,y)\leqslant-d(x,M\setminus U,y)^2/2<-d(x,y)^2/2
$$
and Theorem \ref{NFBT} then gives the desired conclusions.
\end{proof}

In proving Theorem \ref{HTK}, we will use the following estimate of Nagel, Stein \& Wainger \cite{MR793239} for the sub-Riemannian distance.
There is a covering of $M$ by relatively compact charts $\phi:U\to{\mathbb R}^d$ such that, for some constants ${\alpha}(U)>0$ and $C(U)<\infty$, for all $x,y\in U$, we have
$$
C^{-1}|\phi(x)-\phi(y)|\leqslant d(x,y)\leqslant C|\phi(x)-\phi(y)|^{\alpha}.
$$
If follows that, for any positive $C^\infty $ measure $\nu$ on $M$, we can choose $C(U)$ so that moreover, 
for all $x\in U$ and all $r\in(0,\infty)$ such that $B(x,r)=\{y\in M:d(x,y)<r\}\subseteq U$, we have
\begin{equation}\label{VE}
C^{-1}r^{d/{\alpha}}\leqslant\nu(B(x,r))\leqslant Cr^d
\end{equation}
In particular, $C(U)$ may be chosen so that the following local volume-doubling inequality holds: 
for all $x\in U$ and all $r\in(0,\infty)$ such that $B(x,2r)\subseteq U$, we have
\begin{equation}\label{DP}
\nu(B(x,2r))\leqslant C\nu(B(x,r)).
\end{equation}
We will also need the local Poincar\'e inequality proved by Jerison \cite{MR850547}, which also relies on the fact that $a$ has a sub-Riemannian structure.
There is a covering of $M$ by open sets $U$ such that, for some constant $C(U)<\infty$, 
for all $x\in U$ and all $r\in(0,\infty)$ such that $B(x,2r)\subseteq U$, for all $f\inC^\infty _c(M)$, we have
\begin{equation}\label{PI}
\int_{B(x,r)}|f-\<f\>_{B(x,r)}|^2d\nu\leqslant Cr^2\int_{B(x,2r)}a(\nabla f,\nabla f)d\nu
\end{equation}
where $\<f\>_B=\int_Bfd\nu/\nu(B)$ is the average value of $f$ on $B$.
The sector condition \eqref{SECT} gives us uniform parabolicity in the sense of Sturm, so \eqref{DP} and \eqref{PI} allow us to apply \cite[Theorem 2.1]{MR1355744}
to obtain the following parabolic mean-value estimate.
For simplicity, we will state a simple version adequate for our needs.
Let $x\in M$ and let $r_0\in(0,\infty)$.
Let $D$ be a domain in $M$ which compactly contains the ball $B(x,r_0)$.
Suppose that $u$ is a non-negative weak solution of the heat equation $(\partial/\partial t)u_t={\cal L} u_t$ on $(0,\infty)\times D$.
Then there is constant $C(D,r_0)<\infty$ such that, for all $t\in(0,\infty)$ and all $r\in(0,r_0]$ with $r^2\leqslant t/2$, we have
\begin{equation}\label{MVE}
u_t(x)^2\leqslant C\fint_{t-r^2}^t\fint_{B(x,r)}u_s^2d\nu ds.
\end{equation}
See Lierl and Saloff-Coste \cite[Theorem 2.6]{1205.6493} for an alternative proof, which clarifies some points in the non-symmetric case.

The sub-Riemannian distance has a dual formulation, proved in Jerison \& Sanchez-Calle \cite{MR922334}, 
which we adapt to the case of a general sub-Riemannian manifold, without completeness, and to the distance through $K$.
\begin{proposition}\label{DUAL}
For all $x,y\in M$ and any closed subset $K$ of $M$, we have
$$
d(x,K,y)=\sup\{w^+(y)-w^-(x):w^-,w^+\in{\cal F}\text{ with $w^+=w^-$ on $K$}\}
$$
where ${\cal F}$ is the set of all locally Lipschitz functions $w$ on $M$ such that $a(\nabla w,\nabla w)\le1$ almost everywhere.
\end{proposition}
\begin{proof}
Denote the right hand side by ${\delta}(x,K,y)$ for now.
Suppose that ${\omega}\inH^{x,y}$ has driving path $\xi$ and that ${\omega}_t\in K$.
Let $w^-,w^+\in{\cal F}$, with $w^+=w^-$ on $K$.
It will suffice to consider the case where ${\omega}$ is simple, when there exists a relatively compact chart $U$ for $M$ containing ${\omega}$.
Then, given ${\varepsilon}>0$, since $a$ is continuous, we can find $C^\infty $ functions $f^-,f^+$ on $U$ such that $|f^\pm(z)-w^\pm(z)|\leqslant{\varepsilon}$ and $a(\nabla f^\pm,\nabla f^\pm)(z)\le1+{\varepsilon}$ for all $z\in U$.
Then
$$
w^+(y)-w^-(x)=w^+(y)-w^+({\omega}_t)+w^-({\omega}_t)-w^-(x)\leqslant f^+(y)-f^+({\omega}_t)+f^-({\omega}_t)-f^-(x)+4{\varepsilon}
$$
and
\begin{align*}
&f^+(y)-f^+({\omega}_t)+f^-({\omega}_t)-f^-(x)\\
&\quad\q=\int_0^t\<\nabla f^-({\omega}_s),\dot{\omega}_s\>ds+\int_t^1\<\nabla f^+({\omega}_s),\dot{\omega}_s\>ds\\
&\quad\q=\int_0^t\<\nabla f^-({\omega}_s),a({\omega}_s)\xi_s\>ds+\int_t^1\<\nabla f^+({\omega}_s),a({\omega}_s)\xi_s\>ds\\
&\quad\q\leqslant\left(\int_0^ta(\nabla f^-,\nabla f^-)({\omega}_s)ds+ \int_t^1a(\nabla f^+,\nabla f^+)({\omega}_s)ds\right)^{1/2}\left(\int_0^1a(\xi_s,\xi_s)ds\right)^{1/2}\\
&\quad\q\leqslant(1+{\varepsilon})\sqrt{I({\omega})}.
\end{align*}
Hence $w^+(y)-w^-(x)\leqslant\sqrt{I({\omega})}$.
On taking the supremum over $w^\pm$ and the infimum over ${\omega}$, we deduce that ${\delta}(x,K,y)\leqslant d(x,K,y)$.

Now we prove the reverse inequality.
Choose a $C^\infty $ quadratic form ${\alpha}$ on $T^*M$ which is everywhere positive-definite.
Consider the energy function $I_{\alpha}$ and the distance functions $d_{\alpha}$ and ${\delta}_{\alpha}$ obtained by replacing $a$ by $a+{\alpha}$ in the definitions of $d$ and ${\delta}$.
Set $w^+(z)=d_{\alpha}(x,K,z)$ and $w^-(z)=d_{\alpha}(x,z)$.
Note that $w^+=w^-$ on $K$.
Since $a+{\alpha}$ is positive-definite, the functions $w^-$ and $w^+$ are locally Lipschitz, and their weak gradients $\nabla w^\pm$ satisfy $(a+{\alpha})(\nabla w^\pm,\nabla w^\pm)\le1$ almost everywhere.
Hence 
$$
d_{\alpha}(x,K,y)=w^+(y)-w^-(x)\leqslant{\delta}_{\alpha}(x,K,y)\leqslant{\delta}(x,K,y)
$$
and we can complete the proof by finding, for each ${\varepsilon}\in(0,1]$, a choice of ${\alpha}$ so that $d(x,K,y)\leqslant(1+{\varepsilon})(d_{\alpha}(x,K,y)+{\varepsilon})+{\varepsilon}$.

Choose $\phi^{\alpha}\in\Omega^{x,y}$, passing through $K$ and such that $\sqrt{I_{\alpha}(\phi^{\alpha})}\leqslant d_{\alpha}(x,K,y)+{\varepsilon}$.
There exist $p\in{\mathbb N}$ and $C^\infty $ vector fields $Y_1,\dots,Y_p$ on $M$ such that ${\alpha}=\sum_{i=1}^pY_i\otimes Y_i$.
Then there exist $h\inH^0(\R^m)$ and $k\inH^0(\R^p)$ such that $\phi^{\alpha}$ has weak derivative\footnote{%
See the discussion preceding Proposition \ref{2.1}.}
$$
\dot\phi_t^{\alpha}=\sum_{\ell=1}^mX_\ell(\phi^{\alpha}_t)\dot h^\ell_t+\sum_{i=1}^pY_i(\phi_t^{\alpha})\dot k_t^i.
$$
By reparametrizing $\phi^{\alpha}$ if necessary, we may assume that $|\dot h_t|^2+|\dot k_t|^2=I_{\alpha}(\phi^{\alpha})$ for almost all $t$.
It will be convenient to assume that $d(x,\infty)\ge1/2$, which we may do without loss of generality by a scaling argument.
Define $S_{-1}=S_0=\emptyset$ and consider for $n\ge1$ the compact set
$$
S_n=\{z\in M:d(z,x)\leqslant 2^n\text{ and }d(z,\infty)\ge2^{-n}\}.
$$
Define, recursively, sequences of times $t_0,t_1,\dots,t_J$, points $x_1,\dots,x_J$ and $y_1,\dots,y_J$, and positive integers $n_1,\dots,n_J$ as follows.  
Set $t_0=0$ and $n_1=1$.
For $j\ge1$, set $x_j=\phi^{\alpha}_{t_{j-1}}$, define $(\phi^j_t)_{t\geqslant t_{j-1}}$ by the differential equation
$$
\dot\phi_t^j=\sum_{\ell=1}^mX_\ell(\phi_t^j)\dot h^\ell_t,\quad\phi_{t_{j-1}}^j=x_j
$$
and set $y_j=\phi_{t_j}^j$, where
$$
t_j=\inf\{t\geqslant t_{j-1}:\phi^j_t\in\partial S_{n_j-1}\cup\partial S_{n_j+1}\}\wedge1.
$$
If $t_j<1$, then set 
$$
n_{j+1}=
\begin{cases}
n_j-1,&\text{if }y_j\in\partial S_{n_j-1},\\
n_j+1,&\text{if }y_j\in\partial S_{n_j+1}.
\end{cases}
$$
Thus $y_j\in\partial S_{n_{j+1}}$.
If $t_j=1$, then set $J=j$ and stop.
We will see below that the recursion does stop.
Note that, for all $j$ and all $s,t\in[t_{j-1},t_j]$ with $s\leqslant t$, we have
$$
d(\phi_s^j,\phi_t^j)\leqslant(t-s)\|\dot h\|_\infty.
$$
Note also that $x_1=\phi^{\alpha}_0=x\in S_{n_1}\subseteq S_{n_1+1}\setminus S_{n_1-1}$ and $\phi_{t_{j-1}}^{\alpha}=\phi_{t_{j-1}}^j=x_j$ for all $j$.
Suppose, inductively for $j\ge1$, that $x_j\in S_{n_j+1}\setminus S_{n_j-1}$.
Then $\phi^j_t$ remains in this set for $t\in[t_{j-1},t_j)$ and, by choosing ${\alpha}$ sufficiently small on $S_{n_j+2}\setminus S_{n_j-2}$, we can ensure that 
$$
d(\phi_t^j,\phi_t^{\alpha})\le2^{-2n_j-2}{\varepsilon}\le2^{-n_{j+1}-2},\quad t\in[t_{j-1},t_j]. 
$$
Then, in particular, we have $d(y_j,x_{j+1})\leqslant 2^{-n_{j+1}-2}$.
Now $y_j\in\partial S_{n_{j+1}}$ so this implies that $x_{j+1}\in S_{n_{j+1}+1}\setminus S_{n_{j+1}-1}$ and the induction proceeds.
Also, since $y_{j+1}\in\partial S_{n_{j+2}}$, we have $d(y_j,y_{j+1})\ge2^{-n_{j+1}-1}$ and so $d(x_{j+1},y_{j+1})\ge2^{-n_{j+1}-2}$.
Certainly $d(x_1,y_1)\ge2^{-n_1-2}$ so, for all $j\ge1$, we have $d(x_j,y_j)\ge2^{-n_j-2}$.
Hence
$$
(t_j-t_{j-1})\|\dot h\|_\infty\ge2^{-n_j-2}
$$
and so, for all $n\ge1$, we have $n_j=n$ at most $2^{n+2}\|\dot h\|_\infty$ times.
In particular, the recursion must stop, or $n_j\to\infty$ as $j\to\infty$ forcing $\phi_{t_j}^{\alpha}$ to leave all compact sets, which is impossible.
Hence
$$
\sum_{j=1}^{J-1}d(y_j,x_{j+1})+d(y_J,y)\leqslant\sum_{n=1}^\infty2^{n+2}\|\dot h\|_\infty.2^{-2n-2}{\varepsilon}={\varepsilon}\|\dot h\|_\infty.
$$
Also, there exist $j$ and $t$ so that $t\in[t_{j-1},t_j]$ and $\phi_t^{\alpha}\in K$, and then
$$
d(x_j,K,y_j)\leqslant d(x_j,\phi_t^j)+2d(\phi_t^j,\phi_t^{\alpha})+d(\phi_t^j,y_j)\leqslant(t_j-t_{j-1})\|\dot h\|_\infty+{\varepsilon}
$$
while, for $k\not=j$, we have the estimate
$$
d(x_k,y_k)\leqslant(t_k-t_{k-1})\|\dot h\|_\infty.
$$
We finally combine these estimates, using the triangle inequality, to obtain
$$
d(x,K,y)\leqslant(1+{\varepsilon})\|\dot h\|_\infty+{\varepsilon}\leqslant(1+{\varepsilon})(d_{\alpha}(x,K,y)+{\varepsilon})+{\varepsilon}.
$$
\end{proof}


\begin{proof}[Proof of Theorem \ref{HTK}]
We will show that the argument used in \cite[Theorem 1.2]{MR1484769}, for the case where $a$ is positive-definite and ${\beta}=0$, generalizes to the present context\footnote{%
The idea is to combine a standard argument for heat kernel upper bounds with a reflection trick. 
In terms of Markov processes, we give a random sign to each excursion of the diffusion process into $D$, viewing it as taking values in $D^-$ or $D^+$.
Then a generalization of the classical reflection principle for Brownian motion allows to express the density for paths from $x$ to $y$ via $K$ in terms of this enhanced process.
In fact the heat kernel $\tilde p$ for this process may be written in terms of $p$ and $p_D$, and we find it technically simpler to define $\tilde p$ in those terms, rather than set up the enhanced process.}
.
Consider the set $\tilde M=M^-\cup M^+$, where $M^\pm=K\cup D^\pm$ and $D^-,D^+$ are disjoint copies of $D=M\setminus K$. 
Write $\pi$ for the obvious projection $\tilde M\to M$.
For functions $f$ defined on $M$, we will write $f$ also for the function $f\circ\pi$ on $\tilde M$.
Thus we will sometimes consider $a$ as a quadratic form on $T^*D^\pm$ and ${\beta}$ as a $1$-form on $D^\pm$.
Define a measure $\tilde\nu$ on $\tilde M$ by
$$
\tilde\nu(A)=\nu(A\cap K)+\tfrac12\nu(\pi(A\cap D^-))+\tfrac12\nu(\pi(A\cap D^+)).
$$
Note that $\nu=\tilde\nu\circ\pi^{-1}$.
Now define
$$
\tilde p(t,x,y)=
\begin{cases}
p(t,x,y)+p_D(t,x,y),&\text{ if }x,y\in D^\pm,\\
p(t,x,y)-p_D(t,x,y),&\text{ if }x\in D^\pm\text{ and }y\in D^\mp,\\
p(t,x,y),&\text{ if }x\in K\text{ or }y\in K.
\end{cases}
$$
Given functions $f^-,f^+\inC^\infty _c(M)$ with $f^-=f^+$ on $K$, write $f$ for the function on $\tilde M$ such that $f=f^\pm\circ\pi$ on $M^\pm$,
and set $\bar f=(f^-+f^+)/2$ and $f^D=(f^+-f^-)/2$.
Let $\phi^\pm\inC^\infty _c(M)$ with $\phi^-=\phi^+$ on $K$ and define $\phi$ on $\tilde M$ and $\bar\phi$ and $\phi^D$ on $M$ similarly.
Define for $t\in(0,\infty)$ functions $u_t$ on $\tilde M$, $\bar u_t$ on $M$ and $u^D_t$ on $D$ by
$$
u_t(x)=\int_{\tilde M}\tilde p(t,x,y)f(y)\tilde\nu(dy)
$$
and
$$
\bar u_t(x)=\int_Mp(t,x,y)\bar f(y)\nu(dy),\quad u^D_t(x)=\int_Mp_D(t,x,y)f^D(y)\nu(dy).
$$
Then $\bar u_t$ and $u^D_t$ solve the heat equation with Dirichlet boundary conditions in $M$ and $D$ respectively.
It is straightforward to check that $u_t=u^\pm_t\circ\pi$ on $M^\pm$, where $u_t^\pm=\bar u_t\pm u_t^D$ and we extend $u_t^D$ by $0$ on $K$.
Hence
$$
\int_{\tilde M}\phi u_td\tilde\nu=\int_M\bar\phi\bar u_td\nu+\int_D\phi^Du_t^Dd\nu
$$
and so
\begin{align}
\notag
&\frac{d}{dt}\int_{\tilde M}\phi u_td\tilde\nu
=\frac{d}{dt}\int_M\bar\phi\bar u_td\nu+\frac{d}{dt}\int_D\phi^D u_t^Dd\nu\\
\notag
&=-\frac12\int_Ma(\nabla\bar\phi,\nabla\bar u_t)d\nu+\int_Ma(\bar\phi{\beta},\nabla\bar u_t)d\nu
-\frac12\int_Da(\nabla\phi^D,\nabla u_t^D)d\nu+\int_Da(\phi^D{\beta},\nabla u_t^D)d\nu\\
\label{WHE}
&=-\frac12\int_{\tilde M}a(\nabla\phi,\nabla u_t)d\tilde\nu+\int_{\tilde M}a(\phi{\beta},\nabla u_t)d\tilde\nu.
\end{align}

Fix $x,y\in D$ and $r_0\in(0,\infty)$ so that the balls $B(x,r_0)$ and $B(y,r_0)$ are compactly contained in $D$.
Write $x^-$ and $y^+$ for the unique points in $D^-$ and $D^+$ respectively such that $\pi(x^-)=x$ and $\pi(y^+)=y$.
Fix $r\in(0,r_0]$ and set
$$
B^-=\{z\in D^-:\pi(z)\in B(x,r)\},\quad B^+=\{z\in D^+:\pi(z)\in B(y,r)\}.
$$
Consider the set ${\cal F}_r$ of pairs of bounded locally Lipschitz functions $(w^-,w^+)$ on $M$ such that $w^-$ is constant on $B(x,r)$, $w^+$ is constant on $B(y,r)$, 
$w^-=w^+$ on $K$ and $a(\nabla w^\pm,\nabla w^\pm)\le1$ almost everywhere.
Set 
$$
d_r(x,K,y)=\sup\{w^+(y)-w^-(x):(w^-,w^+)\in{\cal F}_r\}.
$$
Fix $(w^-,w^+)\in{\cal F}_r$ and define a function $w$ on $\tilde M$ by setting $w=w^\pm\circ\pi$ on $M^\pm$.
Take $f^-=0$ and choose $f^+\ge0$ supported on $B(y,r)$ and such that $\int_Mf^2d\nu=2$.
Then $\int_{\tilde M}\tilde f^2d\tilde\nu=1$.
Fix ${\delta}\in[0,\infty)$ and ${\theta}\in{\mathbb R}$ and set $\psi={\theta} w$.
We deduce from \eqref{WHE} by a standard argument that
\begin{align*}
\frac d{dt}\int_{\tilde M}(e^\psi u_t)^2d\tilde\nu
&=-\int_{\tilde M}a(\nabla(e^{2\psi}u_t),\nabla u_t)d\tilde\nu+2\int_{\tilde M}a({\beta} e^{2\psi}u_t,\nabla u_t)d\tilde\nu\\
&=-\int_{\tilde M}a(\nabla u_t,\nabla u_t)e^{2\psi}d\tilde\nu+2\int_{\tilde M}a(({\beta}-\nabla\psi)u_t,\nabla u_t)e^{2\psi}d\tilde\nu\\
&\leqslant\int_{\tilde M}a({\beta}-\nabla\psi,{\beta}-\nabla\psi)(e^\psi u_t)^2d\tilde\nu\leqslant {\lambda}\int_{\tilde M}(e^\psi u_t)^2d\tilde\nu
\end{align*}
where ${\lambda}=\|a({\beta}-\nabla\psi,{\beta}-\nabla\psi)\|_\infty\leqslant(1+{\delta}){\theta}^2+(1+1/{\delta})\|a({\beta},{\beta})\|_\infty$. 
So, by Gronwall, we have
$$
\int_{\tilde M}(e^\psi u_t)^2d\tilde\nu\leqslant e^{{\lambda} t}\int_{\tilde M}(e^\psi\tilde f)^2d\tilde\nu=e^{{\lambda} t+2\psi(y^+)}
$$
and so
$$
\int_{B^-}u_t^2d\tilde\nu\leqslant e^{2{\theta}(w^+(y)-w^-(x))+{\lambda} t}.
$$
Moreover, since $u^-_t\ge0$ and $(\partial/\partial t)u_t^-={\cal L} u_t^-$ on $(0,\infty)\times D$, by the parabolic mean-value estimate, there is a constant $C(D,r_0)<\infty$ such that, 
for all $t\in(0,\infty)$ and all $r\in(0,r_0]$ with $r^2\leqslant t/2$, we have
$$
u_t(x^-)^2
\leqslant C\fint_{t-r^2}^t\fint_{B^-}u_s^2d\tilde\nu ds
\leqslant C\nu(B(x,r))^{-1}e^{2{\theta}(w^+(y)-w^-(x))+{\lambda} t}.
$$
Set $v_t(z)=p(t,x,K,z)$, then $(\partial/\partial t)v_t=\hat{\cal L} v_t$ on $(0,\infty)\times D$, where $\hat{\cal L} f=\frac12\operatorname{div}(a\nabla f)-a({\beta},\nabla f)$.
So, by the parabolic mean-value estimate again, we can choose $C(D,r_0)$ so that
$$
p(t,x,K,y)^2\leqslant 
C\fint_{t-r^2}^t\fint_{B(y,r)}p(s,x,K,z)^2\nu(dz)ds
=C\fint_{t-r^2}^t\fint_{B^+}\tilde p(s,x^-,z)^2\tilde\nu(dz)ds.
$$
Assume now that $r^2\leqslant t/4$.
For each $s\in[t-r^2,t]$, we may take $f^+=cp(s,x,K,.)1_{B(y,r)}$, where $c$ is chosen so that $\int_{\tilde M}f^2d\tilde\nu=1$.
For this choice of $f^+$, we have
$$
u_s(x^-)^2=\int_{B^+}\tilde p(s,x^-,z)^2\tilde\nu(dz).
$$
Hence
\begin{align*}
p(t,x,K,y)^2
&\leqslant C\nu(B(y,r))^{-1}\fint_{t-r^2}^tu_s(x^-)^2ds\\
&\leqslant C^2\nu(B(x,r))^{-1}\nu(B(y,r))^{-1}e^{2{\theta}(w^+(y)-w^-(x))+{\lambda} t}
\end{align*}
from which we obtain, on optimizing over ${\theta}$ and $w$,
$$
p(t,x,K,y)\leqslant C\nu(B(x,r))^{-1/2}\nu(B(y,r))^{-1/2}
\operatorname{exp}\left\{-\frac{d_r(x,K,y)^2}{2(1+{\delta})t}+\left(1+\frac 1{\delta}\right)\frac{\|a({\beta},{\beta})\|_\infty t}2\right\}.
$$
For $(w^-,w^+)\in{\cal F}_0$ and $r\in(0,r_0]$, set $\bar w_x=\sup_{z\in B(x,r)}w^-(z)$ and $\bar w_y=\inf_{z\in B(x,r)}w^+(z)$ and suppose that $\bar w_x\leqslant\bar w_y$.
We can define $(w_r^-,w_r^+)\in{\cal F}_r$ by setting $w_r^-(z)=\bar w_x\vee w(z)\wedge\bar w_y$ and $w_r^+(z)=\bar w_x\vee w(z)\wedge\bar w_y$.
Note that $|w^-(z)-w^-(x)|\leqslant r$ for all $z\in B(x,r)$ and $|w^+(z)-w^+(y)|\leqslant r$ for all $z\in B(y,r)$, so $w^+(y)-w^-(x)\leqslant w_r^+(y)-w_r^-(x)+2r$.
On optimizing over $w$, we see that
$$
d(x,K,y)\leqslant d_r(x,K,y)+2r.
$$
Hence, using the volume estimate \eqref{VE}, we finally obtain, for all $t\in(0,4r_0^2\vee1]$
$$
p(t,x,K,y)\leqslant Ct^{-d/(2{\alpha})}(1+d(x,K,y)^2/t)^{d/(2{\alpha})}
\operatorname{exp}\left\{-\frac{d(x,K,y)^2}{2(1+{\delta})t}+\left(1+\frac 1{\delta}\right)\frac{\|a({\beta},{\beta})\|_\infty t}2\right\}
$$
which is certainly sufficient to imply the claimed asymptotics.
\end{proof}

\section{Second variation of the energy in a sub-Riemannian manifold of non-constant rank}\label{SUB}
Let $M$ be a connected $C^\infty $ manifold of dimension $d$ and let $a$ be a $C^\infty $ non-negative quadratic form on $T^*M$ having a sub-Riemannian structure, as in Section \ref{NPD}. 
In the section and the next, we assume that $M$ is complete for the sub-Riemannian distance.
Since all questions which we address concern properties determined for a finite energy path by any neighbourhood of that path, this assumption of completeness results in no essential loss of generality.
Recall, for a continuous path ${\omega}$ in $M$, the notions of energy and driving path defined in Section \ref{NPD}.
In Section \ref{SRN}, we reviewed the notion of cut locus and defined the Gaussian measure $\mu_{\gamma}$ in terms of the bicharacteristic flow.
Extrapolating from the Riemannian case, we might hope to characterize these objects instead in terms of the energy function $I$.
This will be done in Section \ref{QG}.
In the present section, in preparation, we show that the set $H^{x,y}$ of finite-energy paths has, at suitably selected paths ${\omega}$, a well-defined set of tangent directions $T_\o H^{x,y}$, 
and that $T_\o H^{x,y}$ has the structure of a Hilbert space.
Then we show that, when ${\omega}$ is minimal and $\xi$ is a driving path for ${\omega}$, the energy has a well-defined second variation $Q_\xi$ in a dense set of tangent directions, 
which allows us to define a continuous non-negative quadratic form $Q_\xi$ on $T_\o H^{x,y}$.
Finally, we show that $Q_\xi$ is minimized over $\xi$ by a unique driving path ${\lambda}$, which is in fact a bicharacteristic.
Given a sub-Riemannian structure $X$ for $a$, there are related constructions on the set of control paths, which are well known, evidently depending on the choice of $X$.
We emphasise that $T_\o H^{x,y}$ and $Q=Q_{\lambda}$ depend on $a$ alone.

Two sub-Riemannian structures $(X_1,\dots,X_m)$ and $(Y_1,\dots,Y_n)$ are equivalent (see \cite{ABB}) if $X_\ell=\sum_{k=1}^nf_{\ell k}Y_k$
and $Y_k=\sum_{\ell=1}^mg_{k\ell}X_\ell$ for all $\ell$ and $k$, for some $C^\infty $ functions $f_{\ell k}$ and $g_{k\ell}$ on $M$.
If $a$ has constant rank, then all sub-Riemannian structures for $a$ are equivalent.
However, this is not true in general, as the following example shows.
In ${\mathbb R}^2$, take 
$$
X_1(x,y)=Y_1(x,y)=y\frac{\partial}{\partial x},\quad X_2(x,y)=Y_2(x,y)=\frac{\partial}{\partial y}
$$
\begin{equation}\label{SREX}
X_3(x,y)=\operatorname{sgn}(x)Y_3(x,y)=e^{-1/|x|}\frac{\partial}{\partial x}.
\end{equation}
Then 
$$
\sum_{\ell=1}^3X_\ell\otimes X_\ell=\sum_{\ell=1}^3Y_\ell\otimes Y_\ell
$$ 
but $X$ and $Y$ define inequivalent sub-Riemannian structures on ${\mathbb R}^2$.
Thus, in the non-constant rank case, 
we cannot establish that an object is intrinsic to $a$ by showing it is intrinsic to an equivalence class of sub-Riemannian structures.
Instead, our approach will be to work directly from $a$, using a sub-Riemannian structure for existence but not uniqueness.

Recall that $H^x$ denotes the set of finite-energy paths starting at $x$ and $H^{x,y}$ denotes the set of such paths terminating at $y$.
Also, $H^0(\R^m)$ denotes the space of absolutely continuous paths $h:[0,1]\rightarrow{\mathbb R}^m$ starting from $0$ such that
$$
\|h\|^2=\int_0^1| \dot h_t|^2dt<\infty.
$$
We fix a sub-Riemannian structure $X=(X_1,\dots,X_m)$ for $a$ and use this to construct some associated objects.
We will make clear which objects depend on the choice of $X$ and which do not.
Given ${\omega}\in H^x$, we define $h({\omega})\inH^0(\R^m)$ by
$$
\dot h_t({\omega})=X({\omega}_t)^*\xi_t
$$
where $\xi$ is any driving path for ${\omega}$.
Then $h({\omega})$ does not depend on the choice of $\xi$.
For $x\in M$ and $h\inH^0(\R^m)$, write $\phi(x,h)$ for the solution $(\phi_t)_{t\in[0,1]}$ of the differential equation
$$
\dot\phi_t=\sum_{\ell=1}^mX_\ell(\phi_t)\dot h_t^\ell,\quad \phi_0=x.
$$
Denote by $p_X(x)$ the orthogonal projection ${\mathbb R}^m\rightarrow(\ker X(x))^\perp$.
For $x\in M$ and $h,k\inH^0(\R^m)$, define $\pi=\pi(x,h)k\inH^0(\R^m)$ by
$$
\dot\pi_t=p_X(\phi_t(x,h))\dot k_t.
$$
The following parametrization of $H^x$ by $H^0(\R^m)$ using $X$ is well known.

\begin{proposition}\label{2.1}
Let $x\in M$ and let ${\gamma}\in H^x$. 
Then ${\gamma}=\phi(x,h({\gamma}))$.
Moreover, for all $h\inH^0(\R^m)$, if we set ${\omega}=\phi(x,h)$ and $\pi=\pi(x,h)h$, then ${\omega}\inH^x$ and $\pi=h({\omega})$ and $I({\omega})=\|\pi\|^2\leqslant\|h\|^2$.
\end{proposition}
\begin{proof}
Since $(\ker X(x))^\perp=\operatorname{im} X(x)^*$, there is a measurable map $\xi:[0,1]\rightarrowT^* M$ over ${\omega}$ such that $\dot\pi_t=X({\omega}_t)^*\xi_t$.
Then $\dot{\omega}_t=X({\omega}_t)X({\omega}_t)^*\xi_t=a({\omega}_t)\xi_t$ and
$$
I({\omega})=\int_0^1\<\xi_t,a({\omega}_t)\xi_t\>dt=\int_0^1 |X({\omega}_t)^*\xi_t|^2dt=\|\pi\|^2.
$$
We leave the remaining details to the reader.
\end{proof}

Let ${\omega}$ be a finite-energy path and fix a chart along ${\omega}$.
We say that a driving path $\xi$ for ${\omega}$ is {\em tame} if
$$
\int_0^1|\xi_t|^2dt<\infty
$$
where we have used the Euclidean norm in the chart.
This condition does not depend on the choice of chart.
We say that ${\omega}$ is {\em tame}%
\footnote{%
It is straightforward to see that, if the diffusivity $a$ has constant rank, then every finite-energy path is tame. 
On the other hand, consider on ${\mathbb R}^2$ the vector fields $X_1=\partial/\partial x_1,X_2=x_1\partial/\partial x_2$ and the path ${\omega}_t=(t,t^2/2)$ for $t\in[0,1]$. 
Since $X_1,X_2$ span ${\mathbb R}^2$ except on $\{x_1=0\}$, the only driving path for ${\omega}$ is $\xi_t=dx_1+(1/t)dx_2$. 
Then $\<\xi_t,a({\omega}_t)\xi_t\>=2$ for all $t$ so ${\omega}$ has finite energy, but $\xi$ is not tame.}
if it has a tame driving path.

Let ${\omega}$ be a tame path in $H^x$.
Fix a chart along ${\omega}$ and a tame driving path $\xi$ for ${\omega}$. 
We will define a space $T_\o H^x$ of {\em finite-energy variations} of ${\omega}$ which, for now, may appear to depend on the choice of chart and of $\xi$. 
Denote by $T_\o H^x$ the set of absolutely continuous maps $v:[0,1]\rightarrow TM$ over ${\omega}$ such that
\begin{equation}\label{v}
\dot v_t=(\nabla_{v_t}a)({\omega}_t)\xi_t+a({\omega}_t)\eta_t,\quad v_0=0
\end{equation}
for some measurable path $\eta:[0,1]\rightarrowT^* M$ over ${\omega}$, with
$$
\|v\|^2_\xi:=\int_0^1\<\eta_t,a({\omega}_t)\eta_t\>dt<\infty.
$$
Note that, in \eqref{v}, the meaning of the derivatives $\dot v_t$ and $\nabla_{v_t}a$ depends on the choice of chart.
Note also that, since $\xi$ is tame, by Gronwall's lemma, there is a constant $C<\infty$ such that $\|v\|_\infty\leqslant C\|v\|_\xi$ for all $v\inT_\o H^x$.
If $\eta$ can be chosen in \eqref{v} so that, in addition,
$$
\int_0^1|\eta_t|^2dt<\infty
$$
then we will call $v$ a {\em tame finite-energy variation}.
Write $\tilde T_\o H^x$ for the set of tame finite-energy variations of ${\omega}$. 
Set $T_\o H^{x,y}=\{v\inT_\o H^x:v_1=0\}$ and $\tilde T_\o H^{x,y}=\{v\in\tilde T_\o H^x:v_1=0\}$.

\begin{proposition}
Let ${\omega}\in H^x$ be tame, with tame driving path $\xi$.
Then $\tilde T_\o H^x$ is dense in $T_\o H^x$.
Moreover $\tilde T_\o H^{x,y}$ is dense in $T_\o H^{x,y}$.
\end{proposition}
\begin{proof}
For ${\varepsilon}>0$, set $\eta^{\varepsilon}_t=(\sqrt{a({\omega}_t)}+{\varepsilon} I)^{-1}\sqrt{a({\omega}_t)}\eta_t$ and define $v^{\varepsilon}$ by
$$
\dot v_t^{\varepsilon}=(\nabla_{v^{\varepsilon}_t}a)({\omega}_t)\xi_t+a({\omega}_t)\eta_t^{\varepsilon},\quad v_0^{\varepsilon}=0.
$$
Note that $\sqrt{a({\omega}_t)}(\eta^{\varepsilon}_t-\eta_t)\to0$ for all $t\in[0,1]$ and 
$$
\<\eta_t^{\varepsilon}-\eta_t,a({\omega}_t)(\eta^{\varepsilon}_t-\eta_t)\>={\varepsilon}^2|\eta_t^{\varepsilon}|^2\leqslant\<\eta_t^{\varepsilon},a({\omega}_t)\eta^{\varepsilon}_t\>+{\varepsilon}^2|\eta_t^{\varepsilon}|^2\leqslant\<\eta_t,a({\omega}_t)\eta_t\>.
$$
Hence
$$
\int_0^1\<\eta_t^{\varepsilon},a({\omega}_t)\eta_t^{\varepsilon}\>dt+{\varepsilon}^2\int_0^1|\eta_t^{\varepsilon}|^2dt\leqslant\int_0^1\<\eta_t,a({\omega}_t)\eta_t\>dt<\infty
$$
so $v^{\varepsilon}\in\tilde T_\o H^x$ for all ${\varepsilon}$, and by dominated convergence, as ${\varepsilon}\to0$,
$$
\|v^{\varepsilon}-v\|^2_\xi=\int_0^1\<\eta_t^{\varepsilon}-\eta_t,a({\omega}_t)(\eta^{\varepsilon}_t-\eta_t)\>dt\to0.
$$
We have proved the first assertion.

Set $V=\{v_1:v\inT_\o H^x\}$. 
There is a finite set $E\subseteqT_\o H^x$ such that $\{e_1:e\in E\}$ is a basis for $V$.
Given ${\delta}>0$, we can find for each $e\in E$ a tame $\tilde e\in\tilde T_\o H^x$ such that $\|e-\tilde e\|_\xi\leqslant{\delta}$.
Then $|e_1-\tilde e_1|\leqslant C{\delta}$.
We can therefore choose ${\delta}$ sufficiently small that $\{\tilde e_1:e\in E\}$ remains a basis for $V$.
Set
$$
A:=\sup\left\{\left\|\sum_{e\in E}{\alpha}_e\tilde e\right\|_\xi:{\alpha}_e\in{\mathbb R},\left|\sum_{e\in E}{\alpha}_e\tilde e_1\right|=1\right\}.
$$
Then $A<\infty$.
Given $v\inT_\o H^{x,y}$ and ${\varepsilon}>0$, there exists $v'\in\tilde T_\o H^x$ such that $\|v-v'\|_\xi<{\varepsilon}/(1+CA)$.
Then $v'_1\in V$ and $|v'_1|\leqslant C{\varepsilon}/(1+CA)$.
We can write $v'_1=\sum_{e\in E}{\alpha}_e\tilde e_1$ for some ${\alpha}_e\in{\mathbb R}$.
Set 
$$
\tilde v=v'-\sum_{e\in E}{\alpha}_e\tilde e
$$
then $\tilde v\in\tilde T_\o H^{x,y}$ and 
$$
\|v-\tilde v\|_\xi\leqslant\|v-v'\|_\xi+\left\|\sum_{e\in E}{\alpha}_e\tilde e\right\|_\xi\leqslant{\varepsilon}
$$
which proves the second assertion.
\end{proof}

By standard arguments, the map $\phi:M\timesH^0(\R^m)\rightarrow H$ is differentiable, in both arguments. 
Fix $x\in M$ and $h\inH^0(\R^m)$.
Set ${\omega}_t=\phi_t(x,h)$ and define
$$
u_t=\phi^*_t(x,h)=\frac{\partial}{\partial x}\phi_t(x,h)\in T_{{\omega}_t}M\otimesT^*_x M.
$$
For $k\inH^0(\R^m)$, define $v=v(k)$ by
$$
v_t=\frac{\partial}{\partial h}\phi_t(x,h)k\in T_{{\omega}_t}M.
$$
Then $u$ and $v$ satisfy the differential equations
\begin{align}\notag
\dot u_t&=\sum_{\ell=1}^m\nabla X_\ell({\omega}_t)u_t\dot h^\ell_t,\quad u_0=I,\\\label{e1}
\dot v_t&=\sum_{\ell=1}^m\nabla X_\ell({\omega}_t)v_t\dot h^\ell_t+X_\ell({\omega}_t)\dot k^\ell_t,\quad v_0=0
\end{align}
where the derivatives $\dot u_t$, $\dot v_t$ and $\nabla X_\ell$ are understood in the chosen chart.
Then, by the variation of constants formula,
\begin{equation}\label{vp}
v_t=\sum_{\ell=1}^mu_t\int_0^t u_s^{-1}X_\ell({\omega}_s)\dot k_s^\ell ds.
\end{equation}

\begin{proposition}\label{2.2}
Let ${\omega}\inH^x$ be tame and let $k\inH^0(\R^m)$.
Fix a chart along ${\omega}$ and a tame driving path $\xi$ for ${\omega}$, and write $T_\o H^x$ for the associated space of finite-energy variations.
Set
$$
v(k)=\frac\partial{\partial h}\phi(x,h({\omega}))k.
$$
Then $v(k)\inT_\o H^x$.
On the other hand, for any $v\inT_\o H^x$, we can define $k(\xi,v)\inH^0(\R^m)$ by
$$
\dot k_t^\ell(\xi,v)=\<\xi_t,\nabla X_\ell({\omega}_t)v_t\>+\<\eta_t,X_\ell({\omega}_t)\>
$$
where $\eta:[0,1]\to T^*M$ is a measurable path over ${\omega}$ satisfying {\rm\eqref{v}}.
Then $k(\xi,v)$ does not depend on the choice of $\eta$, and 
$$
v=\frac\partial{\partial h}\phi(x,h({\omega}))k(\xi,v).
$$
Moreover there is a constant $C<\infty$, depending only on $I({\omega})$, $\int_0^1|\xi_t|^2dt$ and a uniform bound for $X$ and $\nabla X$ along ${\omega}$, such that
$$
\|v(k)\|_\xi\leq C\|k\|,\quad \|k(\xi,v)\|\leqslant C\|v\|_\xi.
$$
Moreover, if $v=v(k)$, then $\pi(x,h({\omega}))(k(\xi,v)-k)=0$.
\end{proposition}
\begin{proof}
Write $h=h({\omega})$.
Then $v(k)$ satisfies \eqref{e1} so, by Gronwall's lemma, there is a constant $C<\infty$ such that $\|v(k)\|_\infty\leqslant C\|k\|$.
Here and below, $C$ is understood to have the dependence claimed in the statement.
Since $\xi$ is tame, we can define $g(k)\inH^0(\R^m)$ by
$$
\dot g_t^\ell(k)=\<\xi_t,\nabla X_\ell({\omega}_t)v_t(k)\>.
$$
There is a constant $C<\infty$ such that $\|g(k)\|\leq C\|v(k)\|_\infty$.
We can find a measurable map $\eta:[0,1]\rightarrowT^* M$ over ${\omega}$, such that
\begin{equation}\label{e3}
a({\omega}_t)\eta_t=\sum_{\ell=1}^mX_\ell({\omega}_t)(\dot k^\ell_t-\dot g^\ell_t(k)).
\end{equation}
Note that
$$
\nabla a(x)=\sum_{\ell=1}^m\nabla X_\ell(x)X_\ell(x)^*+X_\ell(x)\nabla X_\ell(x)^*
$$
so $v(k)$ satisfies 
\begin{equation*}
\dot v_t(k)=(\nabla_{v_t(k)}a)({\omega}_t)\xi_t+a({\omega}_t)\eta_t,\quad v_0(k)=0
\end{equation*}
and, moreover
$$
\|v(k)\|^2_\xi=\int_0^1\<\eta_t,a({\omega}_t)\eta_t\>dt=\|\pi(x,h)(k-g)\|^2\leq C\|k\|^2.
$$
Hence $v(k)\inT_\o H^x$ and $\|v(k)\|_\xi\leqslant C\|k\|$.

On the other hand, for $v\inT_\o H^x$, we see from \eqref{v} that there is a constant $C<\infty$ such that $\|v\|_\infty\leq C\|v\|_\xi$, 
so $\|k(\xi,v)\|\leq C\|v\|_\xi$. Moreover, we can write \eqref{v} in the form
$$
\dot v_t=\sum_{\ell=1}^m\nabla X_\ell({\omega}_t)v_t\dot h^\ell_t+X_\ell({\omega}_t)\dot k^\ell_t(\xi,v)
$$
so $v=v(k(\xi,v))$, as claimed.
Finally, if $v=v(k)$, then $v(k(\xi,v)-k)=0$, so, from \eqref{vp}, 
$$
\sum_{\ell=1}^mX_\ell({\omega}_t)(\dot k^\ell_t(\xi,v)-\dot k^\ell_t)=0
$$ 
for almost all $t$, and so $\pi(x,h)(k(\xi,v)-k)=0$.
\end{proof}

Fix ${\omega}$ and set $\tilde H=\{k\inH^0(\R^m):\pi(x,h({\omega}))k=k\}$. 
Then $\tilde H$ is a closed subspace of the Hilbert space $H^0(\R^m)$. 
Proposition \ref{2.2} shows that, for any choice of chart along ${\omega}$ and any choice of tame driving path $\xi$ for ${\omega}$,
the map $k\mapsto v(k)$ is a linear isomorphism $\tilde H\toT_\o H^x$, which is bounded with bounded inverse, when $T_\o H^x$ is given the norm $\|.\|_\xi$. 
Hence, the space $T_\o H^x$ does not depend on the choice of chart and driving path.
It clearly does not depend either on any choice of sub-Riemannian structure $X$. 
Moreover, the norms $\|.\|_\xi$ are all equivalent, and all make $T_\o H^x$ into a Hilbert space.

A (minimizing) finite-energy path ${\omega}\inH^{x,y}$ is said to be {\em regular} if the linear map
$$
\frac{\partial}{\partial h}\phi_1(x,h({\omega})):H^0(\R^m)\rightarrowT_y M
$$
is onto. 
By \eqref{vp}, this is equivalent to Bismut's condition that the {\em deterministic Malliavin covariance matrix}
\begin{eqnarray}
\label{mat}
C_1({\omega})
=\sum_{\ell=1}^m\int_0^1\left(u_t^{-1}X_\ell({\omega}_t)\right)\otimes\left(u_t^{-1}X_\ell({\omega}_t)\right)dt
=\int_0^1u_t^{-1}a({\omega}_t)(u_t^{-1})^*dt
\end{eqnarray}
is invertible.
In particular, when $a$ is positive-definite, every ${\omega}\inH^{x,y}$ is regular.
In general, these conditions may depend on the choice of sub-Riemannian structure.
By Proposition \ref{2.2}, a tame path ${\omega}\inH^{x,y}$ is regular if and only 
$$
\{v_1:v\inT_\o H^x\}=T_y M
$$
so for tame paths the notion of regularity depends only on $a$ and not on the choice of sub-Riemannian structure.
It is straightforward to see that ${\omega}$ is regular if and only if its time-reversal is regular. 
In \cite[pp. 22--24]{B}, Bismut gives an argument which shows that invertibility of $C_1({\omega})$ is intrinsic to the equivalence class of the sub-Riemannian structure $X$.
As the example \eqref{SREX} shows, this is not the same as being intrinsic to $a$. 

Consider a regular tame finite-energy path ${\omega}\inH^{x,y}$.
We use the sub-Riemannian structure $X$ to define $h({\omega})\inH^0(\R^m)$ as above and we
write $K$ for the kernel of the linear map $(\partial/\partial h)\phi_1(x,h({\omega})):H^0(\R^m)\to T_yM$.
Then
$$
K^\perp=\{k\inH^0(\R^m):\dot k_t^\ell=\<\eta_0,u_t^{-1}X_\ell({\omega}_t)\>,\eta_0\inT^*_x M\}
$$
and $\left.(\partial/\partial h)\phi_1(x,h({\omega}))\right|_{K^\perp}$ is invertible. 
By the implicit function theorem in Hilbert space, there exist ${\delta}>0$ and a $C^\infty $ map ${\theta}:K\to H$ such that, for all $k\in K$, we have
$$
\phi_1(x,h({\omega})+k+{\theta}(k))=y
$$
and such that, for all $k\in K$ and all $k'\in K^\perp$ with $\| k+k'\|<{\delta}$, we have ${\theta}(k)\in K^\perp$ and
$$
\phi_1(x,h+k+k')=y\quad\text{only if}\quad k'={\theta}(k).
$$
Note that ${\theta}(0)=0$. 
For $k\in K$ and ${\varepsilon}$ sufficiently small, we have
\begin{equation*}
\phi_1(x,h+{\varepsilon} k+{\theta}({\varepsilon} k))=y.
\end{equation*}
On differentiating in ${\varepsilon}$ at $0$, we obtain $(\partial/\partial h)\phi_1(x,h)(k+{\theta}'(0)k)=0$, so ${\theta}'(0)k\in K$.
Since ${\theta}$ takes values in $K^\perp$, we deduce that ${\theta}'(0)=0$.
On taking the second derivative, we obtain
\begin{equation}\label{kk}
\frac{\partial^2}{\partial h^2}\phi_1(x,h)(k,k)
+\frac{\partial}{\partial h}\phi_1(x,h){\theta}''(0)(k,k)=0.
\end{equation}
Since ${\theta}''(0)(k,k)\in K^\perp$, this equation determines ${\theta}''(0)$.

We note the following useful identity. 
Let $k\in K^\perp$ and let $k'\inH^0(\R^m)$. 
Then $\dot k_t^\ell=\<\eta_0,u_t^{-1}X_\ell({\omega}_t)\>$ for some $\eta_0\inT^*_x M$.
So
\begin{equation}\label{ui}
\<k,k'\>
=\int_0^1\<\eta_0,u_t^{-1}X_\ell({\omega}_t)\>(\dot k_t')^\ell\, dt=\<\eta_0,u_1^{-1}v_1(k')\>.
\end{equation}

The key arguments and computations, for a sub-Riemannian structure $X$, in the next two results are due to Bismut \cite[Theorem 1.17 and Theorem 1.24]{B}.
Our new contribution is to construct objects and show results which do not depend on the choice of $X$.
\begin{proposition}\label{2.3}
Let ${\omega}\inH^{x,y}$ be tame and regular and let $\xi$ be a tame driving path for ${\omega}$.
Let $v\in\tilde T_\o H^{x,y}$ be a tame finite-energy variation.
There exists a measurable map
$$
({\varepsilon},t)\mapsto\xi^{\varepsilon}_t:(-1,1)\times[0,1]\rightarrowT^* M
$$
such that $\xi^0=\xi$ and
\begin{itemize}
\item[{\rm (i)}] ${\omega}^{\varepsilon}=\pi\xi^{\varepsilon}\inH^{x,y}$ for all ${\varepsilon}$, with $\dot{\omega}_t^{\varepsilon}=a({\omega}_t^{\varepsilon})\xi_t^{\varepsilon}$,
\item[{\rm (ii)}] in any chart along ${\omega}$, there is a constant $C<\infty$ such that, for all ${\varepsilon}\in(-1,1)$,
$$
\sup_{t\in[0,1]}|{\omega}_t^{\varepsilon}-{\omega}_t-{\varepsilon} v_t|\leq C{\varepsilon}^2
$$
and, writing $\xi^{\varepsilon}=({\omega}^{\varepsilon},p^{\varepsilon})$ and $ \eta=({\omega},q)$, 
$$
\int_0^1\<p_t^{\varepsilon}-p_t-{\varepsilon} q_t,a({\omega}_t^{\varepsilon})(p_t^{\varepsilon}-p_t-{\varepsilon} q_t)\>dt\leq C{\varepsilon}^4.
$$
\end{itemize}
Moreover, for any such map ${\varepsilon}\mapsto \xi^{\varepsilon}$, the map ${\varepsilon}\mapsto I({\omega}^{\varepsilon})$ is differentiable at ${\varepsilon}=0$, and is twice differentiable at ${\varepsilon}=0$ if ${\omega}$ is minimal.
Define
$$
L(v)=\left.\frac{\partial}{\partial{\varepsilon}}\right|_{{\varepsilon}=0} I({\omega}^{\varepsilon})
$$
then $L$ extends uniquely to a continuous linear form on $T_\o H^{x,y}$.
In the case where ${\omega}$ is minimal, define
\begin{equation}\label{Q}
Q_\xi(v)=\frac{1}{2}\left.\frac{\partial^2}{\partial{\varepsilon}^2}\right|_{{\varepsilon}=0} I({\omega}^{\varepsilon})
\end{equation}
then $Q_\xi$ extends uniquely to a continuous quadratic form on $T_\o H^{x,y}$.
Finally, given a sub-Riemannian structure $X$ for $a$, for $h=h({\omega})$ and
$k=k(\xi,v)$, we have
$$
L(v)=2\<h,k\>
$$
and, when ${\omega}$ is minimal, then $h\in K^\perp$ and
$$
Q_\xi(v)=\|k\|^2+\<h,{\theta}''(0)(k,k)\>.
$$
\end{proposition}
\begin{proof}
Fix a sub-Riemannian structure $X$ for $a$.
We will use $X$ to construct a map $\xi^{\varepsilon}$ having the claimed properties. 
Set $h=h({\omega})$ and $k=k(\xi,v)$.
Since $v_1=0$, we have $k\in K$.
For ${\varepsilon}\in(-1,1)$, set $h^{\varepsilon}=h+{\varepsilon} k+{\theta}({\varepsilon} k)$ and set ${\omega}^{\varepsilon}=\phi(x,h^{\varepsilon})$.
Then ${\omega}^{\varepsilon}\inH^{x,y}$ for all ${\varepsilon}$.
Moreover, the map ${\varepsilon}\mapsto{\omega}^{\varepsilon}:(-1,1)\to{\Omega}$ is $C^1$ and there is a constant $C<\infty$ such that
$$
\sup_{t\in[0,1]}|{\omega}_t^{\varepsilon}-{\omega}_t-{\varepsilon} v_t|\leq C{\varepsilon}^2.
$$
Moreover, fixing a chart along ${\omega}$ and writing $\xi=({\omega},p)$ and $\eta=({\omega},q)$, we have 
$$
\dot{\omega}_t^{\varepsilon}-a({\omega}_t^{\varepsilon})(p_t+{\varepsilon} q_t)=X_\ell({\omega}_t^{\varepsilon})\dot g_t^{{\varepsilon},\ell}
$$
where
\begin{align*}
\dot g_t^{{\varepsilon},\ell}&=\dot h_t^\ell+{\varepsilon} \dot k_t^\ell+\dot{\theta}_t^\ell({\varepsilon} k)
-\<p_t+{\varepsilon} q_t,X_\ell({\omega}_t^{\varepsilon})\>\\
&=-\left\<p_t,X_\ell({\omega}_t^{\varepsilon})-X_\ell({\omega}_t)-{\varepsilon}\nabla X_\ell({\omega}_t)v_t\right\>
-{\varepsilon}\<q_t,X_\ell({\omega}_t^{\varepsilon})-X_\ell({\omega}_t)\>+\dot{\theta}_t^\ell({\varepsilon} k).
\end{align*}
Since $\xi$ and $\eta$ are tame, there is a constant $C<\infty$ such that, for all ${\varepsilon}$, we have
$$
\int_0^1 |\dot g_t^{\varepsilon}|^2dt\leq C{\varepsilon}^4
$$
so we can find a measurable map $r_t^{\varepsilon}$ such that
$$
X_\ell({\omega}_t^{\varepsilon})\dot g_t^{{\varepsilon},\ell}=a({\omega}_t^{\varepsilon})r_t^{\varepsilon}
$$
and
$$
\int_0^1\<r_t^{\varepsilon},a({\omega}_t^{\varepsilon})r_t^{\varepsilon}\>dt\leq C{\varepsilon}^4.
$$
If we now set $p_t^{\varepsilon}=p_t+{\varepsilon} q_t+r_t^{\varepsilon}$, then $\xi_t^{\varepsilon}=({\omega}_t^{\varepsilon},p_t^{\varepsilon})$
has the required properties.

Suppose now, more generally, that $\xi^{\varepsilon}=({\omega}^{\varepsilon},p^{\varepsilon})$ and $\eta=({\omega},q)$ are maps having the properties described in the statement. 
Define $h^{\varepsilon}\inH^0(\R^m)$ by
$$
\dot h_t^{{\varepsilon},\ell}=\<p_t^{\varepsilon},X_\ell({\omega}_t^{\varepsilon})\>.
$$
Then
\begin{align*}
&\dot h_t^{{\varepsilon},\ell}-\dot h_t^\ell-{\varepsilon} \dot k_t^\ell\\
&\quad=\<p_t^{\varepsilon}-p_t-{\varepsilon} q_t,X_\ell({\omega}_t^{\varepsilon})\>+{\varepsilon}\<q_t,X_\ell({\omega}_t^{\varepsilon})-X_\ell({\omega}_t)\>+\<p_t,X_\ell({\omega}_t^{\varepsilon})-X_\ell({\omega}_t)-{\varepsilon}\nabla X_\ell({\omega}_t)v_t\>
\end{align*}
so, since $\xi$ and $\eta$ are tame, there is a constant $C<\infty$ such that, for all ${\varepsilon}$, we have
$$
\|h^{\varepsilon}-h-{\varepsilon} k\|\leq C{\varepsilon}^2.
$$
Write $k^{\varepsilon}$ for the orthogonal projection of ${\varepsilon}^{-1}(h^{\varepsilon}-h)$ onto $K$. 
Since $k\in K$, we have $\|k^{\varepsilon}-k\|\leqslant\|{\varepsilon}^{-1}(h^{\varepsilon}-h)-k\|\leqslant C{\varepsilon}$. 
Since $\phi_1(x,h^{\varepsilon})=y$, we have, for ${\varepsilon}$ sufficiently small,
$$
h^{\varepsilon}=h+{\varepsilon} k^{\varepsilon}+{\theta}({\varepsilon} k^{\varepsilon}).
$$
Hence, as ${\varepsilon}\rightarrow 0$,
$$
I({\omega}^{\varepsilon})-I({\omega})=\|h+{\varepsilon} k^{\varepsilon}+{\theta}({\varepsilon} k^{\varepsilon})\|^2-\|h\|^2
=2{\varepsilon}\<h,k\>+O({\varepsilon}^2).
$$
Hence ${\varepsilon}\mapsto I({\omega}^{\varepsilon})$ is differentiable at ${\varepsilon}=0$ with derivative $L(v)=2\<h,k\>$.
By Proposition \ref{2.2}, we have $\|k(\xi,v)\|\leqslant C\|v\|_\xi$, so $L$ is continuous on $\tilde T_\o H^{x,y}$, which is dense in $T_\o H^{x,y}$, so $L$ extends uniquely to $T_\o H^{x,y}$.

Now, if ${\omega}$ is minimal, we must have $L(v)=0$ for all $v\in\tilde T_\o H^{x,y}$ and hence for all $v\inT_\o H^{x,y}$.
So, for any $k'\in K$, we have
$$
\<h,k'\>=\<h,k(\xi,v(k'))\>=L(v(k'))=0
$$
and so $h\in K^\perp$.
Then, for ${\omega}^{\varepsilon}$ as above, as ${\varepsilon}\rightarrow 0$, we have
$$
I({\omega}^{\varepsilon})-I({\omega})={\varepsilon}^2\|k^{\varepsilon}\|^2+2\<h,{\theta}({\varepsilon} k^{\varepsilon})\>+\|{\theta}({\varepsilon} k^{\varepsilon})\|^2={\varepsilon}^2\{\|k\|^2+\<h,{\theta}''(0)(k,k)\>\}+O({\varepsilon}^3).
$$
This shows that ${\varepsilon}\mapsto I({\omega}^{\varepsilon})$ is twice differentiable at ${\varepsilon}=0$ with the claimed second derivative, which is then continuous on $\tilde T_\o H^{x,y}$
and extends uniquely to $T_\o H^{x,y}$.
\end{proof}

Note that, since $L(v)$ and $Q_\xi(v)$ can be computed either by choosing a suitable family of paths ${\varepsilon}\mapsto\xi^{\varepsilon}$ or by choosing a sub-Riemannian structure $X$, they depend on neither choice.
For ${\gamma}\inH^{x,y}$ tame, regular and minimal, we define a continuous quadratic form $Q$ on ${T_\g H^{x,y}}$ by
$$
Q(v)=\inf_{\xi}Q_\xi(v)
$$
where the infimum is taken over all tame driving paths $\xi$.

\begin{proposition}\label{2.4}
Let ${\gamma}\in H^{x,y}$ be tame, regular and minimal. 
Then there exists a unique tame driving path ${\lambda}$ such that $Q=Q_{\lambda}$. 
The path ${\lambda}$ is a bicharacteristic and is the only bicharacteristic which is a driving path for ${\gamma}$.
Moreover, given a sub-Riemannian structure $X$ for $a$, we have
$$
Q_\xi(v)=Q(v)+\|k(\xi,v)-k({\lambda},v)\|^2,\quad  Q(v)=q(k({\lambda},v))
$$
where $q$ is the quadratic form on $K=\ker(\partial/\partial h)\phi_1(x,h({\gamma}))$ given by
$$
q(k)=\|k\|^2-\left\<{\lambda}_1,\frac{\partial^2}{\partial h^2}\phi_1(x,h({\gamma}))(k,k)\right\>.
$$
\end{proposition}
\begin{proof}
Choose a sub-Riemannian structure $X$ for $a$ and define as above
$$
h=h({\gamma}),\quad u_t=\frac\partial{\partial x}\phi_t(x,h),\quad K=\ker\frac\partial{\partial h}\phi_1(x,h).
$$
Then, by Proposition \ref{2.3}, we have $h\in K^\perp$, so there exists a unique ${\lambda}_0\in T^*_xM$ such that
$$
\dot h_t^\ell=\<{\lambda}_t,X_\ell({\gamma}_t)\>
$$
for almost all $t$, where ${\lambda}_t=(u_t^{-1})^*{\lambda}_0$. 
Fix a chart along ${\gamma}$ and write ${\lambda}_t=({\gamma}_t,p_t)$.
Then
$$
\dot{\gamma}_t=\sum_{\ell=1}^mX_\ell({\gamma}_t)\dot h_t^\ell,\quad \dot p_t=-\sum_{\ell=1}^m\<p_t,\nabla X_\ell({\gamma}_t)\>\dot h_t^\ell
$$
so ${\lambda}$ is a bicharacteristic.
Suppose, on the other hand, that $\xi=({\gamma},q)$ is a bicharacteristic over ${\gamma}$ and define $k\inH^0(\R^m)$ by $\dot k_t^\ell=\<\xi_t,X_\ell({\gamma}_t)\>$.
Then ${\gamma}=\phi(x,k)$ and $I({\gamma})=\|k\|^2$ and
$$
\dot q_t=-\sum_{\ell=1}^m\<q_t,\nabla X_\ell({\gamma}_t)\>\dot k_t^\ell
$$
so $k=h({\gamma})$ and $q=p$, proving uniqueness.

Write $h$ for $h({\gamma})$ and recall that $\dot h_t^\ell=\<{\lambda}_0,u_t^{-1}X_\ell({\gamma}_t)\>$ and ${\lambda}_1=(u_1^{-1})^*{\lambda}_0$.
Then, by \eqref{kk} and \eqref{ui}, we have
\begin{equation}\label{ix}
q(k)=\|k\|^2+\<h,{\theta}''(0)(k,k)\>.
\end{equation}
So $q(k(\xi,v))=Q_\xi(v)$ by Proposition \ref{2.3}.
Take now $v=v(k)$ and set $k'=k-k({\lambda},v)$.
Then $(\partial/\partial h)\phi(x,h)k'=0$, so $\sum_{\ell=1}^mX_\ell({\gamma}_t)(\dot k_t')^\ell=0$ for almost all $t$, and so $\phi(x,h+{\varepsilon} k')={\gamma}$ for all ${\varepsilon}\in{\mathbb R}$. 
Hence $q(k')=\|k'\|^2$ and it will suffice to show that $q(k({\lambda},v),k')=0$.
Recall that $\dot k_t^\ell({\lambda},v)=\<{\lambda}_t,\nabla X_\ell({\gamma}_t)v_t\>+\<\eta_t,X_\ell({\gamma}_t)\>$.
We differentiate the identity
\begin{equation}
\label{k}
\left(\frac{\partial}{\partial x}\phi_1(x,h)\right)^{-1}
\frac{\partial}{\partial h}\phi_1(x,h)k'
=\sum_{\ell=1}^m\int_0^1\left(\frac{\partial}{\partial x}\phi_t(x,h)\right)^{-1}X_\ell(\phi_t(x,h))(\dot k_t')^\ell\, dt
\end{equation}
in $h$, in the direction $k({\lambda},v)$, to obtain
$$
u_1^{-1}\frac{\partial^2}{\partial h^2}\phi_1(x,h)(k({\lambda},v),k')
=\sum_{\ell=1}^m\int_0^1 u_t^{-1}\nabla X_\ell({\gamma}_t)v_t(\dot k_t')^\ell\, dt.
$$
Hence
$$
\left\<{\lambda}_1,\frac{\partial^2}{\partial h^2}\phi_1(x,h)(k({\lambda},v),k')\right\>
=\sum_{\ell=1}^m\int_0^1\<{\lambda}_t,\nabla X_\ell({\gamma}_t)v_t\>(\dot k_t')^\ell\, dt=\<k({\lambda},v),k'\>.
$$
This shows that $q(k({\lambda},v),k')=0$ as required.
\end{proof}

Note that our extra condition of tameness, under which regularity of finite-energy paths is intrinsic to $a$, does not exclude any normal minimal paths as these all have smooth driving paths.

\section{Characterization of the cut locus and fluctuation measure by the energy function}\label{QG}
Throughout this section, ${\gamma}$ is a regular tame minimal path in $H^{x,y}$, ${\lambda}$ is the unique bicharacteristic projecting to ${\gamma}$ and $Q$ is the quadratic form on ${T_\g H^{x,y}}$ defined in the preceding section. 
Fix $s\in[0,1]$ and $\eta_s\inT^*_{\g_s}M$. 
For $t\in[s,1]$, define $J_{ts}:T^*_{\g_s}M\rightarrowT_{\g_t}M$ by
$$
J_{ts}\eta_s
=\left.\frac{\partial}{\partial{\varepsilon}}\right|_{{\varepsilon}=0}\pi\psi_{t-s}({\lambda}_s+{\varepsilon} \eta_s)
$$
where $(\psi_t)_{t\in{\mathbb R}}$ is the bicharacteristic flow.
We already defined $J_t=J_{t0}$ by equation \eqref{JFE} above.
Set $v_t=J_{ts}\eta_s$.
Given a sub-Riemannian structure $X$ for $a$, the following equations determine $\phi_{ts}(x,h)$ and $\dot h_t$, for $t\in[s,1]$ and $x$ near ${\gamma}_t$, 
as functions of ${\lambda}_s\inT^*_{\g_s}M$
$$
\phi_{ss}(x,h)=x,\quad
\dot\phi_{ts}(x,h)=\sum_{\ell=1}^mX_\ell(\phi_{ts}(x,h))\dot h_t^\ell,\quad
\dot h_t^\ell=\left\<{\lambda}_s,\left(\frac\partial{\partial x}\phi_{ts}({\gamma}_s,h)\right)^{-1}X_\ell(\phi_{ts}({\gamma}_s,h))\right\>.
$$
On differentiating in ${\lambda}_s$ in the direction $\eta_s$, we find that $v_t=(\partial/\partial h)\phi_{ts}({\gamma}_s,h)k$, where $k$ satisfies, for $t\in[s,1]$,
\begin{equation}\label{ak}
\dot k_t^\ell =\<\eta_s,u_{ts}^{-1}X_\ell({\gamma}_t)\>+A_{ts}^\ell({\lambda},k).
\end{equation}
Here, $u_{ts}=(\partial/\partial x)\phi_{ts}({\gamma}_s,h)=u_tu_s^{-1}$ and
$$
A_{ts}^\ell({\lambda},k)
=\left\<{\lambda}_s,\frac{\partial}{\partial h}\left[\left(\frac{\partial}{\partial x}\phi_{ts}({\gamma}_s,h) \right)^{-1}X_\ell(\phi_{ts}({\gamma}_s,h))\right]k\right\>.
$$
By writing a differential equation for $A_{ts}^\ell({\lambda},k)$, we find a constant $C<\infty$ such that, for all $t\in[s,1]$,
$$
|A_{ts}({\lambda},k)|^2\leq C\int_s^t|\dot k_r|^2\, dr.
$$
It follows that the equation \eqref{ak} uniquely determines $k$.
Note that, if we set $v_t=0$ and $\dot k_t=0$ for $t\in[0,s)$, then we have, for all $t\in[0,1]$,
\begin{equation}\label{al}
v_t=\frac\partial{\partial h}\phi_t(x,h)k, \quad 
\dot k_t^\ell =\<u_s^*\eta_s,u_t^{-1}X_\ell({\gamma}_t)\>1_{\{t\geqslant s\}}+A_{t0}^\ell({\lambda},k).
\end{equation}

\begin{proposition}\label{2.5}
Let ${\gamma}\inH^{x,y}$ be tame, regular and minimal and let $v\in{T_\g H^{x,y}}$.
Then $Q(v)=0$ if and only if $v_t=J_t\eta_0$ for all $t\in[0,1]$, for some $\eta_0\inT^*_x M$.
\end{proposition}
\begin{proof}
Choose a sub-Riemannian structure $X$ for $a$.
Suppose that $Q(v)=0$.
Set $k=k({\lambda},v)$, then $q(k)=0$. 
On differentiating the identity \eqref{k} in $h$, in the direction $k$, we obtain
$$
u_1^{-1}\frac{\partial^2}{\partial h^2}\phi_1(x,h)(k,k')
=\sum_{\ell=1}^m\int_0^1 (\dot k_t')^\ell\frac{\partial}{\partial h}\left[\left(\frac{\partial}{\partial x}\phi_t(x,h) \right)^{-1}X_\ell(\phi_t(x,h))\right]k\, dt
$$
so, by Proposition \ref{2.4},
\begin{equation}\label{DIFF}
q(k,k')=\sum_{\ell=1}^m\int_0^1 (\dot k_t')^\ell\left\{\dot k_t^\ell
-\left\<{\lambda}_0,\frac{\partial}{\partial h}\left[\left(\frac{\partial}{\partial x}\phi_t(x,h) \right)^{-1}X_\ell(\phi_t(x,h))\right]k\right\>\right\}\, dt.
\end{equation}
Since ${\gamma}$ is minimal, $q$ is non-negative on $K$, so $q(k,k')=0$ for all $k'\in K$.
Hence there exists an $\eta_0\inT^*_x M$ such that
$$
\dot k_t^\ell=\<\eta_0,u_t^{-1}X_\ell({\gamma}_t)\>+A_{t0}^\ell({\lambda},k).
$$
As we argued above, this forces $v_t=J_t\eta_0$ for all $t$.
On the other hand, if $v_t=J_t\eta_0$ then, by the same calculations, we see that $v=(\partial/\partial h)\phi(x,h)k$ where $k$ satisfies $q(k,k')=0$ for all $k'\in K$.
But $k\in K$, so $Q(v)\leqslant q(k,k)=0$.
\end{proof}

\begin{proposition}\label{2.6}
Let ${\gamma}\inH^{x,y}$ be tame, regular and minimal and let $\eta_0\in T^*_x M$.
Then $J_t\eta_0=0$ for all $t$ only if $\eta_0=0$.
\end{proposition}
\begin{proof}
Fix a sub-Riemannian structure $X$ for $a$ and a chart along ${\gamma}$.
Write $\psi_t({\lambda}_0+{\varepsilon}\eta_0)=({\gamma}^{\varepsilon}_t,p^{\varepsilon}_t)$ and $\dot h^{{\varepsilon},\ell}_t=\<p_t^{\varepsilon},X_\ell({\gamma}^{\varepsilon}_t)\>$.
Set 
$$
v_t=J_t\eta_0=\left.\frac\partial{\partial{\varepsilon}}\right|_{{\varepsilon}=0}{\gamma}^{\varepsilon}_t,\quad
r_t=\left.\frac\partial{\partial{\varepsilon}}\right|_{{\varepsilon}=0}p^{\varepsilon}_t,\quad
\dot k^\ell_t=\left.\frac\partial{\partial{\varepsilon}}\right|_{{\varepsilon}=0}\dot h^{{\varepsilon},\ell}_t.
$$
By differentiating the bicharacteristic equations we obtain, for all $t\in[0,1]$,
$$
\dot v_t=\sum_{\ell=1}^m\nabla X_\ell({\gamma}_t)v_t\dot h_t^\ell+X_\ell({\gamma}_t)\dot k_t^\ell,\quad
\dot k_t^\ell=\<r_t,X_\ell({\gamma}_t)\>+\<p_t,\nabla X_\ell({\gamma}_t)v_t\>.
$$
But $v_t=J_t\eta_0=0$ for all $t$ so $\sum_{\ell=1}^mX_\ell({\gamma}_t)\dot k_t^\ell=0$ and so $|\dot k_t|^2=\sum_{\ell=1}^m\<r_t,X_\ell({\gamma}_t)\>^2=0$ for all $t$.
Now 
\begin{equation}\label{ELL}
\dot h_t^{{\varepsilon},\ell}=\<p_t^{\varepsilon},X_\ell({\gamma}_t^{\varepsilon})\>=\<{\lambda}_0+{\varepsilon}\eta_0,(u_t^{\varepsilon})^{-1}X_\ell({\gamma}_t^{\varepsilon})\>
\end{equation}
where 
$$
\dot u_t^{\varepsilon}=\nabla X_\ell({\gamma}_t^{\varepsilon})u_t^{\varepsilon}\dot h^{{\varepsilon},\ell}_t,\quad u_0=\operatorname{id}. 
$$
By differentiating this equation, we see that $(\partial/\partial{\varepsilon})|_{{\varepsilon}=0}u_t^{\varepsilon}=0$, so on differentiating \eqref{ELL} we obtain
$$
\<\eta_0,u_t^{-1}X_\ell({\gamma}_t)\>=0
$$
for all $t$, and this implies that $\eta_0=0$ since ${\gamma}$ is regular.
\end{proof}

The next result is a characterization of the cut locus by path properties intrinsic to $a$.
An analogous characterization in terms of a given sub-Riemannian structure is given in \cite[Th\'eor\`eme 1.18]{BA}. 
We note the use of Proposition \ref{2.6}. 
\begin{theorem}\label{2.7}
Let $x,y\in M$. The following are equivalent
\begin{itemize}
\item[{\rm (i)}] $(x,y)$ lies outside the cut locus,
\item[{\rm (ii)}] there is a unique minimal path ${\gamma}\inH^{x,y}$, which is tame and regular, and the quadratic form $Q$ is positive-definite on ${T_\g H^{x,y}}$.
\end{itemize}
\end{theorem}
\begin{proof}
In both (i) and (ii) we have a unique minimal path ${\gamma}\inH^{x,y}$.
By Proposition \ref{2.3}, if ${\gamma}$ is tame and regular then ${\gamma}$ is the projection of a bicharacteristic, while any projected bicharacteristic is tame.

Suppose that $Q$ is positive-definite and that $J_1\eta_0=0$ for some $\eta_0\inT^*_x M$.
Then $v_t=J_t\eta_0\in{T_\g H^{x,y}}$ and $Q(v)=0$, so $J_t\eta_0=0$ for all $t$, and so $\eta_0=0$, by Proposition \ref{2.6}.
Hence $J_1$ is invertible.

Suppose, on the other hand, that $J_1$ is invertible. 
Since $\{v_1:v\in{T_\g H^{x,y}}\}$ contains $\{J_1\eta_0:\eta_0\inT^*_x M\}$, we see that ${\gamma}$ is regular.
Suppose further that $Q(v)=0$ for some $v\in{T_\g H^{x,y}}$.
By Proposition \ref{2.5}, we have $v_t=J_t\eta_0$ for some $\eta_0\inT^*_x M$.
Then $J_1\eta_0=0$, so $\eta_0=0$, and so $v=0$.
Hence $Q$ is positive-definite.
\end{proof}

Recall from \eqref{JFE} that we define $K_t:T^*_y M\rightarrowT_{\g_t}M$ by
$$
K_t\eta_1
=\left.\frac{\partial}{\partial{\varepsilon}}\right|_{{\varepsilon}=0}\pi\psi_{-(1-t)}({\lambda}_1+{\varepsilon}\eta_1)
$$
and, for $s,t\in[0,1]$ with $s\leqslant t$, set
$$
C(s,t)=C(t,s)^*=J_sJ_1^{-1}K_t^*\inT_{\g_s}M\otimesT_{\g_t}M.
$$

\begin{proposition}\label{2.8}
Let ${\gamma}\inH^{x,y}$ be tame, regular and minimal.
Suppose that $x$ and $y$ are non-conjugate along ${\gamma}$.
Let $s\in[0,1]$ and ${\beta}\inT^*_{\g_s}M$.
For $t\in[0,1]$, set $v_t^{{\beta},s}=C(t,s){\beta}$.
Then $v^{{\beta},s}\in{T_\g H^{x,y}}$ and, for all $v\in{T_\g H^{x,y}}$,
$$
Q(v,v^{{\beta},s})=\<{\beta},v_s\>.
$$
\end{proposition}
\begin{proof}
Set $\eta_0=J_1^{-1}J_{1s}{\beta}$.
We can define $w=w^1-w^2\in{T_\g H^{x,y}}$ by setting $w^1_t=J_t\eta_0$ for all $t$, and $w^2_t=0$ for $t\leqslant s$ and $w^2_t=J_{ts}{\beta}$ for $t>s$.
Then $w=(\partial/\partial h)\phi(x,h)k$, where $k=k^1-k^2$ with
$$
(\dot k^1_t)^\ell=\<\eta_0,u_t^{-1}X_\ell({\gamma}_t)\>+A_{0t}^\ell({\lambda},k^1),\quad
(\dot k^2_t)^\ell=\<u^*_s{\beta},u^{-1}_tX_\ell({\gamma}_t)\>1_{\{t\geqslant s\}}+A_{0t}^\ell({\lambda},k^2).
$$
Take $v'\in{T_\g H^{x,y}}$ and set $k'=k({\lambda},v')$.
Then, using \eqref{DIFF} for the second equality, we have
\begin{align*}
Q(w,v')
&=q(k,k')=\sum_{\ell=1}^m\int_0^1(\dot k_t')^\ell\left\{\dot k_t^\ell-A_{0t}^\ell({\lambda},k) \right\}\, dt\\
&=\sum_{\ell=1}^m\int_0^1\<\eta_0,u_t^{-1}X_\ell({\gamma}_t)\>(\dot k_t')^\ell\, dt-\int_s^1\<u^*_s{\beta},u_t^{-1}X_\ell({\gamma}_t)\>(\dot k_t')^\ell\, dt=\<{\beta},v_s'\>.
\end{align*}
For $t\leq s$ we have
$$
w_t=J_{t0}J_{10}^{-1}J_{1s}{\beta}=J_tJ_1^{-1}K_s^*{\beta}=C(t,s){\beta}.
$$
Consider now the analogous construction in reverse time.
Set $\eta_1=K_0^{-1}K_{0s}{\beta}$ and define $\hat w=\hat w^1-\hat w^2\in{T_\g H^{x,y}}$ by setting $\hat w_t^1=K_t\eta_1$,
and $\hat w^2_t=0$ for $t\geqslant s$ and $\hat w^2_t=K_{ts}{\beta}$ for $t<s$.
Then, by the same argument, $Q(\hat w,v')=\<{\beta},v'_s\>$ for all $v'\in{T_\g H^{x,y}}$.
But $w-\hat w\in{T_\g H^{x,y}}$, so this implies $Q(w-\hat w,w-\hat w)=0$ and hence $\hat w=w$.
So, for $t\geqslant s$, we have
$$
w_t=\hat w_t=K_{t1}K_{01}^{-1}K_{0s}{\beta}=K_t(J_1^{-1})^*J_s^*{\beta}=C(s,t)^*{\beta}.
$$
Hence $w=v^{{\beta},s}$ and $v^{{\beta},s}$ has the claimed property.
\end{proof}

\begin{theorem}\label{SRMR}
Let $x,y\in M$ and suppose that $(x,y)$ lies outside the cut locus.
Denote the unique minimal path in $H^{x,y}$ by ${\gamma}$.
Recall the definitions \eqref{BWW} of the process $W$ and \eqref{SDEF} of the random variables $S(z)$ for $z\in{\mathbb R}^d$.
Then ${\mathbb E}(e^{pS(z)/2})<\infty$ for all $z\in{\mathbb R}^d$, for some $p>1$.

Define a new probability measure $\tilde\mathbb{P}$ on ${\O^0(\R^m)}$ by $d\tilde\mathbb{P}/d\mathbb{P}\propto e^{S(0)/2}$ 
and write $\mu_\g$ for the law on ${T_\g\oxy}$ of $Y(0)=v(W)$ under $\tilde\mathbb{P}$.
Then $\mu_\g$ is a zero-mean Gaussian probability measure on ${T_\g\oxy}$ with the following properties
\begin{itemize}
\item[{\rm (i)}] for all continuous linear functionals $\phi$ on ${T_\g\oxy}$, 
$$
\int_{T_\g\oxy}\phi(v)^2\mu_\g(dv)=Q(\tilde\phi)
$$
where $\tilde\phi\in{T_\g H^{x,y}}$ is given by $\phi(v)=Q(\tilde\phi,v)$ for all $v\in{T_\g H^{x,y}}$,
\item[{\rm (ii)}] for all $s,t\in[0,1]$ with $s\leqslant t$,
$$
\int_{{T_\g\oxy}} v_s\otimes v_t\,\mu_\g(dv)=J_sJ_1^{-1}K_t^*.
$$
\end{itemize}
Moreover, the properties {\rm (i)} and {\rm (ii)} both characterize $\mu_\g$ uniquely.
\end{theorem}
\begin{proof}
Choose a sub-Riemannian structure $X$ for $a$ and consider the continuous quadratic form on $H^0(\R^m)$ given by
$$
s(k)=\<{\lambda}_1,(\partial/\partial h)^2\phi_1(x,h)(k,k)\>.
$$
Fix an orthonormal basis $(e_n:n\in{\mathbb N})$ for $K$ which diagonalizes $s$ on $K$.
Formally, we have $S(z)=s(W(z))$ so, by a standard calculus for the Wiener chaos, 
\begin{equation}\label{TSZ}
S(z)=\sum_ns(e_n)\<W,e_n\>^2+2\sum_ns(e_n,{\tau}(z))\<W,e_n\>+s({\tau}(z))
\end{equation}
where the sums are understood in $L^2(\mathbb{P})$ and the random variables $\<W,e_n\>$ are independent standard Gaussians.
Set $s_n=-s(e_n)$ and ${\sigma}_n=s(e_n,{\tau}(z))$.
We see from \eqref{TSZ} that the series $\sum_ns_n$ converges and $\sum_n{\sigma}_n^2<\infty$. In particular $s_n\to0$ as $n\to\infty$.
Since $x$ and $y$ are non-conjugate along ${\gamma}$, for all non-zero $k\in K$, we have
$$
\|k\|^2-s(k)=q(k)=Q(v(k))+\|k'\|^2>0
$$
where $k'=k-k({\lambda},v(k))$.
Hence $1+s_n>0$ for all $n$.
Denote by $\nu$ the standard Gaussian distribution on ${\mathbb R}$.
Then, for all $s\in(-1,\infty)$ and all ${\sigma}\in{\mathbb R}$, we have
\begin{equation}\label{TSB}
\int_{\mathbb R} e^{{\sigma} x}e^{-sx^2/2}\nu(dx)=\frac{1}{\sqrt{1+s}}e^{{\sigma}^2/(2(1+s))}.
\end{equation}
Hence, for $p>1$ sufficiently close to $1$, for all $z\in{\mathbb R}^d$, we have
$$
{\mathbb E}(e^{pS(z)/2})=e^{ps({\tau}(z))}\prod_n\int_{\mathbb R} e^{p{\sigma}_n x}e^{-ps_nx^2/2}\nu(dx)=e^{s({\tau}(z))}\prod_n\frac{1}{\sqrt{1+ps_n}}e^{p^2{\sigma}_n^2/(2(1+ps_n))}<\infty.
$$ 
Let $\phi$ be a continuous linear functional on ${T_\g\oxy}$ and write $\tilde\phi=v(k^*)$, where $k^*=k({\lambda},\tilde\phi)$.
Then
$$
\phi(v(e_n))=Q(\tilde\phi,v(e_n))=q(k^*,e_n)=\<k^*,e_n\>q(e_n).
$$
From \eqref{TSB}, we see that under $\tilde\mathbb{P}$, the random variables $\<W,e_n\>$ are independent zero-mean Gaussians of variance $(1+s_n)^{-1}=q(e_n)^{-1}$.
Now
$$
\phi(Y)=\phi(v(W))=\sum_n\<W,e_n\>\phi(v(e_n))
$$
so the law $\mu_\g$ of $Y$ on ${T_\g\oxy}$ under $\tilde\mathbb{P}$ is a zero-mean Gaussian measure with 
$$
\int_{T_\g\oxy}\phi(v)^2\mu_\g(dv)={\mathbb E}(\phi(Y)^2)=\sum_n\phi(v(e_n))^2/q(e_n)=\sum_n\<k^*,e_n\>^2q(e_n)=q(k^*)=Q(\tilde\phi)
$$
as claimed.
Finally, by Proposition \ref{2.8}, for ${\alpha}\inT^*_{\g_s}M$ and ${\beta}\inT^*_{\g_t}M$, we have
$$
\int_{T_\g\oxy}\<{\alpha},v_s\>\<{\beta},v_t\>\mu_\g(dv)=Q(v^{{\alpha},s},v^{{\beta},t})=\<{\alpha},C(s,t){\beta}\>
$$
so $\mu_\g$ also has the claimed covariance. The uniqueness statements are standard.
\end{proof}


