Let $(M,g)$ be a Riemannian manifold, and fix $x_0\in M$. For $x\in M$, let $\tau(x) = d_g(x_0,x)$. For any unit-length geodesic $\gamma$ with $\gamma(0) = x_0$, we have the identity $$ \gamma'(t_0) = \text{grad }\tau|_{\gamma(t_0)} $$ for sufficiently small $t_0> 0$. This can be seen by noting that $\|\text{grad }\tau\|_g\le 1$ essentially by the triangle inequality, while $$\langle\text{grad }\tau,\gamma'(t_0)\rangle_g = \gamma'(t_0)(\tau) = 1$$ since $\tau(\gamma(t)) = t$ for small $t>0$, so the desired equality follows from the fact that equality in Cauchy-Schwarz is attained.
The geodesic flow on $TM$ can also be viewed as a Hamilton flow on $T^*M$ with respect to the Hamiltonian $H(x,\xi) = \frac{1}{2}|\xi|^2_g$, where $|\xi|^2_g$ is the dual metric on $T^*M$ induced by $g$. Then $\gamma(t_0) = x(t_0)$, where $(x(t),\xi(t))$ is the trajectory of some Hamiltonian flow with $x(0) = x_0$ and $\xi(0)$ satisfying $|\xi(0)|^2_g = 1$, and $\gamma'(t_0)$ can be identified with $\xi(t_0)$ via the metric identification. Furthermore $\tau(x)$ denotes the smallest $t\ge 0$ for which $x(0) = x_0$ and $x(t) = x$ for some choice of $\xi(0)$ satisfying $|\xi(0)|^2_g = 1$.
Suppose we have chosen coordinates $(x_1,\dots,x_n)$ on $M$ near $x_0$ (which gives induced coordinates on $T^*M$), with the metric and dual metric written as $g_{ij}(x)$ and $g^{ij}(x)$. If we write $\gamma'(t_0) = \sum_{i=1}^n{\gamma_i'(t_0)\partial_i}$, then we can write $\xi_i(t_0) = \sum_{j=1}^n{g_{ij}(x(t_0))\gamma'_j(t_0)}$, and so for $v = \sum_{i=1}^n{v_i\partial_i}\in T_{x(t_0)}M$ we have $$\langle\gamma'(t_0),v\rangle_g = \sum_{i,j=1}^n{g_{ij}(x(t_0))\gamma'_i(t_0)v_j} = \sum_{j=1}^n{\xi_j(t_0)v_j} = \xi(t_0)\cdot v.$$ On the other hand by definition we have $$ \langle\text{grad }\tau|_{x(t_0)},v\rangle_g = d\tau|_{x(t_0)}\cdot v$$ so we have the identity $$\xi(t_0) = d\tau|_{x(t_0)}.$$
Question: Does a similar property hold for more general Hamiltonian flows, say for Hamiltonians $H(x,\xi)$ which are homogeneous of degree $2$ in $\xi$? Precisely, if we have a general $H$ homogeneous of degree $2$, and we let $$\tau(x) = \inf\{t\ge 0\,:\,x = x(t)\text{ for some trajectory }(x(t),\xi(t))\text{ with }x(0) = x_0, H(x_0,\xi(0)) = \frac{1}{2}\}$$ and we assume, for example, that $\xi\mapsto \partial_{\xi}H(x,\xi)$ is invertible on $T^*_xM\backslash\{0\}$ for all $x$ (to guarantee uniqueness of trajectories, at least locally), then for a trajectory $(x(t),\xi(t))$ with $x(0) = x_0$ and $H(x_0,\xi(0)) = 1/2$, is it true that $\xi(t_0) = d\tau|_{x(t_0)}$ for all sufficiently small $t_0> 0$?
You expectation is true for homogeneous Hamiltonian functions $H$, as in this case your function $\tau$ is closely related to an action functional $S$ for the corresponding Hamiltonian flow. Now, it is one's of Hamilton's greats achievements to have pursued the optical-mechanical analogy beyond "Fermat's least time principle is analogous to Maupertuis' least action principle": Huygens' principle should also have a mechanical counterpart, which is the Hamilton-Jacobi principle. Notice that your Riemannian context could serve as a model for a light ray moving (at constant unit speed) inside a homogeneous medium. Different from -- in fact dual to -- the light ray perspective is Huygens' wave front perspective; a wave front at time $t$ emanating from a point $x_0 \in M$ is the set of points $\tau(x) = t$. Huygens' principle stated essentially that the velocity vector of the light ray at $x \in \tau^{-1}(t)$ should be dual to the covector $\mathrm{d}\tau_x$. The Hamilton-Jacobi principle generalizes this to mechanical systems by considering the action functional $S$ instead of the time functional $\tau$ and establishes that $\mathrm{d}S_x = \xi$ (using your notations and under some specifications). All of this is explained in Arnold's Mathematical Methods of Classical Mechanics, but I shall try to summarize/sketch this below for the present context.
Given a path $\gamma : I = [0, T] \to T^*M$, one define its action as: $$S[\gamma] = \int_I [\lambda(\dot{\gamma}(t)) - H(\gamma(t))]dt = \int_I \left[\sum_j \xi_j \dot{x}^j - H(x, \xi)\right]dt $$ where $\lambda$ is the tautological 1-form, namely in Darboux coordinates $(x^i, \xi_j)$, $\lambda = \sum_j \xi_j dx^j$. If $\gamma$ is a solution to Hamilton's equations, then $\dot{x} = \partial_{\xi}H$; if $H$ is homogeneous in $\xi$ of degree $m$, then Euler's homogeneity theorem implies that $\lambda(\dot{\gamma}(t)) = m H(\gamma(t))$. Since $H$ is constant under the Hamiltonian flow, we obtain $$ (\star) \qquad S[\gamma] = (m-1)H(\gamma(0)) T .$$ Observe that if $(x(t), \xi(t))$ is a solution of Hamilton's equations such that $x(T) = x$, then $(x'(t), \xi'(t)) := (x(tT), T^{1/(m-1)}\xi(tT))$ is a solution which satisfies $x'(1) = x$ and $$(\star\star) \qquad H(x'(0), \xi'(0)) = H(x(0), \xi(0)) T^{m/(m-1)} .$$
Assume that the Legendre transformation $\mathcal{L}_x : T_x^*M \to T_xM : \xi \mapsto v := \partial_{\xi}H(x, \xi)$ is invertible for $x$ near $x_0$. We define the Lagrangian function $$ L : TM \to \mathbb{R} : L(x,v) = \lambda_{\mathcal{L}_x^{-1}(v)}(v) - H(\mathcal{L}_x^{-1}(v)) $$ and we consider a somewhat different action function $S'$ on path $\gamma'$ in (a neighborhood of $x_0$ in) $M$, $$ S'[\gamma'] = \int_I L(\gamma', \dot{\gamma'}) dt . $$ It turns out that $\gamma'$ satisfies Euler-Lagrange equations if and only if the corresponding Legendre transform $\gamma$ satisfies Hamilton's equations, hence $S[\gamma] = S'[\gamma']$. In a small neighborhood $U$ of $x_0$ in $M$, set $$ S : U \to \mathbb{R} : S(x) := \mathrm{inf} \, \left\{ S'[\gamma'] \, | \, \gamma' : [0,1] \to U, \; \gamma'(0) = x_0, \; \gamma'(1) = x \right\} .$$ The infima are attained when $U$ is sufficiently small, the minimising curves satisfying the Euler-Lagrange equations.
Now consider variations of $S'[\gamma']$ among curves $\gamma' : [0,1] \to U$ which begin at $x_0$ and which all satisfy Euler-Lagrange equations (so that $S'[\gamma'] = S(\gamma'(1))$). We easily get that (see for example this derivation) $$ \mathrm{d}S_{\gamma'(1)}\left[\frac{d(\delta \gamma')}{dt}(1)\right] = \delta S'_{\gamma'}[\delta \gamma'] = \frac{\partial L}{\partial v}(\gamma'(1), \dot{\gamma'}(1)) \; (\delta \gamma')(1)$$ The correspondence $T_xM \to T_x^*M : v \mapsto (\partial_v L)(x,v)$ happens to be the inverse of the previous Legendre transform. We deduce $\mathrm{d}S_{x'(1)} = \xi'(1)$, where $(x'(t), \xi'(t)) : [0,1] \to T^*U$ is the corresponding solution of Hamilton's equations.
The path $(x(t), \xi(t)) = (x'(t/T), T^{-1/(m-1)}\xi'(t/T))$ is the solution of Hamilton's equations you are really interested in. Using $(\star)$ on $(x', \xi')$ and $(\star\star)$ to switch to $(x, \xi)$, setting $x = x(T)$ so that $T = \tau(x)$, we compute $$ T^{1/(m-1)}\xi(T) = \xi'(1) = \mathrm{d}\left( (m-1)H(\xi(0))\tau(x)^{m/(m-1)} \right) = m T^{1/(m-1)}H(\xi(0))\mathrm{d}\tau_{x}$$ Hence $\xi(T) = m H(x(0), \xi(0)) \mathrm{d}\tau_{x(T)}$. When $m=2$ and $H(x(0), \xi(0)) = 1/2$, we recover the relation you were expecting.