Let me establish that the computation of \rho is not in question here, and q_{N} - q_{1} and is not relevant to discussion at all.
The current sampler is not trying to estimate or otherwise compute q_{N} - q_{1} for exact trajectories or numerical trajectories. The difference in positions appeared in the original No-U-Turn sampler as a heuristic based on the special geometry of the real numbers. The difference between the end points of exact trajectories doesn’t generalize beyond the real numbers. Moreover the difference between the end points of a numerical trajectory isn’t well-defined because of its dependence on a particular discretization.
The integrated momenta,
\rho = \int_{0}^{T} \mathrm{d} t \, p(z(t)),
however, is a proper geometric object. Take care as this is not a regular integral, both because of the implicit z(t) dependence but also because p and \rho are not real numbers but rather more general geometric objects. On Euclidean manifolds, where the No-U-Turn heuristic is well-defined, and exact trajectories these two definitions give the same termination criteria. The integrated momenta, however, generalizes to any manifold and hence is better suited for Hamiltonian Monte Carlo. In other words the original No-U-Turn criterion always approximated the new one, not the other way around.
For any function over phase space the natural estimator of its temporal expectation is the empirical average of evaluations over each discrete state,
\int_{0}^{T} \mathrm{d} t \, f(z(t)) = \frac{1}{N} \sum_{n = 1}^{N} f(z_{n}).
These estimators are unbiased as N \rightarrow \infty and enjoy nice convergence properties for finite N due to the nature of symplectic integrators.
To approximate the integrated momenta we just utilize the usual estimators for temporal expectations giving the current algorithm. For the very technically minded when computing \rho in practice we compute its components in a coordinate system, which are given by D functions defined as the compositions of the canonical momenta operator p: T^{*}Q \rightarrow T_{q}^{*} Q with local cotangent space coordinate functions, p_{i} : T_{q}^{*} Q \rightarrow \mathbb{R}^{i}.
The difference between the estimator for \rho and the value of q_{N} - q_{0} for a numerical trajectory is irrelevant because we are not trying to estimate the distance between the end points anymore. The difference between q_{N} - q_{0} for exact and numerical trajectories also doesn’t matter, nor does the difference between \rho or estimates of \rho and those two Euclidean distances. The distances are an artifact of a the original heuristic from which we have long since moved passed.
The current dynamic termination criterion is based on the geometric calculation p^{\sharp} \cdot \rho, and we approximate that second term using a numerical trajectory generated by a symplectic integrator per standard practice.