Proof sketch. Writing down the iteration 3 in two consecutive steps and eliminating the momentum difference yields, for \(\Delta q_k \mathrel{\vcenter{:}}= x_k - y_k\), the three-term recurrence \[\label{eqn-three-term-recurr-in-x-gaussian} \Delta q_{k+1} = (2I-h^2 \Sigma^{-1})\Delta q_k - \Delta q_{k-1}\,,\quad k \geq 1.\tag{20}\] Deriving the corresponding recurrence for the propagator matrices reveals the Chebyshev polynomials. ◻

Proof of Example 1. We assume \(m\) is chosen so that \(mh\) avoids the discrete set of resonant values at which some eigenmode fails to contract. Away from these, the leap-frog iteration contracts in every mode and the chain converges to \(\pi_h\). As these values are exceptional, we omit the condition from the example. According to [50], the invariant distribution for a target distribution \(\pi=\mathcal{N}(\mu, \Sigma)\) is given by \(\pi_h = \mathcal{N}(\mu, \Sigma_{\infty})\), where \(\Sigma_\infty^{-1} = \Sigma^{-1}(I- \frac{h^2}{4} \Sigma^{-1})\). Consider the coupling between \(\pi_h\) and \(\pi\) as \(X = \Sigma^{1/2}Z + \mu\) and \(Y = \Sigma_\infty^{1/2}Z+\mu\) with \(Z \sim \mathcal{N}(0,I)\). Then, \[W_{2,\ell^\infty}^2(\pi_h, \pi) \leq \mathbb{E}[|X-Y|_{\ell^\infty}^2] = \mathbb{E}[|(\Sigma^{1/2} - \Sigma_\infty^{1/2})Z|_{\ell^\infty}^2].\] Entries of \((\Sigma^{1/2}-\Sigma_\infty^{1/2})Z\) are \(|\Sigma^{1/2}-\Sigma_\infty^{1/2}|_{\ell^2}^2\)-subgaussian. And \[|\Sigma^{1/2}-\Sigma_\infty^{1/2}|_{\ell^2} = \max_{1\leq i \leq d}\left|\frac{1}{\omega_i} - \frac{1}{\hat{\omega}_i}\right| = O(\sqrt{\beta}h^2)\,.\] By Lemma 4, we have \(W_{2,\ell^\infty}(\pi_h, \pi) = O(\sqrt{\beta}h^2 \sqrt{\log(2d)})=O(h\sqrt{\log(2d)})\) as \(h < 1/\sqrt{\beta}\). ◻

Proof of Proposition 8. Part (a). We derive the correspondence with multivariate damped matrix polynomials from three-term recurrence relations, similar to Proposition 5 but with varying Hessians. With the same noise applied to both trajectories, the BAOAB iteration 4 for the differences \(\Delta q_k \mathrel{\vcenter{:}}= x_k - y_k\) and \(\Delta p_k \mathrel{\vcenter{:}}= p_k - p'_k\) is \[\label{eqn-inner-loop-difference} \begin{align} \Delta q_{k+1} &= \Delta q_k + \tfrac{h}{2}(1+\eta)\Delta p_k - \tfrac{h^2}{4}(1+\eta)H_k \Delta q_k\,,\\ \Delta p_{k+1} &= \eta\big(\Delta p_k - \tfrac{h}{2}H_k\Delta q_k\big) - \tfrac{h}{2}H_{k+1} \Delta q_{k+1}\,. \end{align}\tag{30}\] The leap-frog iteration 3 is the special case \(\eta = 1\), so it suffices to treat 30 for general \(\eta\). We write down 30 also for \((\Delta q_k, \Delta p_k)\) and eliminate the momentum \(\Delta p_k\) in two consecutive steps by matching its coefficients. We then obtain \[\label{eqn-three-term-recurrence-baoab-general} \Delta q_{k+1} = (1+\eta)(I - \tfrac{h^2}{2}H_k)\Delta q_k - \eta \Delta q_{k-1}\,, \qquad k \geq 1\,.\tag{31}\] At \(\eta = 1\), the recurrence recovers the relation for the leap-frog scheme, \(\Delta q_{k+1} = (2I - h^2 H_k)\Delta q_k - \Delta q_{k-1}\).

We prove 23 by induction. It holds for \(k=0,1\) by the initial conditions of \(\tilde{T}_k^\eta\) and \(\tilde{U}_k^\eta\) in Definition 1. For \(k \geq 1\), the induction hypothesis and 31 give \[\begin{align} \Delta q_{k+1} &= \big((1+\eta)(I -\tfrac{h^2}{2}H_k)\tilde{T}_k^\eta - \eta \tilde{T}_{k-1}^\eta\big)\Delta q_0\\ & \quad + \big((1+\eta)(I - \tfrac{h^2}{2}H_k)h\tilde{U}_{k-1}^\eta - \eta h\tilde{U}_{k-2}^\eta\big)\Delta p_0\,, \end{align}\] and the three-term recurrences for \(\tilde{T}_k^\eta\) and \(\tilde{U}_k^\eta\) in 21 and 22 identify the right-hand side as 23 , completing the induction.

Part (b). We iterate the single-scheme identity at \(\eta = 1\) across the HMC outer loops. For \(N \geq 1\), since both trajectories share the refreshment \(\xi_{N-1}\), the identity gives \[\begin{align} \Delta q_{Nm+k} &= \tilde{T}_m^1\big(\{I -\tfrac{h^2}{2}H_i\}_{i=(N-1)m+k}^{Nm+k-1}\big)\,\Delta q_{(N-1)m + k}\\ &= \prod_{j \in\{k + im\}_{i=0}^{N-1}}^{\longleftarrow} \tilde{T}_m^1\big(\{I - \tfrac{h^2}{2}H_l\}_{l=j}^{j+m-1}\big)\,\Delta q_k\,, \end{align}\] and the conclusion follows by substituting the expression for \(\Delta q_k\) from 23 . ◻

Proof of Lemma 1. The definitions of \(\tilde{T}_k^\eta\) and \(\tilde{U}_k^\eta\) in 21 and 22 give the matrix recurrences \[\begin{align} \tag{34} \tilde{T}_{k+1}^\eta &= (1+\eta)(I - \tfrac{h^2}{2}H_{k})\tilde{T}_k^\eta - \eta \tilde{T}_{k-1}^\eta\,,\\ \tag{35} \tilde{U}_{k+1}^\eta &= (1+\eta)(I - \tfrac{h^2}{2}H_{k+1})\tilde{U}_k^\eta - \eta \tilde{U}_{k-1}^\eta\,. \end{align}\] We establish the expansions 32 and 33 by induction, the cases \(k=0,1\) following from the initial conditions.

Constant terms. The constant coefficients satisfy \(t_0^{(k+1)} = (1+\eta) t_0^{(k)} - \eta t_0^{(k-1)}\) and \(u_0^{(k+1)} = (1+\eta) u_0^{(k)} - \eta u_0^{(k-1)}\), with \(t_0^{(0)} = t_0^{(1)} = 1\), \(u_{0}^{(-1)} = 0\), and \(u_0^{(0)} = (1+\eta)/2\). Solving these gives \(t_0^{(k)} = 1\) and \(u_0^{(k)} = \tfrac12 (1+\eta)\frac{\eta^{k+1}-1}{\eta-1}\), the latter read as \(k+1\) when \(\eta = 1\).

Higher-order terms. Assume 32 and 33 hold up to index \(k\). Passing to \(k+1\) via 34 and 35 , the only new matrix products arise from \(-(1+\eta)\tfrac{h^2}{2}H_k\tilde{T}_k^\eta\) and \(-(1+\eta)\tfrac{h^2}{2}H_{k+1}\tilde{U}_k^\eta\), which prepend \(H_k\) (resp. \(H_{k+1}\)) to the existing terms. Applied to the order-\(j\) terms of \(\tilde{T}_k^\eta\) and \(\tilde{U}_k^\eta\), this produces the order-\((j+1)\) contributions \[(-1)^{j+1}h^{2(j+1)}\tfrac{1+\eta}{2}\!\!\sum_{k-1\geq i_1>\dots>i_j \geq 0}\!\! \tilde{t}_{i_1,\dots,i_j}^{(k)} H_k H_{i_1}\cdots H_{i_j}\,,\] \[(-1)^{j+1}h^{2(j+1)}\tfrac{1+\eta}{2}\!\!\sum_{k\geq i_1>\dots>i_j\geq1}\!\!\tilde{u}_{i_1,\dots,i_j}^{(k)}H_{k+1}H_{i_1}\cdots H_{i_j}\,,\] together with the new order-one terms \(-(1+\eta)\tfrac{h^2}{2}t_0^{(k)}H_k\) and \(-(1+\eta)\tfrac{h^2}{2}u_0^{(k)}H_{k+1}\). In each new product the prepended factor carries a strictly larger index than all of \(i_1 > \dots > i_j\), so the matrices remain in descending index order, and 32 and 33 hold at \(k+1\).

Nonnegativity. We match matrix products across the recurrences and induct on \(j\), the number of factors. The case \(j=0\) was settled above. A product \(H_{i_1}\cdots H_{i_j}\) first appears in \(\tilde{T}_{i_1+1}^\eta\) and \(\tilde{U}_{i_1}^\eta\), with coefficients, for \(i_1\geq 1\), \[\label{eqn-coeff-first-emerge} \tilde{t}_{i_1,\dots,i_j}^{(i_1+1)} = \tfrac{1+\eta}{2}\,\tilde{t}_{i_2,\dots,i_j}^{(i_1)}\,, \qquad \tilde{u}_{i_1,\dots,i_j}^{(i_1)} = \tfrac{1+\eta}{2}\,\tilde{u}_{i_2,\dots,i_j}^{(i_1-1)}\,,\tag{36}\] where for \(j=1\) the right-hand factors are read as the constant coefficients \(t_0^{(i_1)}\) and \(u_0^{(i_1-1)}\). Besides the induction, we also note that the first-kind base case is \(\tilde{t}_0^{(1)}=\frac{1+\eta}{4} \geq 0\), the coefficient of \(-h^2H_0\) in \(\tilde{T}_1^\eta\). By the induction hypothesis on \(j-1\), these initial coefficients are nonnegative. Tracking the coefficient of this same product \(H_{i_1}\cdots H_{i_j}\) in the higher-index polynomials, it satisfies \[\tilde{t}_{i_1,\dots,i_j}^{(i_1+2)} = (1+\eta)\tilde{t}_{i_1,\dots,i_j}^{(i_1+1)}\,, \qquad \tilde{u}_{i_1,\dots,i_j}^{(i_1+1)} = (1+\eta)\tilde{u}_{i_1,\dots,i_j}^{(i_1)}\,,\] and, for \(k \geq i_1+2\), \[\tilde{t}_{i_1,\dots,i_j}^{(k+1)} = (1+\eta)\tilde{t}_{i_1,\dots,i_j}^{(k)} - \eta\tilde{t}_{i_1,\dots,i_j}^{(k-1)}\,, \qquad \tilde{u}_{i_1,\dots,i_j}^{(k)} = (1+\eta)\tilde{u}_{i_1,\dots,i_j}^{(k-1)} - \eta\tilde{u}_{i_1,\dots,i_j}^{(k-2)}\,.\] Solving these recurrences gives, for \(k \geq i_1+1\), \[\tilde{t}_{i_1,\dots,i_j}^{(k)} = \frac{\eta^{k-i_1}-1}{\eta-1}\,\tilde{t}_{i_1,\dots,i_j}^{(i_1+1)}\,, \qquad \tilde{u}_{i_1,\dots,i_j}^{(k-1)} = \frac{\eta^{k-i_1}-1}{\eta-1}\,\tilde{u}_{i_1,\dots,i_j}^{(i_1)}\,,\] which are nonnegative, with the \(\eta=1\) values read as limits as before. ◻

Proof of Proposition 17. Taking \(H_i \mathrel{\vcenter{:}}= I\) in Lemma 1 collapses the expansions onto their aggregate coefficients, \[\begin{align} \tilde{T}_k^\eta(\{(1- \tfrac{h^2}{2})I\}_{i=0}^{k-1}) &= t_0^{(k)} I + \sum_{j=1}^k (-1)^j h^{2j}t_j^{(k)}I\,, \\ \tilde{U}_k^\eta(\{(1 - \tfrac{h^2}{2})I\}_{i=1}^{k}) &= u_0^{(k)}I + \sum_{j=1}^k (-1)^jh^{2j}u_j^{(k)}I\,. \end{align}\] At \(\eta = 1\) with equal arguments, the multivariate polynomials reduce to the classical Chebyshev polynomials \(T_k\) and \(U_k\) (as noted after Definition 1), so \[\label{eqn-reduced-matrix-polynomial-coefficients} T_k(1 - \tfrac{h^2}{2})\,I = \sum_{j=0}^k(-1)^j h^{2j}t_j^{(k)}I\,, \qquad U_k(1-\tfrac{h^2}{2})\,I = \sum_{j=0}^k(-1)^jh^{2j}u_j^{(k)}I\,.\tag{37}\] The hypergeometric forms of \(T_k\) and \(U_k\) [51] give explicit expressions for \(T_k(1-x)\) and \(U_k(1-x)\) at \(x = h^2/2\), and matching coefficients yields the stated formulas. ◻

Proof of Proposition 18. We derive iterative inequalities for the coefficients from their recurrences and conclude by induction. We treat \(t_j^{(k)}\), and the argument for \(u_j^{(k)}\) is analogous. For \(j \geq 1\), Lemma 1 gives \[t_j^{(k)} = \sum_{k-1 \geq i_1 > \dots > i_j \geq 0}\tilde{t}_{i_1,\dots,i_j}^{(k)} = \sum_{k-1\geq i_1 > \dots > i_j \geq 0}\frac{\eta^{k-i_1}-1}{\eta - 1}\,\tilde{t}_{i_1,\dots,i_j}^{(i_1+1)}\,.\] By 36 and the initial conditions \(\tilde{t}_0^{(1)} = \frac{1+\eta}{4}\) (which is the coefficient of \(-h^2H_0\)) and \(t_0^{(0)}= 1\), we have \(\tilde{t}_{i_1,\dots,i_j}^{(i_1+1)} \leq \tfrac12(1+\eta)\,\tilde{t}_{i_2,\dots,i_j}^{(i_1)}\), so \[\begin{align} t_j^{(k)} &\leq \tfrac12(1+\eta)\sum_{k-1\geq i_1\geq j-1}\frac{\eta^{k-i_1}-1}{\eta-1}\left(\sum_{i_1-1\geq i_2 > \dots > i_j \geq 0} \tilde{t}_{i_2,\dots,i_j}^{(i_1)}\right)\\ & = \tfrac12(1+\eta)\sum_{k-1\geq i_1 \geq j-1}\frac{\eta^{k-i_1}-1}{\eta-1}\,t_{j-1}^{(i_1)} \leq \frac{1}{1-\eta}\sum_{k-1\geq i_1 \geq j-1} t_{j-1}^{(i_1)}\,, \end{align}\] using the definition of \(t_{j-1}^{(i_1)}\) in the equality. We now induct on \(j\). The case \(j = 0\) holds since \(t_0^{(k)} = 1\) for all \(k\). Assuming ?? for some \(j\) and all \(k \geq j\), then for \(k \geq j+1\), \[t_{j+1}^{(k)}\leq \frac{1}{1-\eta}\sum_{k-1\geq i_1 \geq j} \Big(\frac{1}{1-\eta}\Big)^j\binom{i_1}{j} = \Big(\frac{1}{1-\eta}\Big)^{j+1}\binom{k}{j+1}\,,\] where the last step uses Pascal’s rule \(\sum_{k-1\geq i_1 \geq j}\binom{i_1}{j} = \binom{k}{j+1}\). This proves the bound for \(t_j^{(k)}\).

For \(u_j^{(k)}\), the same steps give the iterative inequality \(u_j^{(k)} \leq \frac{1}{1-\eta}\sum_{k \geq i_1\geq j} u_{j-1}^{(i_1-1)}\). Combined with the initial condition \(u_0^{(k)} = \tfrac12(1+\eta)\frac{\eta^{k+1}-1}{\eta-1}\) from Lemma 1 and the identity \(\sum_{k\geq i_1\geq j}\binom{i_1-1}{j-1} = \binom{k}{j}\), this yields the stated bound. ◻

Proof of Proposition 19. Alternating signs. In each factor \(\tilde{T}_k^1\) and \(\tilde{U}_k^1\) with arguments \(I - \tfrac{h^2}{2}H_i\), a product of \(j\) matrices carries the sign \((-1)^j\). This sign is multiplicative across the product of factors in ?? and ?? , so every order-\(j\) product carries \((-1)^j\), and the coefficients are nonnegative.

Reduction to Chebyshev polynomials. Taking \(H_i \mathrel{\vcenter{:}}= I\) collapses the expansions onto their aggregate coefficients, and at \(\eta = 1\) the multivariate polynomials reduce to the classical \(T_k\) and \(U_k\), giving \[\begin{align} \big(T_m(1-\tfrac{h^2}{2})\big)^N T_k(1-\tfrac{h^2}{2})\,I &= a_0^{(Nm+k)}I + \sum_{j=1}^{Nm+k}(-1)^j h^{2j}a_j^{(Nm+k)}I\,,\\ \big(T_m(1-\tfrac{h^2}{2})\big)^N U_{k-1}(1-\tfrac{h^2}{2})\,I &= b_0^{(Nm+k-1)}I + \sum_{j=1}^{Nm+k-1}(-1)^j h^{2j}b_j^{(Nm+k-1)}I\,. \end{align}\] Writing \(\theta = \arccos(1-h^2/2)\), so that \(T_k(1-\tfrac{h^2}{2}) = \cos(k\theta)\) and \(U_{k-1}(1-\tfrac{h^2}{2}) = \sin(k\theta)/\sin\theta\), and expanding \[\label{eqn-cos-power-N} \big(\cos(m\theta)\big)^N = \frac{1}{2^N}\sum_{i=0}^N\binom{N}{i}\exp(\mathrm{i}(2i-N)m\theta)\,,\tag{38}\] we pair conjugate exponentials into cosines and apply \(2\cos x\cos y = \cos(x-y)+\cos(x+y)\) and \(2\cos x\sin y = \sin(x+y) - \sin(x-y)\). Rewriting the resulting trigonometric sums back in terms of Chebyshev polynomials expresses \(\big(T_m\big)^N T_k\) as a binomially weighted sum of \(T_{\,\cdot}\) at shifted indices \(2im \pm k\) (for \(N\) even) or \((2i-1)m \pm k\) (for \(N\) odd), all with positive weights. For \(\big(T_m\big)^N U_{k-1}\) the \(\sin(x+y) - \sin(x-y)\) identity instead produces a signed sum of \(U_{\,\cdot}\) at indices \(2im \pm k - 1\) (for \(N\) even) or \((2i-1)m \pm k - 1\) (for \(N\) odd), the \(x-y\) terms entering with a minus sign. Matching coefficients with \(t_j^{(\cdot)}\) and \(u_j^{(\cdot)}\) then writes each \(a_j^{(Nm+k)}\) and \(b_j^{(Nm+k-1)}\) as a \(2^{-N}\)-weighted sum of inner-loop aggregate coefficients \(t_j^{(i)}\) and \(u_j^{(i)}\), with binomial weights and indicator factors \(\mathbf{1}_{j \leq i}\). The sum for \(a_j^{(Nm+k)}\) has all positive weights, while that for \(b_j^{(Nm+k-1)}\) inherits the minus signs above. Explicit formulas can be read off in this way, but the main results need only the bound below.

Bound. Each \(a_j^{(Nm+k)}\) is a sum of the form \(\tfrac{1}{2^N}\sum_i w_i\, t_j^{(i)}\mathbf{1}_{j \leq i}\) with nonnegative binomial weights \(w_i\) satisfying \(\sum_i w_i = 2^N\). Since the \(t_j^{(i)}\) are nonnegative by Lemma 1, bounding each term by \(t_j^{(i)}\mathbf{1}_{j \leq i} \leq \max_{j \leq i \leq Nm+k} t_j^{(i)}\) gives \[a_j^{(Nm+k)} \leq \frac{1}{2^N}\Big(\sum_i w_i\Big)\max_{j \leq i \leq Nm+k} t_j^{(i)} = \max_{j \leq i \leq Nm+k} t_j^{(i)}\,.\] The same argument bounds \(b_j^{(Nm+k-1)}\). Here some \(u_j^{(i)}\) enter with negative signs, but since the \(u_j^{(i)}\) are nonnegative, discarding the subtracted terms only increases the bound, and the remaining terms are controlled as above by \(\max_{j \leq i \leq Nm+k-1} u_j^{(i)}\). ◻

Proof of Proposition 20. We interpret the polynomials with general \(\{H_i\}\) as the propagator of a leap-frog iteration with varying quadratic potentials, which lets us bound their \(\ell^2\) norms through the contraction and boundedness of the discrete-time dynamics, analyzed via their continuous-time interpolation.

Connection with dynamics. Consider the time-inhomogeneous linear dynamics driven by \(\{H_i\}\), \[\label{eqn-linear-inner-loop} \begin{align} q_{k+1} &= q_k + hp_k - \tfrac{h^2}{2}H_k q_k\,,\\ p_{k+1} &= p_k - \tfrac{h}{2}H_kq_k - \tfrac{h}{2}H_{k+1}q_{k+1}\,, \end{align}\tag{39}\] which is the leap-frog scheme 3 with the quadratic potential \(V_k(q) = \tfrac12 q^\top H_k q\) at step \(k\). For two initial states \((x_0, p_0)\), \((y_0, p_0')\), the differences \(\Delta q_k \mathrel{\vcenter{:}}= x_k - y_k\), \(\Delta p_k \mathrel{\vcenter{:}}= p_k - p'_k\) satisfy \[\label{eqn-linear-inner-loop-difference} \begin{align} \Delta q_{k+1} &= (I- \tfrac{h^2}{2} H_k)\Delta q_k + h \Delta p_k\,,\\ \Delta p_{k+1} &= \Delta p_k - \tfrac{h}{2}H_k \Delta q_k - \tfrac{h}{2}H_{k+1}\Delta q_{k+1}\,. \end{align}\tag{40}\] This matches 30 in the proof of Proposition 8, now with general \(\{H_i\}\) in place of the integrated Hessians. As that proof does not use the specific Hessian structure, the same argument gives the three-term recurrence \(\Delta q_{k+1} = (2I - h^2 H_k)\Delta q_k - \Delta q_{k-1}\) for \(k \geq 1\) and, by induction, \[\label{eqn-prop-proof-q-dependence-q0-p0} \Delta q_{k} = \tilde{T}_{k}^1\,\Delta q_0 + h\tilde{U}_{k-1}^1\,\Delta p_0\,.\tag{41}\] Thus \(\tilde{T}_k^1\) and \(h\tilde{U}_{k-1}^1\) are the propagators of 39 . We bound their \(\ell^2\) norms by choosing initial conditions that isolate each, setting \(\Delta p_0 = 0\) for \(\tilde{T}_k^1\), and \(\Delta q_0 = 0\) for \(h\tilde{U}_{k-1}^1\).

\(\ell^2\) norm of \(\tilde{T}_k^1\). Adapting the proof of [42] to a fixed matrix sequence \(\{H_i\}\) rather than the leap-frog scheme of the target potential, we show that the position difference still contracts when the potential is decoupled from the trajectory. We include it both for completeness and to set up the \(\tilde{U}_k^1\) bound. Introduce the continuous-time interpolation of 39 , \[\label{eqn-linear-inner-loop-interpolation} \begin{align} \frac{{\rm d} q_t}{{\rm d }t} &= p_{\lfloor t \rfloor_h} - (t - \lfloor t \rfloor_h)H_{\lfloor t \rfloor_h}q_{\lfloor t \rfloor_h}\,,\\ \frac{{\rm d}p_t}{{\rm d}t} &= -\tfrac12(H_{\lfloor t \rfloor_h} q_{\lfloor t \rfloor_h} + H_{\lceil t \rceil_h} q_{\lceil t \rceil_h})\,, \end{align}\tag{42}\] with \(\lfloor t \rfloor_h \mathrel{\vcenter{:}}= \lfloor t/h\rfloor h\) and \(\lceil t \rceil_h \mathrel{\vcenter{:}}= \lceil t/h \rceil h\). Take two trajectories from \((x_0, p_0)\) and \((y_0, p_0')\) with \(p_0 = p_0'\), set \(z_t \mathrel{\vcenter{:}}= x_t - y_t\), \(w_t \mathrel{\vcenter{:}}= p_t - p_t'\), and at each grid point \(t_k = kh\) define \(\Phi_{kh} \mathrel{\vcenter{:}}= H_{k} z_{kh}\) and \(\varphi_{kh} \mathrel{\vcenter{:}}= \Phi_{kh}\cdot z_{kh}\), with \(\cdot\) the inner product. Then \(w_0 = 0\) and \(|\Phi_t|_{\ell^2}^2 \leq \beta \varphi_t\), and 42 gives \[\label{eqn-linear-inner-loop-interpolation-difference} \begin{align} \frac{{\rm d}z_t}{{\rm d}t} &= w_{\lfloor t \rfloor_h} - (t - \lfloor t \rfloor_h) \Phi_{\lfloor t \rfloor_h}\,,\\ \frac{{\rm d}w_t}{{\rm d}t} &= -\tfrac12(\Phi_{\lfloor t \rfloor_h}+\Phi_{\lceil t \rceil_h})\,. \end{align}\tag{43}\] Let \(a_t = |z_t|_{\ell^2}^2\) and \(b_t = 2 z_t \cdot w_t\), so \(b_0 = 0\) and \(\alpha a_{kh} \leq \varphi_{kh} \leq \beta a_{kh}\) since \(\alpha I \preceq H_k \preceq \beta I\). We track \(a_t, b_t\) to bound \(a_t\). From 43 , \[\label{eqn-linear-inner-loop-interpolation-ab} \begin{align} \frac{{\rm d} a_t}{{\rm d}t} &= b_t + \delta_t\,,\\ \frac{{\rm d}b_t}{{\rm d}t} &= - C_1\alpha a_t + \eta_t\,, \end{align}\tag{44}\] where \(C_1>0\) is chosen later and the perturbations are \[\begin{align} \delta_t &= 2z_t \cdot(w_{\lfloor t \rfloor_h}-w_t - (t - \lfloor t \rfloor_h)\Phi_{\lfloor t \rfloor_h}) = (t - \lfloor t \rfloor_h)z_t \cdot (\Phi_{\lceil t\rceil_h} - \Phi_{\lfloor t \rfloor_h})\,,\\ \eta_t &= C_1\alpha a_t + 2|w_{\lfloor t \rfloor_h} - (t - \lfloor t \rfloor_h)\Phi_{\lfloor t \rfloor_h}|_{\ell^2}^2 - z_t \cdot(\Phi_{\lfloor t \rfloor_h} + \Phi_{\lceil t \rceil_h})\\ &\quad -(t - \lfloor t \rfloor_h)(w_{\lfloor t \rfloor_h} - (t - \lfloor t \rfloor_h)\Phi_{\lfloor t \rfloor_h})\cdot(\Phi_{\lceil t \rceil_h}-\Phi_{\lfloor t \rfloor_h})\,. \end{align}\] The term \(\delta_t\) is piecewise smooth, with \[\frac{{\rm d} \delta_t}{{\rm d}t} = z_t \cdot(\Phi_{\lceil t \rceil_h}-\Phi_{\lfloor t\rfloor_h})+(t - \lfloor t \rfloor_h)(w_{\lfloor t \rfloor_h} - (t - \lfloor t\rfloor_h)\Phi_{\lfloor t\rfloor_h})\cdot(\Phi_{\lceil t\rceil_h}-\Phi_{\lfloor t\rfloor_h})\] in the interior of each step, and at grid points \(\delta_{kh^+} = 0\) and \[\delta_{kh^-} = h(\varphi_{kh} - \varphi_{(k-1)h}) - h(hw_{(k-1)h}-\tfrac12h^2\Phi_{(k-1)h})\cdot \Phi_{(k-1)h}\,.\] Solving the perturbed system 44 by variation of parameters, with \(s_{t} = \frac{\sin(\sqrt{C_1\alpha}t)}{\sqrt{C_1\alpha}}\) and \(c_t = \cos(\sqrt{C_1\alpha}t)\) (so \(s_t\) is increasing and nonnegative on \(0\leq t \leq \frac{\pi}{2\sqrt{C_1\alpha}}\)), gives \[\label{eqn-solution-a} \begin{align} a_t & = c_t a_0 + \int_0^t c_{t-r}\delta_r {\rm d}r + \int_0^t s_{t-r}\eta_r {\rm d}r\\ & = c_t a_0 + \sum_{k=1}^{\lceil t/h\rceil-1}[s_{t-r}\delta_r]|^{r=kh^+}_{r=kh^-} - [s_{t-r}\delta_r]|^{r=t^-}_{r=0^+} + \int_0^t s_{t-r}(\eta_r + \tfrac{{\rm d}\delta_r}{{\rm d}r}) {\rm d} r\\ & = c_t a_0 -\sum_{k=1}^{\lceil t/h\rceil -1}s_{t-kh}\delta_{kh^-}+\int_0^t s_{t-r}(-2\varphi_{\lfloor r\rfloor_h}+C_1\alpha a_{\lfloor r \rfloor_h} + \varepsilon_r) {\rm d} r\,, \end{align}\tag{45}\] where we integrated by parts and set \(\varepsilon_t \mathrel{\vcenter{:}}= \varepsilon_t^1 + \varepsilon_t^2 + \varepsilon_t^3\), \[\begin{align} \varepsilon_t^1 &= 2|w_{\lfloor t \rfloor_h} - (t - \lfloor t \rfloor_h)\Phi_{\lfloor t \rfloor_h}|_{\ell^2}^2\,,\\ \varepsilon_t^2 &= -2(z_t - z_{\lfloor t \rfloor_h})\cdot \Phi_{\lfloor t \rfloor_h}\,,\\ \varepsilon_t^3 &= C_1\alpha(z_t - z_{\lfloor t \rfloor_h}) \cdot (z_t + z_{\lfloor t \rfloor_h})\,. \end{align}\] We aim to show \(a_t \leq c_t a_0\) for a suitable \(C_1\), which requires bounding the remaining terms in 45 . Using \(|\Phi_{\lfloor t \rfloor_h}|_{\ell^2}^2 \leq \beta \varphi_{\lfloor t \rfloor_h}\), \[\varepsilon_t^1 \leq 6|w_{\lfloor t \rfloor_h}|_{\ell^2}^2 + 3(t-\lfloor t \rfloor_h)^2|\Phi_{\lfloor t \rfloor_h}|_{\ell^2}^2 \leq 6|w_{\lfloor t \rfloor_h}|_{\ell^2}^2 + 3\beta h^2\varphi_{\lfloor t \rfloor_h}\,.\] Substituting \(z_t - z_{\lfloor t \rfloor_h}\) from 43 into \(\varepsilon_t^2\), \[\begin{align} \varepsilon_t^2 &= -2 (t - \lfloor t \rfloor_h) w_{\lfloor t \rfloor_h}\cdot \Phi_{\lfloor t \rfloor_h} + (t - \lfloor t \rfloor _h)^2|\Phi_{\lfloor t\rfloor_h}|_{\ell^2}^2\\ & \leq 2|w_{\lfloor t \rfloor_h}|_{\ell^2}^2+ \tfrac32 h^2 |\Phi_{\lfloor t \rfloor_h}|_{\ell^2}^2\leq 2|w_{\lfloor t \rfloor_h}|_{\ell^2}^2 + \tfrac32\beta h^2\varphi_{\lfloor t \rfloor_h}\,, \end{align}\] and similarly for \(\varepsilon_t^3\), \[\begin{align} \varepsilon_t^3 &= C_1 \alpha|(t-\lfloor t \rfloor_h)w_{\lfloor t \rfloor_h} - \tfrac12(t - \lfloor t \rfloor_h)^2 \Phi_{\lfloor t\rfloor_h}|_{\ell^2}^2\\ & \quad +2 C_1\alpha z_{\lfloor t \rfloor_h} \cdot ((t - \lfloor t \rfloor_h) w_{\lfloor t \rfloor_h} - \tfrac12(t - \lfloor t \rfloor_h)^2\Phi_{\lfloor t\rfloor_h})\,. \end{align}\] The first term is bounded by \(2 C_1\alpha h^2|w_{\lfloor t \rfloor_h}|_{\ell^2}^2 + \tfrac{1}{2}C_1\alpha h^4\beta \varphi_{\lfloor t\rfloor_h}\). For the second, \(|z_{\lfloor t \rfloor_h}|_{\ell^2}^2 = a_{\lfloor t \rfloor_h} \leq \varphi_{\lfloor t \rfloor_h}/\alpha\) gives \[2 C_1\alpha z_{\lfloor t \rfloor_h} \cdot (t - \lfloor t \rfloor_h)w_{\lfloor t\rfloor_h} \leq h^2 C_1^2\alpha^2|z_{\lfloor t \rfloor_h}|_{\ell^2}^2 + |w_{\lfloor t \rfloor_h}|_{\ell^2}^2 \leq h^2C_1^2\alpha \varphi_{\lfloor t\rfloor_h}+|w_{\lfloor t \rfloor_h}|_{\ell^2}^2\,.\] Since \(z_{\lfloor t \rfloor_h} \cdot \Phi_{\lfloor t \rfloor_h} \geq 0\), \[\varepsilon_t^3 \leq (1 + 2C_1\alpha h^2)|w_{\lfloor t\rfloor_h}|_{\ell^2}^2 + C_1\alpha(C_1 h^2 + \tfrac12 \beta h^4)\varphi_{\lfloor t \rfloor_h}\,,\] and combining the three estimates, \[\label{eqn-bound-varepsilont} \varepsilon_t \leq (9 + 2C_1\alpha h^2)|w_{\lfloor t\rfloor_h}|_{\ell^2}^2 + \left[ \tfrac92\beta h^2+C_1\alpha(C_1 h^2 + \tfrac12 \beta h^4) \right]\varphi_{\lfloor t \rfloor_h}\,.\tag{46}\] For the boundary sum in 45 , when \(t/h\in \mathbb{Z}^+\), summation by parts gives \[\label{eqn-bound-boundary-terms-original} \begin{align} -\sum_{k=1}^{t/h-1}s_{t-kh}\delta_{kh^-} = & \sum_{k=2}^{t/h-1}h(s_{t-kh} - s_{t-(k-1)h})\varphi_{(k-1)h}-h s_h\varphi_{t-h} + h s_{t-h}\varphi_0 \\ & + \sum_{k=1}^{t/h-1}hs_{t-kh}(hw_{(k-1)h} -\tfrac12h^2\Phi_{(k-1)h})\cdot \Phi_{(k-1)h}\,. \end{align}\tag{47}\] Since \(s_{t-kh} \leq s_{t-r}\) for \(r \in [(k-1)h, kh)\), \(\varphi_{kh} \geq 0\), and \[h^2 w_{(k-1)h} \cdot \Phi_{(k-1)h} \leq \tfrac12h|w_{(k-1)h}|_{\ell^2}^2 + \tfrac12h^3|\Phi_{(k-1)h}|_{\ell^2}^2\leq \tfrac12h|w_{(k-1)h}|_{\ell^2}^2 + \tfrac12\beta h^3\varphi_{(k-1)h}\,,\] we obtain, for \(t/h \in \mathbb{Z}^+\), \[\label{eqn-bound-boundary-terms-same-momenta} -\sum_{k=1}^{t/h-1} s_{t-kh} \delta_{kh^-} \leq \int_0^h s_{t-r}\varphi_{\lfloor r \rfloor_h}{\rm d}r+ \tfrac12\int_0^t s_{t-r}(|w_{\lfloor r \rfloor_h}|_{\ell^2}^2 + \beta h^2 \varphi_{\lfloor r \rfloor_h}){\rm d} r\,.\tag{48}\] From 43 with \(w_0=0\), Cauchy–Schwarz and \(|\Phi_{\lfloor \tau\rfloor_h}|_{\ell^2}^2 \leq \beta\varphi_{\lfloor \tau\rfloor_h}\) give \[\label{eqn-w-integration-no-initial-bound} |w_{\lfloor t \rfloor_h}|_{\ell^2}^2 = \left|\tfrac12\int_0^{\lfloor t \rfloor_h}(\Phi_{\lfloor \tau \rfloor_h} + \Phi_{\lceil \tau \rceil_h}){\rm d}\tau\right|_{\ell^2}^2 \leq \tfrac12\beta \lfloor t \rfloor_h \left(\int_0^{\lfloor t\rfloor_h} \varphi_{\lfloor \tau\rfloor_h} {\rm d} \tau + \int_0^{\lfloor t \rfloor_h} \varphi_{\lceil \tau \rceil_h}{\rm d}\tau\right)\,.\tag{49}\] To bound \(\int_0^t s_{t-r}|w_{\lfloor r \rfloor_h}|_{\ell^2}^2{\rm d}r\), we note that for \(t/h\in \mathbb{Z}^+\), \[\label{eqn-sw-bound-prep} \int_0^t s_{t-r} \lfloor r \rfloor_h \left(\int_0^{\lfloor r \rfloor_h}\varphi_{\lceil \tau \rceil_h}{\rm d}\tau \right) {\rm d}r = \int_h^{t} \left(\int_{\lfloor \tau \rfloor_h}^{t}s_{t-r}\lfloor r \rfloor_h {\rm d}r\right) \varphi_{\lfloor \tau \rfloor_h}{\rm d}\tau\,.\tag{50}\] As \(s_t\) is increasing, we decompose the integral in \(r\) into sub-intervals of length \(h\), and direct calculation of the sum of \(\lfloor r \rfloor_h\) terms gives \[\int_{\lfloor\tau\rfloor_h}^{t} s_{t-r}\lfloor r\rfloor_h {\rm d}r \leq \left(\int_{\lfloor \tau \rfloor_h}^{\lfloor \tau \rfloor_h + h}s_{t-r}{\rm d}r\right) \frac{t^2}{2h}\,.\] Plugging into 50 and switching the order of integration yield an upper bound \[\int_h^t \left(\int_{\lfloor \tau \rfloor_h}^{\lfloor \tau \rfloor_h + h}s_{t-r}{\rm d}r\right) \frac{t^2}{2h} \varphi_{\lfloor \tau\rfloor_h}{\rm d}\tau = \frac{t^2}{2}\int_h^{t} s_{t-r}\varphi_{\lfloor r\rfloor_h} {\rm d}r\,.\] We can similarly derive for the other term in 49 that \[\int_0^ts_{t-r}\lfloor r \rfloor_h \left( \int_0^{\lfloor r \rfloor_h}\varphi_{\lfloor \tau \rfloor_h}{\rm d}\tau\right){\rm d}r \leq \frac{t^2}{2}\int_0^{t-h}s_{t-r}\varphi_{\lfloor r\rfloor_h}{\rm d}r\,,\] and combining these two bounds yields \[\label{eqn-sw-bound} \int_0^t s_{t-r}|w_{\lfloor r \rfloor_h}|_{\ell^2}^2 {\rm d}r \leq \tfrac12\beta t^2\int_0^t s_{t-r}\varphi_{\lfloor r \rfloor_h}{\rm d}r\,.\tag{51}\] Hence the last two terms of 45 are bounded by \[\begin{align} & -\sum_{k=1}^{\lceil t/h\rceil -1}s_{t-kh}\delta_{kh^-}+\int_0^t s_{t-r}(-2\varphi_{\lfloor r\rfloor_h}+C_1\alpha a_{\lfloor r \rfloor_h} + \varepsilon_r) {\rm d} r\\ &\leq \int_0^h s_{t-r}\varphi_{\lfloor r\rfloor_h}{\rm d}r + \int_0^t s_{t-r}(-2\varphi_{\lfloor r\rfloor_h} + C_1\alpha a_{\lfloor r\rfloor_h}) {\rm d} r\\ & \quad +\int_0^t s_{t-r}((5\beta h^2 + C_1^2\alpha h^2 + \tfrac12 C_1 \alpha \beta h^4)\varphi_{\lfloor r\rfloor_h} + (\tfrac{19}{2} + 2C_1\alpha h^2)|w_{\lfloor r\rfloor_h}|_{\ell^2}^2){\rm d}r\,. \end{align}\] Applying 51 and \(\alpha a_{kh} \leq \varphi_{kh}\), \[\begin{align} & -\sum_{k=1}^{\lceil t/h\rceil -1}s_{t-kh}\delta_{kh^-}+\int_0^t s_{t-r}(-2\varphi_{\lfloor r\rfloor_h}+C_1\alpha a_{\lfloor r \rfloor_h} + \varepsilon_r) {\rm d} r\\ \leq & \int_0^t s_{t-r}\left[-1 + C_1 + 5\beta h^2 + C_1^2\alpha h^2 + \tfrac12 C_1 \alpha \beta h^4 + (\tfrac{19}{4}+ C_1\alpha h^2)\beta t^2\right]\varphi_{\lfloor r \rfloor_h}{\rm d}r\\ \leq & \int_0^t s_{t-r}\left[-1 + C_1 + (\tfrac{39}{4}+C_1^2 + \tfrac32C_1\beta t^2)\beta t^2 \right]{\varphi_{\lfloor r \rfloor_h}}{\rm d}r\,, \end{align}\] using \(t/h\in \mathbb{Z}^+\) and \(\alpha \leq \beta\) in the last line. Taking \(C_1 = 4/9\) and \(\beta t^2 \leq 1/20\) makes the bracketed coefficient negative, so 45 gives \(a_t \leq c_t a_0\). The same choice gives \(t \leq \pi/(2\sqrt{C_1\alpha})\), ensuring \(s_t\) is monotone on the interval, and \[\begin{align} c_t &\leq 1 - \tfrac12C_1\alpha t^2 + \tfrac16C_1^2\alpha^2 t^4 \leq 1 - \tfrac12 C_1\alpha t^2 + \tfrac{1}{120}C_1^2\alpha t^2(\alpha/\beta)\\ & \leq 1 - (\tfrac12 C_1 - \tfrac{1}{120}C_1^2)\alpha t^2 \leq 1 - \tfrac15\alpha t^2\,. \end{align}\] Thus for \(t = kh\) with \(p_0 = p_0'\), \(|x_t - y_t|_{\ell^2}^2 \leq (1 - \tfrac{\alpha}{5}(kh)^2)|x_0 - y_0|_{\ell^2}^2\), and by 41 , \[|\tilde{T}_{k}^1|_{\ell^2} \leq 1 - \frac{\alpha}{10}(kh)^2\,.\]

\(\ell^2\) norm of \(h\tilde{U}_{k-1}^1\). By 41 , bounding \(|h\tilde{U}_{k-1}^1|_{\ell^2}\) again reduces to bounding \(a_t\), now with initial conditions \(z_0 = 0\), \(w_0 = p_0 - p_0'\), so \(a_0 = b_0 = 0\) and \[\label{eqn-solution-a-same-positions} \begin{align} a_t &= \int_0^t c_{t-r}\delta_r {\rm d}r + \int_0^t s_{t-r}\eta_r {\rm d}r\\ &= -\sum_{k=1}^{\lceil t/h\rceil -1}s_{t-kh}\delta_{kh^-}+\int_0^t s_{t-r}(-2\varphi_{\lfloor r\rfloor_h}+C_1\alpha a_{\lfloor r \rfloor_h} + \varepsilon_r) {\rm d} r\,. \end{align}\tag{52}\] The boundary decomposition 47 still holds, and since now \(\varphi_0=0\), it refines 48 to \[-\sum_{k=1}^{t/h-1}s_{t-kh}\delta_{kh^-} \leq \tfrac12 \int_0^t s_{t-r}(|w_{\lfloor r\rfloor_h}|_{\ell^2}^2 + \beta h^2 \varphi_{\lfloor r \rfloor_h}){\rm d}r\,.\] The bound 46 on \(\varepsilon_r\) still applies, so combining with \(\alpha a_{\lfloor r \rfloor_h} \leq \varphi_{\lfloor r \rfloor_h}\), \[\label{eqn-same-initial-pos-at-intermediate-bound} \begin{align} a_t & \leq \int_0^t s_{t-r}(-2+C_1 +5\beta h^2 + C_1^2\alpha h^2 + \tfrac12 C_1\alpha \beta h^4)\varphi_{\lfloor r \rfloor_h}{\rm d}r\\ & \quad + \int_0^t s_{t-r} (\tfrac{19}{2}+2C_1\alpha h^2)|w_{\lfloor r\rfloor_h}|_{\ell^2}^2 {\rm d}r\,. \end{align}\tag{53}\] For \(w_{\lfloor r \rfloor_h}\), repeating the derivation of 49 but with general \(w_0\), for any \(C_2 > 0\), \[\begin{align} |w_{\lfloor t \rfloor_h}|_{\ell^2}^2 &= \left| w_0 - \tfrac12 \int_0^{\lfloor t\rfloor_h}(\Phi_{\lfloor s\rfloor_h} + \Phi_{\lceil s\rceil_h}) {\rm d} s\right|_{\ell^2}^2\\ &\leq (1+C_2)|w_0|_{\ell^2}^2 + (1+\tfrac{1}{C_2}) \tfrac12\beta \lfloor t \rfloor_h \left(\int_0^{\lfloor t\rfloor_h} \varphi_{\lfloor \tau\rfloor_h} + \varphi_{\lceil \tau \rceil_h}{\rm d} \tau \right)\,. \end{align}\] With the derivation of 51 , this bounds 53 by \[\begin{align} a_t & \leq \int_0^t s_{t-r}\left[-2 + C_1 + 5\beta h^2 + C_1^2\beta h^2+ \tfrac12 C_1 \beta^2 h^4 \right.\\ &\qquad\left.+ (1 + \tfrac{1}{C_2})(\tfrac{19}{4}+C_1\beta h^2)\beta t^2 \right]\varphi_{\lfloor r \rfloor_h}{\rm d}r + (\tfrac{19}{2}+2C_1\alpha h^2)(1+C_2)\int_0^t s_{t-r}|w_0|_{\ell^2}^2 {\rm d}r\,. \end{align}\] With \(C_1 = 4/9\), \(\beta t^2 \leq 1/20\), and \(C_2 = 1/4\), the coefficient of \(\varphi_{\lfloor r \rfloor_h}\) is negative. Since \(\alpha h^2 \leq \beta t^2\) gives \(\tfrac{19}{2} + 2C_1\alpha h^2 \leq 10\), for \(t/h \in \mathbb{Z}^+\), \[a_t \leq 16 |w_0|_{\ell^2}^2 \int_0^t s_{t-r} {\rm d}r = \frac{16(1 - \cos(\sqrt{C_1\alpha}t))}{C_1\alpha} |w_0|_{\ell^2}^2 \leq 8t^2|w_0|_{\ell^2}^2\leq \frac{2}{5\beta}|w_0|_{\ell^2}^2\,,\] using \(\cos t\geq 1-\tfrac12t^2\). By 41 , \[|h\tilde{U}_{k-1}^1|_{\ell^2} \leq \sqrt{\frac{2}{5\beta}}\,.\] ◻

Proof of Proposition 9. This follows from Proposition 20, since the integrated Hessian satisfies \(\alpha I \preceq H_i \preceq \beta I\) by Assumption 1. ◻

Proof of Proposition 11. Following the decomposition of \(\nabla^2 V\) in Assumption 2, set \[H_i^{(D)} \mathrel{\vcenter{:}}= \int \nabla^2 V^{(D)} {\rm d} \nu_i\,,\qquad H_i^{(O)} \mathrel{\vcenter{:}}= \int \nabla^2 V^{(O)} {\rm d} \nu_i\,,\] so that \(H_i = H_i^{(D)}+H_i^{(O)}\). We begin with \(|\tilde{T}_k^1|_{\ell^\infty}\), which we split into a diagonal contribution and a perturbation. By the sign structure of the coefficients in Lemma 1, \[\label{eqn-propagator-infty-norm-decomp-hmc-weak} \begin{align} |\tilde{T}_k^1|_{\ell^\infty} &\leq \underbrace{|\tilde{T}_k^1(\{I - \tfrac{h^2}{2}H_i^{(D)}\}_{i=0}^{k-1})|_{\ell^\infty}}_{(a)}\\ &\quad +\underbrace{\sum_{j=1}^kh^{2j}\sum_{k-1 \geq i_1 > \dots > i_j \geq 0} \tilde{t}_{i_1,\dots,i_j}^{(k)} |H_{i_1}\cdots H_{i_j} - H_{i_1}^{(D)}\cdots H_{i_j}^{(D)}|_{\ell^\infty}}_{(b)}\,. \end{align}\tag{54}\] The matrix in \((a)\) is diagonal, so its \(\ell^\infty\) norm equals its \(\ell^2\) norm. The \(\ell^2\) bound of Proposition 9 gives \[(a) = |\tilde{T}_k^1(\{I - \tfrac{h^2}{2}H_i^{(D)}\}_{i=0}^{k-1})|_{\ell^2} \leq 1 - {\frac{\alpha}{10}(kh)^2}\,.\] For \((b)\), expanding \(H_{i_1}\cdots H_{i_j} = (H_{i_1}^{(D)}+H_{i_1}^{(O)}) \cdots (H_{i_j}^{(D)}+H_{i_j}^{(O)})\) and bounding each term yield \[\begin{align} |H_{i_1} \cdots H_{i_j} - H_{i_1}^{(D)}\cdots H_{i_j}^{(D)}|_{\ell^\infty} & \leq \sum_{l = 1}^j \binom{j}{l}(\sup_k |H_k^{(D)}|_{\ell^\infty})^{j-l}(\sup_k |H_k^{(O)}|_{\ell^\infty})^l\\ & \leq (\beta + \tfrac{\alpha}{50})^j - \beta^j\,, \end{align}\] so that \[\label{eqn-bound-b-combine-binomials} \begin{align} (b) &\leq \sum_{j=0}^k h^{2j} \sum_{k-1 \geq i_1 > \ldots > i_j \geq 0} \tilde{t}_{i_1, \ldots, i_j}^{(k)}\Big((\beta + \tfrac{\alpha}{50})^j - \beta^j\Big)\\ &= \sum_{j=0}^k h^{2j} t_j^{(k)} \Big((\beta + \tfrac{\alpha}{50})^j - \beta^j\Big) = T_k(1 + \tfrac{h^2}{2}(\beta + \tfrac{\alpha}{50})) - T_k(1 + \tfrac{h^2}{2} \beta)\,, \end{align}\tag{55}\] using that the \(t_j^{(k)}\) are the coefficients of \(T_k(1+\tfrac{1}{2}x) = \sum_{j=0}^k t_j^{(k)}x^j\), by 37 . Combining the estimates, \[\label{eqn-propagator-infty-norm-decomp-intermediate} |\tilde{T}_k^1|_{\ell^\infty} \leq 1 - {\frac{\alpha}{10}(kh)^2} + T_k(1 + \tfrac{h^2}{2} (\beta + \tfrac{\alpha}{50})) - T_k(1 + \tfrac{h^2}{2}\beta)\,.\tag{56}\] It remains to show the perturbation in \(T_k\) is controlled, so the contraction factor stays of the same order. For \(x\geq 1\), \(T_k(x)= \cosh(k\log(x + \sqrt{x^2 - 1}))\), with \[\frac{{\rm d}}{{\rm d}x} T_k(x) = \sinh(k\log(x + \sqrt{x^2 - 1})) \frac{k}{\sqrt{x^2 - 1}}\,.\] Substituting \(x = 1 + \tfrac{h^2}{2}z\) with \(z > 0\), \[\left|\frac{{\rm d}}{{\rm d} x} T_k(x)\right|\Big|_{x = 1+ \frac{h^2}{2}z} \leq \tfrac12\left(1 + \tfrac{h^2}{2}z + h \sqrt{z}\sqrt{1 + \tfrac{h^2z}{4}}\right)^k \frac{k}{h\sqrt{z}}\,,\] so by the mean value theorem, for \(z \in [\beta, \beta + \tfrac{\alpha}{50}]\), \[\label{eqn-difference-Tk-bound} \begin{align} &\left|T_k\left(1 + \tfrac{h^2}{2}\left(\beta + \tfrac{\alpha}{50}\right)\right) - T_k\left(1 + \tfrac{h^2}{2}\beta\right)\right|\leq \frac{h^2\alpha}{100} \sup_{\substack{x = 1 + \frac{h^2}{2}z,\\ z \in [\beta,\beta + \frac{\alpha}{50}]}} \left|\frac{\mathrm{d}}{\mathrm{d}x} T_k(x)\right|\\ \leq & \frac{kh\alpha}{200\sqrt{\beta}} \exp\left( kh\Big(\tfrac{51}{100}h\beta + \sqrt{\tfrac{51}{50}\beta}\sqrt{1 + \tfrac{51}{200}h^2\beta}\Big)\right)\leq {\frac{kh\alpha}{100\sqrt{\beta}}}\,, \end{align}\tag{57}\] where the last step uses \(h \leq kh \leq 1/\sqrt{20\beta}\). Substituting into 56 , \[|\tilde{T}_k^1|_{\ell^\infty} \leq 1 - {\frac{\alpha}{10}(kh)^2 + \frac{kh\alpha}{100\sqrt{\beta}}}\,.\] This remains a contraction as in the \(\ell^2\) norm for large enough \(kh\). Plugging in \(kh = \frac{1}{\sqrt{20\beta}}\) gives the upper bound \(1 - \frac{\alpha}{400\beta}\).

An analogous argument bounds \(h \tilde{U}_{k-1}^1\). Lemma 1 gives \[\label{eqn-weak-interaction-Bk-infty-norm-decomp} \begin{align} |h \tilde{U}_{k-1}^1|_{\ell^\infty} &\leq \underbrace{|h \tilde{U}_{k-1}^1(\{I - \tfrac{h^2}{2}H_i^{(D)}\}_{i=1}^{k-1})|_{\ell^\infty}}_{(a)}\\ &\quad + \underbrace{\sum_{j=1}^{k-1}h^{2j+1}\sum_{k-1\geq i_1 >\ldots >i_j \geq 1} \tilde{u}_{i_1,\ldots,i_j}^{(k)}|H_{i_1}\cdots H_{i_j} - H_{i_1}^{(D)}\cdots H_{i_j}^{(D)}|_{\ell^\infty}}_{(b)}\,, \end{align}\tag{58}\] where, as before, \((a) = |h \tilde{U}_{k-1}^1(\{I-\tfrac{h^2}{2}H_i^{(D)}\}_{i=1}^{k-1})|_{\ell^2} \leq \sqrt{\tfrac{2}{5\beta}}\) by Proposition 9, and \[\begin{align} (b) &\leq \sum_{j=0}^{k-1} h^{2j+1}\sum_{k-1\geq i_1 > \ldots > i_j \geq 1}\tilde{u}_{i_1,\ldots, i_j}^{(k)}\Big((\beta + \tfrac{\alpha}{50})^j - \beta^j \Big)\\ & = h U_{k-1}(1 + \tfrac{h^2}{2}(\beta + \tfrac{\alpha}{50})) - h U_{k-1}(1 + \tfrac{h^2}{2}\beta)\leq h U_{k-1}(1 + \tfrac{h^2}{2}(\beta + \tfrac\alpha{50}))\,, \end{align}\] since \(U_{k-1}(x) \geq 0\) for \(x \geq 1\). Using \(U_{k-1}(x) = \sinh(k\log(x+\sqrt{x^2-1}))/\sqrt{x^2-1}\) for \(x \geq 1\), \[\label{eqn-Chebyshev-U-b-bound} (b) \leq \frac{1}{2\sqrt\beta}\exp\left( kh\Big(\tfrac{51}{100}h\beta + \sqrt{\tfrac{51}{50}\beta}\sqrt{1 + \tfrac{51}{200}h^2\beta}\Big)\right)\leq \frac{13}{20\sqrt\beta}\,.\tag{59}\] Combining these in 58 , \[|h\tilde{U}_{k-1}^1|_{\ell^\infty} \leq \frac{3}{2\sqrt{\beta}}\,.\] ◻

Proof of Lemma 2. The \((i,j)\) entry of the product is nonzero only if \(j \in \mathsf{N}_r(i)\), or if \(i \in \mathsf{N}_r(k)\) and \(j \in \mathsf{N}_r(k)\) for some \(k\). Since \(\max_{1\leq i \leq d} |\mathsf{N}_{2r}(i)| = s_r\), each row has at most \(s_r\) nonzero entries, so the matrix is \(s_r\)-sparse. ◻

Proof of Proposition 21. The bounds follow a common framework, so we give the details for \(A_{Nm+k}\). The strategy is to split the matrix expansion by order in \(h\) at a threshold \(r\) that balances sparsity against magnitude, as for a tridiagonal matrix, low-order powers stay sparse (so their \(\ell^\infty\) and \(\ell^2\) norms are comparable), while the growing density of higher-order powers is offset by the decay in \(h^{2j}\). The sparsity here is controlled by \(s_r\).

For any threshold \(r\), Lemma 2 gives that the entries of each order-\(r\) product are supported on \(\mathsf{N}_{2r}(i)\), and the absolute row sum decomposes as \[\label{eqn-inf-bound-two-parts-a} \begin{align} \sum_{j=1}^d|(A_{Nm+k})_{ij}| = & \sum_{j\in \mathsf{N}_{2r}(i)}|(A_{Nm+k})_{ij}| + \sum_{j\notin \mathsf{N}_{2r}(i)}|(A_{Nm+k})_{ij}|\\ \leq &\sqrt{|\mathsf{N}_{2r}(i)|}\sqrt{\sum_{j\in\mathsf{N}_{2r}(i)}|(A_{Nm+k})_{ij}|^2} + \sum_{j \notin \mathsf{N}_{2r}(i)}|(A_{Nm+k})_{ij}|\\ \leq & \sqrt{s_{r}}|A_{Nm+k}|_{\ell^2} + \sum_{j \notin \mathsf{N}_{2r}(i)}|(A_{Nm+k})_{ij}|\,. \end{align}\tag{60}\] By Proposition 19, \[\label{eqn-decomp-neighbor-matrix-A} \begin{align} A_{Nm+k} &= a_0^{(Nm+k)} I + \sum_{j=1}^{r} (-1)^jh^{2j}\sum_{i_1>\dots >i_j}a_{i_1,\dots,i_j}^{(Nm+k)} H_{i_1}\cdots H_{i_j}\\ & \quad + \sum_{j=r+1}^{Nm+k}(-1)^jh^{2j}\sum_{i_1>\dots>i_j}a_{i_1,\dots,i_j}^{(Nm+k)}H_{i_1}\cdots H_{i_j}\,. \end{align}\tag{61}\] For \(j \notin \mathsf{N}_{2r}(i)\), the entry \((A_{Nm+k})_{ij}\) can only come from the second line, so \[\label{eqn-far-neighbor-elements-bound} \sum_{j \notin \mathsf{N}_{2r}(i)}|(A_{Nm+k})_{ij}| \leq \Big|\sum_{j=r+1}^{Nm+k}(-1)^j h^{2j}\sum_{i_1,\dots,i_j} a_{i_1,\dots,i_j}^{(Nm+k)}H_{i_1}\dots H_{i_j}\Big|_{\ell^\infty}\,.\tag{62}\] We control this \(\ell^\infty\) norm through the \(\ell^2\) norm, \[\label{eqn-tail-matrices} \Big|\sum_{j=r+1}^{Nm+k}(-1)^j h^{2j}\sum_{i_1,\dots,i_j} a_{i_1,\dots,i_j}^{(Nm+k)}H_{i_1}\dots H_{i_j}\Big|_{\ell^2} \leq \sum_{j=r+1}^{Nm+k}h^{2j}\beta^j a_j^{(Nm+k)}\,.\tag{63}\] By Proposition 19, \(a_j^{(Nm+k)}\) is bounded by \(\max_{j\leq i\leq Nm+k} t_j^{(i)}\), the Chebyshev coefficients whose explicit formulas are given in Proposition 17. It therefore suffices to bound \(h^{2j}\beta^j t_j^{(i)}\). For \(i > 0\), \[\begin{align} t_j^{(i)} &= i\frac{(i+j-1)!}{(i-j)!(2j)!} = \frac{i(i+j-1) \cdots (i+1)}{(2j)\cdots(j+1)} \frac{i(i-1)\cdots(i-j+1)}{j!}\\ & \leq \frac{(i+j)^j}{j^j} \frac{i^j}{j!} \leq \frac{(i+j)^j (ie)^j}{j^{2j}} = \left(e\left(\tfrac{i^2}{j^2}+\tfrac{i}{j} \right) \right)^j\,, \end{align}\] using \(j! \geq j^j/\exp(j)\) (from \(\exp(j) \geq j^j/j!\)). Hence for \(j \geq (Nm+k)\sqrt{\beta}eh\) and \(i \leq Nm+k\), \[h^{2j}\beta^j t_j^{(i)}\leq \left( h^2\beta e\left(\tfrac{i^2}{j^2}+\tfrac{i}{j} \right) \right)^j \leq \left(e^{-1} + \sqrt\beta h\right)^j \leq \left(\tfrac{3}{5}\right)^j\,,\] where we use \(h \leq mh \leq 1/\sqrt{20\beta}\). Since this holds uniformly over \(i \leq Nm+k\), Proposition 19 gives \(h^{2j} \beta^j a_j^{(Nm+k)} \leq \left(\tfrac35\right)^j\), so 63 yields \[\Big|\sum_{j=r+1}^{Nm+k}(-1)^j h^{2j}\sum_{i_1,\dots,i_j} a_{i_1,\dots,i_j}^{(Nm+k)}H_{i_1}\dots H_{i_j}\Big|_{\ell^2} \leq \sum_{j=r+1}^{Nm+k}\left(\tfrac35\right)^j \leq {3} \left(\tfrac35\right)^{r+1}\,.\] Applying the norm inequality \(|M|_{\ell^\infty} \leq \sqrt{d}\,|M|_{\ell^2}\), valid for any matrix \(M \in \mathbb{R}^{d\times d}\), to 62 , \[\sum_{j \notin \mathsf{N}_{2r}(i)}|(A_{Nm+k})_{ij}| \leq {3}\sqrt{d}\left( \tfrac35\right)^{r+1}\,.\] Substituting into 60 with \(|A_{Nm+k}|_{\ell^2} \leq 1\) (Proposition 9) and taking \(r_{Nm+k} = \lceil (Nm+k)h\sqrt\beta e + \tfrac{\log \sqrt d}{\log (5/3)} \rceil\), \[|A_{Nm+k}|_{\ell^\infty} \leq \sqrt{s_{r_{Nm+k}}} + {3} \left( \tfrac35\right)^{(Nm+k)h\sqrt\beta e}\,,\] which gives the first bound in ?? since \(s_{r_{Nm+k}} \geq 1\).

For \(A_{Nm+k}H\) the decomposition is the same as 60 with \(H\) appended, giving \(|A_{Nm+k}H|_{\ell^2} \leq \beta |A_{Nm+k}|_{\ell^2}\) and \[\Big|\sum_{j=r}^{Nm+k}(-1)^jh^{2j}\sum_{i_1>\dots>i_j}a_{i_1,\dots,i_j}^{(Nm+k)}H_{i_1}\cdots H_{i_j} H\Big|_{\ell^2} \leq \beta \sum_{j=r}^{Nm+k}\left(\tfrac35\right)^j \leq {3}\beta\left( \tfrac{3}{5}\right)^r\,,\] so taking \(r_{Nm+k}\) as before, \[|A_{Nm+k} H|_{\ell^\infty} \leq \beta \sqrt{s_{r_{Nm+k}}} + {3}\beta\left(\tfrac35\right)^{(Nm+k)h\sqrt\beta e} \leq 4\beta \sqrt{s_{r_{Nm+k}}}\,.\]

For \(B_{Nm+k-1}\), the same decomposition gives, for any \(r\), \[\label{eqn-inf-bound-two-parts-b} \sum_{j=1}^d|(B_{Nm+k-1})_{ij}| \leq \sqrt{s_r}|B_{Nm+k-1}|_{\ell^2} + \sqrt{d}\sum_{j=r+1}^{Nm+k-1}h^{2j+1}\beta^jb_j^{(Nm+k-1)}\,.\tag{64}\] The difference from \(A_{Nm+k}\) is that \(b_j^{(Nm+k-1)}\) depends on \(u_j^{(i)}\). By Proposition 17, \[u_j^{(i)} = \frac{(i+j+1)!}{(i-j)!(2j+1)!} \leq \frac{(i+j+1)^j}{j^j}\frac{(i+1)^{j+1}e^{j+1}}{j^{j+1}}\,,\] so by Proposition 19, for \(j \geq (Nm+k)\sqrt\beta eh\) and \(i \leq Nm+k-1\), \[h^{2j+1}\beta^jb_j^{(Nm+k-1)} \leq \frac{1}{\sqrt\beta}\left(h\sqrt\beta\left(\tfrac{i+1}{j}+1\right) \right)^j\left(\tfrac{h\sqrt\beta(i+1)e}{j}\right)^{j+1} \leq \frac{1}{\sqrt\beta}\left(\tfrac35\right)^j\,,\] and hence \[\sqrt{d}\sum_{j=r+1}^{Nm+k-1}h^{2j+1}\beta^jb_j^{(Nm+k-1)} \leq \frac{{3}\sqrt{d}}{\sqrt{\beta}}\left(\tfrac35\right)^{r+1}\,.\] Taking \(r = r_{Nm+k}\) and using \(|B_{Nm+k-1}|_{\ell^2} \leq \sqrt{\tfrac{2}{5\beta}}\) (Proposition 9), 64 gives \[|B_{Nm+k-1}|_{\ell^\infty} \leq \sqrt{\tfrac{2}{5\beta}}\sqrt{s_{r_{Nm+k}}}+ \frac{{3}}{\sqrt\beta}\left(\tfrac35\right)^{(Nm+k)h\sqrt\beta e} \leq \frac{4}{\sqrt\beta} \sqrt{s_{r_{Nm+k}}}\,.\] The bounds for \(B_{Nm+k-1}H\) and \(B_{Nm+k-1}HH'\) follow from similar decompositions: \[|B_{Nm+k-1}H|_{\ell^\infty} \leq \sqrt{\tfrac{2\beta}{5}}\sqrt{s_{r_{Nm+k}}} + {3}\sqrt\beta \left(\tfrac35\right)^{(Nm+k)h\sqrt\beta e} \leq 4 \sqrt\beta\sqrt{s_{r_{Nm+k}}}\,,\] \[|B_{Nm+k-1}H H'|_{\ell^\infty} \leq \sqrt{\tfrac{2\beta^3}{5}}\sqrt{s_{r_{Nm+k}}} + \frac{25}{6} \beta^{\frac{3}{2}}\left(\tfrac{3}{5} \right)^{(Nm+k)\sqrt\beta eh} \leq 5\beta^{\frac{3}{2}}\sqrt{s_{r_{Nm+k}}}\,.\] ◻

Proof of Proposition 22. Consider the time-inhomogeneous linear dynamics driven by \(\{H_i\}\), \[\label{eqn-baoab-linear} \begin{align} q_{k+1} &= q_k + \tfrac{h}{2}(1+\eta)p_k - \tfrac{h^2}{4}(1+\eta)H_k q_k+ \tfrac{h}{2}\sqrt{1-\eta^2}\xi_{k}\,,\\ p_{k+1} &= \eta(p_k - \tfrac{h}{2}H_kq_k) + \sqrt{1-\eta^2}\xi_{k}-\tfrac{h}{2}H_{k+1}q_{k+1}\,, \end{align}\tag{65}\] which is the BAOAB scheme 4 with the varying potential \(V_k(q) = \tfrac12q^\top H_k q\) at step \(k\). For two trajectories coupled through the same noise \(\{\xi_i\}\), the differences \(\Delta q_k, \Delta p_k\) obey \[\label{eqn-leap-frog-varying-potential} \begin{align} \Delta q_{k+1} &= (I - \tfrac{h^2}{4}(1+\eta)H_k) \Delta q_k + \tfrac{h}{2}(1+\eta)\Delta p_k\,,\\ \Delta p_{k+1} &= \eta(\Delta p_k - \tfrac{h}{2}H_k\Delta q_k) -\tfrac{h}{2}H_{k+1}\Delta q_{k+1}\,. \end{align}\tag{66}\] These yield the same three-term recurrence \(\Delta q_{k+1} = (1+\eta)(I - \tfrac{h^2}{2}H_k)\Delta q_k - \eta \Delta q_{k-1}\) as in the proof of Proposition 8, so the same induction gives \[\label{eqn-q-on-initial-q0-p0-baoab-linear} \Delta q_{k} = \tilde{T}_{k}^\eta \Delta q_0 + h\tilde{U}_{k-1}^\eta \Delta p_0\,,\tag{67}\] identifying \(\tilde{T}_k^\eta\) and \(h\tilde{U}_{k-1}^\eta\) as the position-difference propagators of 65 .

Although Lemma 3 is stated for a fixed potential, its proof uses only the spectral bounds \(\alpha I \preceq H_i \preceq \beta I\). [29] employs the decomposition \[\label{eqn-baoab-decomposition-abao} (\mathsf{B}\mathsf{A}\mathsf{O}\mathsf{A}\mathsf{B})^k = \mathsf{B}\mathsf{A}\mathsf{O}(\mathsf{A}\mathsf{B}\mathsf{A}\mathsf{O})^{k-1}\mathsf{A}\mathsf{B}\tag{68}\] and shows that the core operator \(\mathsf{A}\mathsf{B}\mathsf{A}\mathsf{O}\) is a weighted \(\ell^2\) contraction (giving the \((1-c(h))\) factors) while the boundary terms \(\mathsf{B}\mathsf{A}\mathsf{O}\) and \(\mathsf{A}\mathsf{B}\) stay bounded, all under those eigenvalue bounds. Since each \(H_i\) satisfies them, the result carries over to the varying-potential case, giving under 65 \[\label{eqn-baoab-weighted-2-norm-varying-potential} |(\Delta q_k, \Delta p_k)|_{\ell^2_w} \leq 7(1-c(h))^{\frac{k-1}{2}}|(\Delta q_0, \Delta p_0)|_{\ell^2_w}\,.\tag{69}\]

We now derive the operator-norm bounds. Since \(M_{\text{\normalfont\scshape baoab},k}^w = W M_{\text{\normalfont\scshape baoab},k} W^{-1}\) maps \(W(\Delta q_0, \Delta p_0)\) to \(W(\Delta q_k, \Delta p_k)\), 69 is precisely ?? . For the position blocks, the norm equivalence 16 gives \[|\Delta q_k|_{\ell^2} \leq \sqrt{|\Delta q_k|_{\ell^2}^2 + a |\Delta p_k|_{\ell^2}^2} \leq \sqrt{2}|(\Delta q_k, \Delta p_k)|_{\ell^2_w}\,,\] so \[\label{eqn-useful-cor-baoab-pre-bound} |\Delta q_k|_{\ell^2} \leq 7\sqrt{2}(1- c(h))^{\frac{k-1}{2}}|W(\Delta q_0, \Delta p_0)|_{\ell^2}\,.\tag{70}\] Writing 67 as \(\Delta q_k = \begin{bmatrix} \tilde{T}_{k}^\eta & h\tilde{U}_{k-1}^\eta \end{bmatrix} W^{-1}\, W(\Delta q_0, \Delta p_0)\), the bound 70 gives \[\left|\begin{bmatrix} \tilde{T}_{k}^\eta & h\tilde{U}_{k-1}^\eta \end{bmatrix} W^{-1} \right|_{\ell^2} \leq 7\sqrt{2}(1-c(h))^{\frac{k-1}{2}}\,,\] and since \[\begin{bmatrix} \tilde{T}_k^\eta & h \tilde{U}_{k-1}^\eta \end{bmatrix} W^{-1} = \begin{bmatrix} \tilde{T}_k^\eta & - \frac{b}{\sqrt{a-b^2}} \tilde{T}_k^\eta + \frac{1}{\sqrt{a-b^2}}h\tilde{U}_{k-1}^\eta \end{bmatrix}\,,\] the block-wise bounds ?? follow. ◻

Proof of Corollary 1. Inequality 71 follows from 69 and the matrix norm equivalence \(|M|_{\ell^\infty} \leq \sqrt{2d}|M|_{\ell^2}\) for \(M \in \mathbb{R}^{2d \times 2d}\). For 72 , 70 gives \[|\Delta q_k|_{\ell^\infty} \leq |\Delta q_k|_{\ell^2} \leq 7\sqrt{2}(1-c(h))^{\frac{k-1}{2}}|(\Delta q_0, \Delta p_0)|_{\ell^2_w} \leq 14\sqrt{d}(1-c(h))^{\frac{k-1}{2}}|(\Delta q_0, \Delta p_0)|_{\ell^\infty_w}\,.\] ◻

Proof of Proposition 10. This follows from Proposition 22, since the integrated Hessians satisfy the assumption. ◻

Proof of Proposition 12. Estimate for \(M_{\text{\normalfont\scshape abao}}^w\). The matrix is \[M_{\text{\normalfont\scshape abao}}^w(H_0) = \begin{bmatrix} I - \tfrac12 h^2 H_0 - \tfrac{\eta}{1-\eta}h^2 H_0 & h^3\tfrac{\frac{1}{4} + \frac{\eta}{(1-\eta)^2}}{\sqrt{a-b^2}}H_0\\ - \eta h \sqrt{a-b^2}H_0 & \eta (I + (-\tfrac12 + \tfrac{1}{1-\eta})h^2H_0) \end{bmatrix}\,.\] For the second block of rows, the \(\ell^1\) norm of each row is at most \[\label{eqn-abao-row-2-bound} \eta h \sqrt{a-b^2} \cdot \tfrac{21}{20}\beta + \eta (1 + (-\tfrac12 + \tfrac{1}{1-\eta}) h^2 \cdot \tfrac{21}{20}\beta)\,.\tag{73}\] This increases in \(\eta \in [0,1]\), while the target bound \(1 - \frac{h^2\alpha}{8(1-\eta)}\) decreases in \(\eta\). Comparing them at the largest admissible value \(\eta = 1-2\sqrt\beta h\) confirms the target bound for the second block.

For the first block of rows, split off the diagonal and off-diagonal parts \(H_0^{(D)}, H_0^{(O)}\) (defined as in the proof of Proposition 11): \[\begin{align} &\begin{bmatrix} I - \tfrac12 h^2 H_0^{(D)} - \tfrac{\eta}{1-\eta}h^2 H_0^{(D)} & h^3\tfrac{\frac{1}{4} + \frac{\eta}{(1-\eta)^2}}{\sqrt{a-b^2}} H_0^{(D)} \end{bmatrix}\\ & \quad + \begin{bmatrix} - \tfrac12 h^2 H_0^{(O)} - \tfrac{\eta}{1-\eta}h^2 H_0^{(O)} & h^3\tfrac{\frac{1}{4} + \frac{\eta}{(1-\eta)^2}}{\sqrt{a-b^2}} H_0^{(O)} \end{bmatrix}\,. \end{align}\] In the first matrix, \(H_0^{(D)}\) is diagonal, so each row reduces to a scalar \(\lambda \in [\alpha, \beta]\). Since \(\tfrac12 h^2\beta + \tfrac{\eta}{1-\eta}h^2\beta\) increases in \(\eta\) and is below \(1\) at \(\eta = 1 - 2\sqrt{\beta} h\), it is below \(1\) throughout, and the \(\ell^1\) norm of the row is \[\label{eqn-abao-row-1-bound} 1 - \tfrac12 h^2\lambda - \tfrac{\eta}{1-\eta}h^2\lambda + h^3 \tfrac{\frac{1}{4} + \frac{\eta}{(1-\eta)^2}}{\sqrt{a-b^2}}\lambda\,.\tag{74}\] The condition \(h \leq \frac{1-\eta}{2\sqrt\beta}\) gives \(\sqrt{a-b^2} \geq \frac{\sqrt{3}}{2 \sqrt{\beta}}\), so the last term is at most \(\frac{1}{4\sqrt{3}}h^2 \lambda + \frac{1}{\sqrt{3}}\frac{\eta}{1-\eta}h^2 \lambda\), and 74 is bounded by \(1 - c(h)\). Adding the off-diagonal contribution, the \(\ell^\infty\) norm of the first block of rows is \[1 - c(h) + \left((\tfrac12 + \tfrac{\eta}{1-\eta})h^2 + h^3\tfrac{\frac{1}{4} + \frac{\eta}{(1-\eta)^2}}{\sqrt{a-b^2}}\right)| H_0^{(O)}|_{\ell^\infty} \leq 1 - \tfrac12c(h)\,.\]

Estimate for \(M_{\text{\normalfont\scshape bao}}^w\). By the scheme definition in 14 , \[M_{\text{\normalfont\scshape bao}}^w(H_0) = \begin{bmatrix} I - \tfrac{h^2}{4}H_0 - \tfrac{h^2\eta}{2(1-\eta)}H_0 & \tfrac{1}{\sqrt{a-b^2}}(-\tfrac{h}{2}I + \tfrac{h^3(1+\eta)}{4(1-\eta)^2}H_0)\\ - \tfrac{h}{2}\eta \sqrt{a-b^2}H_0 & \eta I + \tfrac{h^2\eta}{2(1-\eta)}H_0 \end{bmatrix}\,.\] We bound each \(2\times2\) block in \(\ell^\infty\). The \((1,1)\) block is bounded by \(1 + (\tfrac{h^2}{4} + \tfrac{h}{2} \cdot \tfrac{1}{2\sqrt\beta})\tfrac{21}{20}\beta \leq \tfrac{383}{320}\), the \((1,2)\) block by \(\frac{2\sqrt{\beta}}{\sqrt{3}}(\tfrac{h}{2} + \tfrac{h}{4}\cdot 2 \cdot \tfrac{1}{4\beta} \cdot \tfrac{21}{20}\beta) \leq \frac{101}{160\sqrt{3}}\), the \((2,1)\) block by \(\tfrac12 \tfrac{1}{\sqrt{\beta}} h \cdot \tfrac{21}{20}\beta \leq \tfrac{21}{80}\), and the \((2,2)\) block by \(\tfrac12 \cdot \tfrac{1}{2\sqrt\beta}h \cdot\tfrac{21}{20}\beta + 1 \leq \tfrac{181}{160}\). Summing the block bounds over each row gives the result.

Estimate for \(M_{\text{\normalfont\scshape ab}}^w\). By 14 , \[M_{\text{\normalfont\scshape ab}}^w(H_0) = \begin{bmatrix} I - \tfrac{h^2}{2(1-\eta)}H_0 & \tfrac{1}{\sqrt{a-b^2}}(\tfrac{h}{2}I + \tfrac{h^3(1+\eta)}{4(1-\eta)^2}H_0)\\ - \tfrac{h}{2}\sqrt{a-b^2}H_0 & I + (\tfrac{h^2}{2(1-\eta)}-\tfrac{h^2}{4})H_0 \end{bmatrix}\,.\] By the same block-wise bounding, the \((1,1)\) block is at most \(1 + \tfrac{h}{4\sqrt\beta}\tfrac{21}{20}\beta \leq \tfrac{181}{160}\), the \((1,2)\) block at most \(\frac{2\sqrt{\beta}}{\sqrt{3}}(\tfrac{h}{2}+\tfrac{h}{2}\tfrac{1}{4\beta}\tfrac{21}{20}\beta) \leq \frac{101}{160\sqrt{3}}\), the \((2,1)\) block at most \(\tfrac{h}{2\sqrt{\beta}}\tfrac{21}{20}\beta \leq \tfrac{21}{80}\), and the \((2,2)\) block at most \(1 + \tfrac{h}{2}\tfrac{1}{2\sqrt{\beta}}\tfrac{21}{20}\beta \leq \tfrac{181}{160}\). A row sum concludes the proof. ◻

Proof of Proposition 14. Write \(A_k \mathrel{\vcenter{:}}= \tilde{T}_{k}^\eta\) and \(B_{k-1}\mathrel{\vcenter{:}}= h\tilde{U}_{k-1}^\eta\). As in 60 , for each threshold \(r\), \[\label{eqn-linf-sparse-baoab-decomp} \sum_{j=1}^d|(A_k)_{ij}| \leq \sqrt{s_r}|A_k|_{\ell^2} + \sum_{j \notin \mathsf{N}_{2r}(i)}|(A_k)_{ij}|\,,\tag{75}\] where \(|A_k|_{\ell^2} \leq 7\sqrt{2}(1-c(h))^{\frac{k-1}{2}}\) by Proposition 10, and the far-neighbor sum comes from the terms of order greater than \(r\), \[\sum_{j \notin \mathsf{N}_{2r}(i)}|(A_k)_{ij}| \leq \Big|\sum_{j=r+1}^k(-1)^jh^{2j}\sum_{k-1 \geq i_1 > \dots > i_j \geq 0} \tilde{t}_{i_1, \dots, i_j}^{(k)}H_{i_1}\dots H_{i_j}\Big|_{\ell^\infty}\,,\] whose \(\ell^2\) norm is at most \(\sum_{j=r+1}^k h^{2j} \beta^j t_j^{(k)}\). By Proposition 18, \[t_j^{(k)} \leq \left(\tfrac{1}{1-\eta}\right)^j \binom{k}{j} \leq \left( \tfrac{k e}{(1-\eta)j} \right)^j\,,\] using \(j! \geq j^j/\exp(j)\). Hence for \(r > \frac{ke^2\beta}{1-\eta}h^2\), \[\begin{align} \Big|\sum_{j=r+1}^k(-1)^jh^{2j}\sum_{k-1 \geq i_1 > \dots > i_j \geq 0} \tilde{t}_{i_1, \dots, i_j}^{(k)}H_{i_1}\dots H_{i_j}\Big|_{\ell^2} &\leq \sum_{j=r+1}^k \left( h^2 \tfrac{ke}{(1-\eta)j}\beta\right)^j\\ &\leq \sum_{j=r+1}^k \exp(-j) \leq \exp(-r)\,. \end{align}\] Substituting into 75 with \(|M|_{\ell^\infty} \leq \sqrt d |M|_{\ell^2}\) and taking \(r_k \mathrel{\vcenter{:}}= \lceil \frac{ke^2\beta}{1-\eta}h^2+\log\sqrt{d} \rceil\), \[\begin{align} |A_k|_{\ell^\infty} &\leq 7\sqrt2\,\sqrt{s_{r_k}}(1- c(h))^{\frac{k-1}{2}}+ \exp\left(-\tfrac{ke^2\beta}{1-\eta}h^2\right) \\&\leq (7\sqrt 2 + 1) \sqrt{s_{r_k}} \exp\left(-\tfrac{\alpha h^2}{8(1-\eta)}(k-1)\right)\,. \end{align}\] For \(B_{k-1}\), Proposition 18 gives \[u_j^{(k-1)} \leq \left(\tfrac{1}{1-\eta}\right)^{j+1}\binom{k-1}{j} \leq \left( \tfrac{1}{1-\eta}\right)^{j+1} \left(\tfrac{(k-1)e}{j}\right)^j\,,\] so the terms of order greater than \(r\) have \(\ell^2\) norm at most \(\sum_{j=r+1}^{k-1}h^{2j+1}\beta^ju_j^{(k-1)} \leq \frac{h}{1-\eta} \exp(-r) = b \exp(-r)\). Combining the tail estimates for \(A_k\) and \(B_{k-1}\), \[\sum_{j\notin \mathsf{N}_{2r}(i)}\left| \left(-\tfrac{b}{\sqrt{a-b^2}} A_k + \tfrac{1}{\sqrt{a-b^2}}B_{k-1} \right)_{ij}\right| \leq 2\sqrt{d}\, \tfrac{b}{\sqrt{a-b^2}}\exp(-r)\,,\] while Proposition 10 gives \(|-\tfrac{b}{\sqrt{a-b^2}} A_k + \tfrac{1}{\sqrt{a-b^2}}B_{k-1}|_{\ell^2} \leq 7\sqrt{2}(1-c(h))^{\frac{k-1}{2}}\). Hence, simplifying \(\frac{b}{\sqrt{a-b^2}}\) via \(h\leq \frac{1-\eta}{2\sqrt{\beta}}\), \[\begin{align} &\left|-\tfrac{b}{\sqrt{a-b^2}} A_k + \tfrac{1}{\sqrt{a-b^2}}B_{k-1}\right|_{\ell^\infty} \leq \sqrt{s_{r_k}}7\sqrt{2}(1-c(h))^{\frac{k-1}{2}} + \tfrac{2b}{\sqrt{a-b^2}}\exp\left(-\tfrac{ke^2\beta}{1-\eta}h^2\right)\\ &\qquad \leq \left(7\sqrt{2}+\tfrac{2}{\sqrt{3}}\right)\sqrt{s_{r_k}}\exp\left(-\tfrac{\alpha h^2}{8(1-\eta)}(k-1)\right)\,. \end{align}\] ◻

Proof of Proposition 23. Let \((q_t, p_t) = \mathsf{U}_{\text{\normalfont\scshape hmc}}^t(q^*, p^*)\). The continuous and discrete position updates have the integral representations \[\label{eqn-integral-hmc-position} \mathsf{Q}^h(q^*,p^*) = q^* + p^*h - \int_0^h\int_0^t \nabla V(q_s){\rm d} s {\rm d} t\,,\quad \mathsf{Q}_h^1(q^*,p^*) = q^* + p^*h - \tfrac12 h^2 \nabla V(q^*)\,,\tag{78}\] whose difference is \(\mathsf{Q}_h^1(q^*,p^*)-\mathsf{Q}^h(q^*,p^*) = \int_0^h \int_0^t (\nabla V(q_s)-\nabla V(q^*)) {\rm d}s {\rm d}t\). Writing \(\nabla V(q_s) - \nabla V(q^*) = (\int_0^1 \nabla^2 V(\tau q_s + (1-\tau)q^*){\rm d}\tau)(q_s - q^*)\) and using the hypothesis on \(M\), \[|M(\mathsf{Q}_h^1(q^*,p^*)-\mathsf{Q}^h(q^*,p^*))|_{\ell^{\infty}}^2 \leq \tfrac12 h^2 C_3^2 \int_0^h\int_0^t|q_s-q^*|_{\ell^\infty}^2{\rm d}s {\rm d}t\,.\] Expanding \(q_s\) bounds the remaining integral, \[\label{eqn-qs-q-double-integral} \begin{align} \int_0^h\int_0^t|q_s-q^*|_{\ell^\infty}^2 {\rm d}s{\rm d}t &\leq 2\int_0^h\int_0^t|p^*s|_{\ell^\infty}^2 {\rm d}s{\rm d}t + 2 \int_0^h\int_0^t\Big|\int_0^s\int_0^\tau \nabla V(q_r){\rm d}r {\rm d}\tau\Big|_{\ell^\infty}^2{\rm d}s{\rm d}t\\ &\leq \tfrac16h^4|p^*|_{\ell^\infty}^2 + \int_0^h\int_0^ts^2\int_0^s\int_0^\tau|\nabla V(q_r)|_{\ell^\infty}^2{\rm d}r{\rm d}\tau {\rm d}s{\rm d}t\,. \end{align}\tag{79}\] Since \((q^*,p^*) \sim \pi \otimes \mathcal{N}(0,I)\), stationarity gives \(q_t \sim \pi\) for all \(t\), so by Lemma 4, \[\label{eqn-qs-q-double-integral-expectation} \mathbb{E}\int_0^h\int_0^t|q_s-q^*|_{\ell^\infty}^2{\rm d}s{\rm d}t\leq \tfrac23h^4\log(2d) + \tfrac{1}{60}h^6 \mathbb{E}_\pi[|\nabla V|_{\ell^\infty}^2]\,.\tag{80}\] Combining these with Lemma 5 and \(h^2\beta \leq 1/20\), \[\begin{align} |M(\mathsf{Q}_h^1(q^*,p^*)-\mathsf{Q}^h(q^*,p^*))|_{2,\ell^\infty} &\leq \tfrac1{\sqrt{3}} h^3C_3\sqrt{\log(2d)} + \tfrac{1}{\sqrt{120}}h^4C_3\sqrt{\mathbb{E}_{\pi}[|\nabla V|_{\ell^\infty}^2]}\\ &\leq \tfrac1{\sqrt{2}}h^3C_3\sqrt{\log(2d)}\,, \end{align}\] which is ?? . For the momentum component, \[\mathsf{P}_h^1(q^*,p^*) = p^*-\tfrac12h\nabla V(q^*)-\tfrac12h\nabla V(\mathsf{Q}_h^1(q^*,p^*))\,,\quad \mathsf{P}^h(q^*,p^*) = p^*-\int_0^h\nabla V(q_t) {\rm d}t\,,\] so \[\label{eqn-decomp-MPhP-a-b} \begin{align} &|M(\mathsf{P}_h^1(q^*,p^*)-\mathsf{P}^h(q^*,p^*))|_{2,\ell^\infty}^2 \\ &\leq \tfrac12 h C_3^2 \underbrace{\mathbb{E}\int_0^h|q_t-q^*|_{\ell^\infty}^2 {\rm d}t}_{(a)} + \tfrac12\underbrace{\mathbb{E}\Big|\int_0^h \tilde{M}(q_t - \mathsf{Q}_h^1(q^*,p^*)){\rm d}t\Big|_{\ell^\infty}^2}_{(b)}\,, \end{align}\tag{81}\] where \(\tilde{M} \mathrel{\vcenter{:}}= \int_0^1 M \nabla^2 V(\tau q_t + (1-\tau)\mathsf{Q}_h^1(q^*,p^*)){\rm d}\tau\). The bounds used in 79 and 80 give \[\label{eqn-bound-a} (a) \leq \tfrac83h^3\log(2d)+ \tfrac{1}{10} h^5 \mathbb{E}_\pi[|\nabla V|_{\ell^\infty}^2]\,.\tag{82}\] Using the position differences in 78 , \[\begin{align} &q_t - \mathsf{Q}_h^1(q^*,p^*) = p^*(t-h) - \int_0^t \int_0^s (\nabla V(q_\tau)-\nabla V(q^*)){\rm d}\tau {\rm d}s+\tfrac12(h^2-t^2)\nabla V(q^*)\\ = & p^*(t-h) - \int_0^t \int_0^s \Big(\int_0^1\nabla^2V(rq_\tau + (1-r)q^*) {\rm d}r \Big) (q_\tau - q^*){\rm d}\tau {\rm d}s+\tfrac12(h^2-t^2)\nabla V(q^*)\,, \end{align}\] so \((b)\) splits as \[\label{eqn-b-decomp-b1-b2-b3} \begin{align} (b) &\leq 3h C_3^2\underbrace{\mathbb{E}\int_0^h |p^*(t-h)|_{\ell^\infty}^2 {\rm d}t}_{(b_1)} + 3h C_3^2 \underbrace{\mathbb{E}\int_0^h \Big| \tfrac12 (h^2-t^2)\nabla V(q^*)\Big|_{\ell^\infty}^2 {\rm d}t}_{(b_2)}\\ & \quad + \tfrac12 h^3 C_4^2 \underbrace{\mathbb{E}\int_0^h\int_0^t\int_0^s |q_\tau - q^*|_{\ell^\infty}^2{\rm d}\tau {\rm d}s {\rm d}t}_{(b_3)}\,. \end{align}\tag{83}\] Direct computation of \((b_1), (b_2)\) and the approach of 80 for \((b_3)\) give \[\label{eqn-bounds-b1-b2} (b_1) \leq \tfrac{4}{3}h^3\log(2d)\,,\quad (b_2) \leq \tfrac{2}{15}h^5 \mathbb{E}_\pi[|\nabla V|_{\ell^\infty}^2]\,,\quad (b_3) \leq \tfrac{2}{15}h^5\log(2d) + \tfrac{1}{420}h^7 \mathbb{E}_\pi[|\nabla V|_{\ell^\infty}^2]\,.\tag{84}\] Combining these via 81 and applying Lemma 5 gives ?? . ◻

Proof of Proposition 24. We use integral representations of the updates. Let \(\mathsf{U}_{\text{\normalfont\scshape ho},h}^\xi\) be the composition of a half-step of Hamiltonian dynamics \(\mathsf{H}= \mathsf{U}_{\text{\normalfont\scshape hmc}}^{h/2}\), followed by a full Ornstein–Uhlenbeck map \(\mathsf{O}\) with noise \(\xi\). We write \((\bar x, \bar p) \mathrel{\vcenter{:}}= \mathsf{U}_{\text{\normalfont\scshape ho},h}^\xi(x,p)\). With \((x(t), p(t)) = \mathsf{U}_{\text{\normalfont\scshape hmc}}^t(x,p)\) the continuous Hamiltonian trajectory 1 , we have \((x(t), p(t))\sim \pi \otimes \mathcal{N}(0,I)\) and \[\bar x = x + \tfrac h 2 p - \int_0^{h/2} \nabla V(x(t))(\tfrac h 2 - t){\rm d}t\,, \quad \bar p = \eta\Big(p - \int_0^{h/2}\nabla V(x(t)){\rm d}t\Big) + \sqrt{1-\eta^2}\xi\,.\]

Position component. Let \((\bar x(t), \bar p(t)) = \mathsf{U}_{\text{\normalfont\scshape hmc}}^t(\bar x, \bar p)\), so \((\bar x(t), \bar p(t)) \sim \pi \otimes \mathcal{N}(0, I)\). The position updates are \[\begin{align} q_{\mathsf{HOH}} \mathrel{\vcenter{:}}= & x + \tfrac{h}{2} p(1+\eta) - \int_0^{h/2}\nabla V(x(t))(\tfrac{h}{2}-t){\rm d}t - \int_0^{h/2}\nabla V(\bar x(t))(\tfrac{h}{2} - t) {\rm d}t\\ &+\tfrac{h}{2}\Big(\eta\Big(-\tfrac{h}{2}\nabla V(x) - \int_0^{h/2}\nabla^2 V(x(t))p(t)(\tfrac{h}{2}-t){\rm d}t \Big) + \sqrt{1-\eta^2}\xi \Big)\,, \end{align}\] \[q_{\mathsf{BAOAB}} \mathrel{\vcenter{:}}= y + \tfrac{h}{2} p' (1+\eta) - \tfrac{h^2}{4}(1+\eta) \nabla V(y) + \tfrac{h}{2}\sqrt{1-\eta^2}\xi\,,\] and their difference satisfies \[\begin{align} &|\Delta_q|_{2,\ell^\infty}\\ \leq &(1 + \tfrac{h^2}{4}\eta C_5) |x - y|_{2,\ell^\infty} + \tfrac{h}{2}(1+\eta)|p - p'|_{2,\ell^\infty}+ \Big| \int_0^{h/2}(\nabla V(x(t)) - \nabla V(y))(\tfrac{h}{2}-t){\rm d}t\\ & + \eta \tfrac{h}{2}\int_0^{h/2}\nabla^2 V(x(t)) p(t) (\tfrac{h}{2}-t){\rm d}t + \int_0^{h/2}(\nabla V(\bar x(t)) - \nabla V(y))(\tfrac{h}{2}-t){\rm d}t\Big|_{2,\ell^\infty}\,. \end{align}\] The last term is at most \[C_5\int_0^{h/2}(|x(t) - y|_{2,\ell^\infty}+|\bar x(t) - y|_{2,\ell^\infty})(\tfrac{h}{2}-t){\rm d}t + \tfrac{h\eta C_5}{2} \int_0^{h/2}|p(t)|_{2,\ell^\infty}(\tfrac{h}{2}-t){\rm d}t\,.\] Expanding the integral representation of \(x(t)\) gives \(|x(t) - y|_{2, \ell^\infty} \leq |x-y|_{2, \ell^\infty} + |tp - \int_0^t \nabla V(x(s))(t-s){\rm d}s|_{2, \ell^\infty}\), and by Lemmas 4 and 5, \[|tp - \int_0^t \nabla V(x(s))(t-s){\rm d}s|_{2, \ell^\infty} \leq h\sqrt{\log (2d)} + \tfrac{h^2}{4}\sqrt{\beta\log(2d)}\,.\] Similarly, with \(\bar x = x(h/2)\), \(|\bar x(t) - y|_{2, \ell^\infty} \leq |\bar x(t) - \bar x|_{2, \ell^\infty} + |\bar x - y|_{2, \ell^\infty}\), where \[\begin{align} &|\bar x(t) - \bar x|_{2, \ell^\infty} \leq h \sqrt{\log(2d)} + \tfrac{h^2}{4}\sqrt{\beta \log(2d)}\,,\\ & |\bar x - y|_{2, \ell^\infty} \leq |x - y|_{2, \ell^\infty} + h\sqrt{\log(2d)}+\tfrac{h^2}{4}\sqrt{\beta \log(2d)}\,. \end{align}\] Combining these and applying \(h \leq 1/(2\sqrt{\beta})\) gives ?? .

Momentum component. The momentum updates are \[\begin{align} p_{\mathsf{HOH}} &\mathrel{\vcenter{:}}= \eta\Big(p -\tfrac{h}{2}\nabla V(x) - \int_0^{h/2}\nabla^2 V(x(t))p(t)( \tfrac h 2 - t){\rm d}t \Big) + \sqrt{1-\eta^2}\xi\\ &\quad - \tfrac{h}{2}\nabla V(\bar x) - \int_0^{h/2}\nabla^2V(\bar x(t)) \bar p(t) ( \tfrac h 2 - t){\rm d}t\,, \end{align}\] \[p_{\mathsf{BAOAB}} \mathrel{\vcenter{:}}= \eta\Big(p' - \tfrac h 2 \nabla V(y) \Big) + \sqrt{1-\eta^2}\xi - \tfrac h 2 \nabla V(\hat{y})\,,\] where \(\hat{y} = y + \tfrac h2 (\eta+1)p' - \tfrac{h^2}4(\eta+1)\nabla V(y) + \tfrac h2 \sqrt{1-\eta^2}\xi\), so \[\begin{align} \Delta_p &= \eta(p - p') - \tfrac{h}{2}\eta(\nabla V(x) - \nabla V(y)) - \eta \int_0^{h/2}\nabla^2V(x(t))p(t)( \tfrac h 2 - t){\rm d}t\\ &\quad - \tfrac h 2 (\nabla V(\bar x) - \nabla V(\hat{y})) - \int_0^{h/2}\nabla^2V(\bar x(t)) \bar p(t) ( \tfrac h 2 - t) {\rm d}t\,. \end{align}\] Grouping the principal terms into \[\alpha_p \mathrel{\vcenter{:}}= \eta (p - p') - \tfrac h 2 \eta (\nabla V(x) - \nabla V(y)) - \tfrac h 2 (\nabla V(\hat{x}_c) - \nabla V(\hat{y}))\,,\] with \(\hat{x}_c = x + \tfrac h2(\eta+1)p - \tfrac{h^2}4(\eta+1)\nabla V(x) + \tfrac h2\sqrt{1-\eta^2}\xi\), \[\alpha_p = (\eta I - \tfrac{h^2(\eta+1)}{4}Q_2)(p-p')+(-\tfrac h 2 \eta Q_1 - \tfrac h 2 Q_2 + \tfrac{h^3(\eta+1)}{8}Q_2 Q_1)(x - y)\,,\] where \(Q_1 = \nabla^2 V([y, x])\), \(Q_2 = \nabla^2 V([\hat{x}_c, \hat{y}])\), and \(\nabla^2 V([v_1, v_2]) \mathrel{\vcenter{:}}= \int_0^1 \nabla^2 V(v_1 + s(v_2 - v_1)) {\rm d}s\). The remainder \(\varepsilon_p \mathrel{\vcenter{:}}= \Delta_p - \alpha_p\) is \(I_1 + I_2 + I_3 + I_4\), with \[\begin{align} I_1 &\mathrel{\vcenter{:}}= -\eta \int_0^{h/2}\nabla^2 V(x(t))p ( \tfrac h 2 -t){\rm d}t + \eta \tfrac{h^2}{4}\nabla^2 V([\hat{x}_c, \bar x])p- \eta \int_0^{h/2}\nabla^2 V(\bar x(t)) \tilde{p}( \tfrac h 2 - t){\rm d}t\\ & \quad + \tfrac{h^2}{4}\nabla^2 V([\hat{x}_c, x])\sqrt{1-\eta^2}\xi - \sqrt{1-\eta^2}\int_0^{h/2}\nabla^2 V(\bar x(t)) \xi ( \tfrac h 2 - t){\rm d}t\,, \end{align}\] where \(\tilde{p} \mathrel{\vcenter{:}}= p - \tfrac h2\nabla V(x) - \int_0^{h/2}\nabla^2 V(x(t)) p(t) (\tfrac h2 - t){\rm d}t\); \[I_2 \mathrel{\vcenter{:}}= \eta \int_0^{h/2} \nabla^2 V(x(t)) (p(t) - p) (\tfrac h 2 - t){\rm d} t\,,\] \[I_3 \mathrel{\vcenter{:}}= \tfrac h 2 \nabla^2 V([\hat{x}_c, \bar x])(-\tfrac{h^2}{4}(\eta+1)\nabla V(x) + \int_0^{h/2}\nabla V(x(t))( \tfrac h 2 - t){\rm d}t)\,,\] \[I_4 \mathrel{\vcenter{:}}= - \int_0^{h/2}\nabla^2 V(\bar x(t))(-t \nabla V(\bar x) - \int_0^t\nabla^2 V(\bar x(s))\bar p(s)(t-s){\rm d}s)( \tfrac h 2 - t){\rm d}t\,.\] These terms are bounded by \[\begin{align} &|I_1|_{2,\ell^\infty} \leq (\tfrac{21}{8}h^2 C_5 + \tfrac18h^4C_6) \sqrt{\log(2d)}\,,\quad |I_2|_{2,\ell^\infty} \leq ({\tfrac14}h^2 C_5 + \tfrac18h^4 C_6)\sqrt{\log(2d)}\,,\\ &|I_3|_{2, \ell^\infty} \leq {\tfrac{3}{8}}h^2 C_5\sqrt{\log(2d)}\,,\quad |I_4|_{2,\ell^\infty} \leq (\tfrac18h^2 C_5 + \tfrac18h^4C_6)\sqrt{\log(2d)}\,, \end{align}\] using \(h \leq 1/(2\sqrt{\beta})\), and combining them gives ?? . ◻

Proof of Proposition 25. By Proposition 24, unrolling the position error recurrence to \(\Delta_q^0 = 0\) gives \[\begin{align} |\Delta_q^l|_{2,\ell^\infty} &\leq |\Delta_q^{l-1}|_{2,\ell^\infty} + h (hC_5|\Delta_q^{l-1}|_{2,\ell^\infty} + |\Delta_p^{l-1}|_{2,\ell^\infty}) + \tfrac{27}{64}C_5h^3\sqrt{\log(2d)}\\ &\leq \sum_{i=1}^{l-1} h(hC_5|\Delta_q^i|_{2, \ell^\infty} + |\Delta_p^i|_{2,\ell^\infty}) + \tfrac{27}{64}C_5lh^3\sqrt{\log(2d)}\,. \end{align}\] Likewise, unrolling the momentum error recurrence and bounding the geometric sum of the constant term by \((1-\eta)^{-1}\) yields \[\begin{align} |\Delta_p^l|_{2,\ell^\infty} &\leq (\eta + \tfrac{h^2}{2}C_5)|\Delta_p^{l-1}|_{2,\ell^\infty} + (h C_5 + \tfrac{h^3}{4}C_6)|\Delta_q^{l-1}|_{2,\ell^\infty} + |\varepsilon_p|_{2,\ell^\infty}\\ &\leq \sum_{i=1}^{l-1}\eta^{i-1}\tfrac{h^2}{2}C_5|\Delta_p^{l-i}|_{2,\ell^\infty} + \sum_{i=1}^{l-1}\eta^{i-1}(hC_5 + \tfrac{h^3C_6}{4})|\Delta_q^{l-i}|_{2,\ell^\infty} + \tfrac{1}{1-\eta}|\varepsilon_p|_{2,\ell^\infty}\,. \end{align}\] Combining the two, with \(h \leq 1/(2\sqrt{\beta})\) and \(\eta \leq 1\), we obtain \[|\Delta_q^l|_{2,\ell^\infty} + \sqrt a |\Delta_p^l|_{2,\ell^\infty} \leq \sum_{i=1}^{l-1}h\Big(\tfrac32\tfrac{C_5}{\sqrt\beta}+ \tfrac{1}{16}\tfrac{C_6}{\beta^{3/2}} + \sqrt\beta\Big)(|\Delta_q^{l-i}|_{2,\ell^\infty}+\sqrt a |\Delta_p^{l-i}|_{2,\ell^\infty}) + E_l\,.\] Applying Grönwall’s inequality leads to ?? . ◻

Proof of Proposition 3. Following the coupling and one-step split 10 into the contraction term \((a)\) and discretization term \((b)\), we bound each.

Contraction term. By Proposition 8, \[(a) = |\tilde{T}_m^1(\{I - \tfrac{h^2}{2}H_i\}_{i=0}^{m-1})(X_k - Y_k)|_{2,\ell^\infty}\,,\] with \(H_i = \int_0^1 \nabla^2 V(\tau x_i + (1-\tau) y_i){\rm d}\tau\) and \(x_i, y_i\) the internal positions after \(i\) leap-frog steps, \(x_i \mathrel{\vcenter{:}}= \mathsf{\Pi}_1\mathsf{U}_{\text{\normalfont\scshape hmc},h}^i(X_k, \xi_k)\), \(y_i \mathrel{\vcenter{:}}= \mathsf{\Pi}_1\mathsf{U}_{\text{\normalfont\scshape hmc}, h}^i(Y_k, \xi_k)\). The contractive \(\ell^\infty\) bound of Proposition 11 gives \[(a) \leq \left(1 - \frac{\alpha}{400\beta}\right) |X_k - Y_k|_{2,\ell^\infty}\,.\]

Discretization error. The telescoping decomposition 12 and the component split 29 express \((b)\) through single integration-step errors. Representing each propagated factor as a matrix polynomial via Proposition 8, with \((q^*, p^*) \sim \pi \otimes \mathcal{N}(0,I)\), \[\begin{align} (b) & \leq\sum_{i=0}^{m-1} \underbrace{|\tilde{T}_{m-i-1}^1(\{I - \tfrac{h^2}{2}H_j\}_{j=0}^{m-i-2})(\mathsf{Q}_h^1 - \mathsf{Q}^h)(q^*, p^*)|_{2,\ell^\infty}}_{(b_1)}\\ & \quad +\sum_{i=0}^{m-1}\underbrace{|h \tilde{U}_{m-i-2}^1(\{I -\tfrac{h^2}{2}H_j\}_{j=1}^{m-i-2})(\mathsf{P}_h^1-\mathsf{P}^h)(q^*, p^*)|_{2,\ell^\infty}}_{(b_2)}\,, \end{align}\] where \(H_j = \int_0^1\nabla^2 V(\tau x_j + (1-\tau) y_j) {\rm d} \tau\) and \(x_j, y_j\) are the \(j\)-th internal leap-frog positions started at \(\mathsf{U}_{\text{\normalfont\scshape hmc},h}^1(q^*, p^*)\) and \(\mathsf{U}_{\text{\normalfont\scshape hmc}}^h(q^*,p^*)\), i.e.we redefine \(x_j \mathrel{\vcenter{:}}= \mathsf{\Pi}_1\mathsf{U}_{\text{\normalfont\scshape hmc},h}^j \mathsf{U}_{\text{\normalfont\scshape hmc},h}^1(q^*, p^*)\), \(y_j \mathrel{\vcenter{:}}= \mathsf{\Pi}_1\mathsf{U}_{\text{\normalfont\scshape hmc},h}^j\mathsf{U}_{\text{\normalfont\scshape hmc}}^h(q^*,p^*)\). By Proposition 11, almost surely \(|\tilde{T}_{m-i-1}^1|_{\ell^\infty} \leq 1 + \frac{\alpha}{400\beta} \leq \frac{401}{400}\) and \(|h\tilde{U}_{m-i-2}^1|_{\ell^\infty} \leq \tfrac{3}{2\sqrt\beta}\).

To bound \((b_1), (b_2)\) we apply Proposition 23, instantiating its general matrix \(M\) to the specific polynomials appearing here and determining the corresponding constants \(C_3, C_4\). For \((b_1)\), we take \(M = \tilde{T}_{m-i-1}^1\), which satisfies \(|M|_{\ell^\infty} \leq {\frac{401}{400}}\). Under the weak interaction assumption the integrated Hessian obeys \(|\int\nabla^2 V{\rm d}\nu|_{\ell^\infty} \leq \beta + \tfrac{\alpha}{50} \leq \tfrac{51}{50}\beta\), so \(|M\int\nabla^2 V{\rm d}\nu|_{\ell^\infty} \leq {\sqrt{2}\beta}\) and we may take \(C_3 = \sqrt{2}\beta\). For \((b_2)\), we take \(M = h\tilde{U}_{m-i-2}^1\) with \(|M|_{\ell^\infty} \leq \tfrac{3}{2\sqrt\beta}\), giving \(C_3 = \tfrac{3}{2\sqrt\beta}\cdot\tfrac{51}{50}\beta = \tfrac{153}{100}\sqrt\beta\) and, for a product of two Hessians, \(C_4 = \tfrac{3}{2\sqrt\beta}\cdot(\tfrac{51}{50}\beta)^2 = \tfrac{3}{2}(\tfrac{51}{50})^2\beta^{\frac{3}{2}}\). With \(h \leq 1/\sqrt{20\beta}\), Proposition 23 yields \[(b_1) \leq h^3\beta\sqrt{\log(2d)} \leq h^2\sqrt{\beta\log(2d)}\,, \quad (b_2) \leq 4h^2\sqrt{\beta\log(2d)}\,,\] and summing over the \(m\) steps, \[(b) \leq 5mh^2\sqrt{\beta \log(2d)} \leq 2h\sqrt{\log(2d)}\,,\] using \(mh = 1/\sqrt{20\beta}\) in the last step. Thus, \[|X_{k+1}-Y_{k+1}|_{2,\ell^\infty} \leq \left( 1 - \frac{\alpha}{400\beta} \right) |X_k - Y_k|_{2,\ell^\infty} + 2h\sqrt{\log(2d)}\,.\] Taking \(X_k, Y_k\) optimally coupled so that \(W_{2,\ell^\infty}(\rho_k, \pi) = |X_k - Y_k|_{2,\ell^\infty}\), we obtain \[W_{2,\ell^\infty}(\rho_{k+1},\pi) \leq \left( 1 - \frac{\alpha}{400\beta} \right) W_{2,\ell^\infty}(\rho_k, \pi) + 2h\sqrt{\log(2d)}\,.\] ◻

Proof of the HMC case of Theorem 1. Under \(C^{(O)} = 1/50\) and \(mh = 1/\sqrt{20\beta}\), taking \(k \to \infty\) in the inequality of Proposition 3 gives \[W_{2,\ell^\infty}(\pi_h, \pi) = O\left(\frac{\beta}{\alpha}h\sqrt{\log(2d)}\right)\,.\] ◻

Proof of Proposition 4. With \(X_k, Y_k\) as before, the multi-step split 11 bounds \(|X_{k+N} - Y_{k+N}|_{2,\ell^\infty}\) by its contraction term \((a)\) and discretization term \((b)\).

Contraction term. Proposition 8 represents the propagators as matrix polynomials, bounded in \(\ell^2\) by Proposition 9 and in \(\ell^\infty\) by Proposition 21, which extends Proposition 13. The \(\ell^2\) bound, converted to \(\ell^\infty\) at a cost of \(\sqrt d\), gives \[\label{eqn-part-a-proof-prop-sketch-sparse-bound-hmc} (a) \leq (1 - \tfrac{\alpha}{200\beta} )^N\sqrt{d}\, |X_k - Y_k|_{2,\ell^\infty}\,.\tag{86}\]

Discretization error. The telescoping decomposition 13 and the component split 29 , with the matrix-polynomial representation of Proposition 8, give \[\begin{align} (b) & \leq \sum_{l=0}^{N-1}\sum_{i=0}^{m-1} \underbrace{\Big|A_{(N-l-1)m+(m-i-1)}(\mathsf{Q}_h^1-\mathsf{Q}^h)(q^*, p^*)\Big|_{\ell^\infty}}_{(b_1)}\\ & \quad + \sum_{l=0}^{N-1}\sum_{i=0}^{m-1} \underbrace{\Big|B_{(N-l-1)m+(m-i-2)}(\mathsf{P}_h^1-\mathsf{P}^h)(q^*, p^*)\Big|_{\ell^\infty}}_{(b_2)}\,, \end{align}\] where, writing \(P_{l,i} \mathrel{\vcenter{:}}= \prod_{k \in\{m-i-1+jm\}_{j=0}^{N-l-2}}^{\longleftarrow}\tilde{T}_m^1(\{I - \tfrac{h^2}{2}H_{k+s} \}_{s=0}^{m-1})\) for the product of full-loop propagators over the \(N-l-1\) subsequent outer loops, \[A_{(N-l-1)m+(m-i-1)} = P_{l,i}\,\tilde{T}_{m-i-1}^1\,,\qquad B_{(N-l-1)m+(m-i-2)} = P_{l,i}\, h\tilde{U}_{m-i-2}^1\,.\] The Hessians \(H_j = \int_0^1 \nabla^2 V(\tau x_j + (1-\tau) y_j){\rm d} \tau\) are indexed by the total internal leap-frog step count: \[\begin{align} &0 \leq j \leq m-i-1: \quad x_{j} = \mathsf{\Pi}_1 \mathsf{U}_{\text{\normalfont\scshape hmc},h}^j \mathsf{U}_{\text{\normalfont\scshape hmc},h}^1(q^*, p^*)\,,\; y_j = \mathsf{\Pi}_1 \mathsf{U}_{\text{\normalfont\scshape hmc},h}^j \mathsf{U}_{\text{\normalfont\scshape hmc}}^h(q^*, p^*)\,;\\ & 0 \leq \hat{l} \leq N-l-1,\;0 \leq j < m-1:\\ &\quad x_{\hat{l} m+j+(m-i-1)} = \mathsf{\Phi}_{h,\xi_{k+l+\hat{l} + 1}}^j \mathsf{\Phi}_{h,\xi_{k+l+\hat{l}}}^m\cdots \mathsf{\Phi}_{h,\xi_{k+l+1}}^m\mathsf{\Pi}_1\mathsf{U}_{\text{\normalfont\scshape hmc},h}^{m-i-1}\mathsf{U}_{\text{\normalfont\scshape hmc},h}^1(q^*, p^*)\,, \\ &\quad y_{\hat{l} m+j+(m-i-1)} = \mathsf{\Phi}_{h,\xi_{k+l+\hat{l} + 1}}^j \mathsf{\Phi}_{h,\xi_{k+l+\hat{l}}}^m\cdots \mathsf{\Phi}_{h,\xi_{k+l+1}}^m\mathsf{\Pi}_1\mathsf{U}_{\text{\normalfont\scshape hmc},h}^{m-i-1}\mathsf{U}_{\text{\normalfont\scshape hmc}}^h(q^*, p^*)\,, \end{align}\] with the convention that for \(\hat{l} = 0\) only \(\mathsf{\Phi}_{h,\xi_{k+l+\hat{l} + 1}}^j\) is applied. By Propositions 21 and 23, for \(r = r_{(N-l-1)m+(m-i-1)}\), \[(b_1) \leq 3h^2 \sqrt{s_{r}\beta\log(2d)}\,,\quad (b_2) \leq 9h^2\sqrt{s_{r}\beta\log(2d)}\,,\] so summing over the \(Nm\) steps, \[\label{eqn-part-b-proof-prop-sketch-sparse-bound-hmc} (b) \leq 12h^2\sqrt{\beta\log(2d)}\sum_{i=0}^{Nm-1}\sqrt{s_{r_{i}}}\,.\tag{87}\] Combining 86 and 87 with \(X_k, Y_k\) optimally coupled completes the proof. ◻

Proof of the HMC case of Theorem 2. Take \(N = \lceil \frac{200\beta}{\alpha}\log(2\sqrt{d})\rceil\), so that the contraction factor \((1-\tfrac{\alpha}{200\beta})^N\sqrt{d} \leq \tfrac12\). Under \(s_k \leq C(k+1)^n\), an integral comparison gives \[\label{eqn-integral-sum-sparsity-hmc} \begin{align} &\sum_{i=0}^{Nm-1}\sqrt{s_{r_{i}}} \leq \sum_{i=0}^{Nm-1} \sqrt{C}\left(ih\sqrt\beta e + \tfrac{\log\sqrt d}{\log (5/3)} + 2\right)^{\frac{n}{2}}\\ &\quad \leq \int_0^{Nm}\sqrt{C}\left(yh\sqrt{\beta}e + \tfrac{\log\sqrt{d}}{\log(5/3)} + 2 \right)^{\frac{n}{2}}{\rm d}y \leq \sqrt{C}\frac{\left(Nmh\sqrt{\beta}e + \tfrac{\log\sqrt{d}}{\log(5/3)} + 2 \right)^{\frac{n}{2}+1}}{(n/2+1)h\sqrt{\beta}e}\,. \end{align}\tag{88}\] With this \(N\) and \(mh \leq 1/\sqrt{20\beta}\), Proposition 4 gives the recursion \[W_{2,\ell^\infty}(\rho_{k+N},\pi) \leq \tfrac12 W_{2,\ell^\infty}(\rho_k,\pi) + h \sqrt{\log(2d)}\left(O(\tfrac{\beta}{\alpha}\log(2d))\right)^{\frac{n}{2} + 1}\,,\] and taking \(k\to \infty\) gives \[W_{2,\ell^\infty}(\pi_h, \pi) = h\sqrt{\log(2d)}\left( O\left( \frac{\beta}{\alpha}\log(2d)\right)\right)^{\frac{n}{2} + 1}\,.\] ◻

Proof of the UL case of Theorem 1. We follow the coupling framework 18 with the weighted matrix representation of Section 4.1. Denote the initial distribution of \(X_0\) by \(\rho_0\). The BAOAB chain \((X_k, P_k)\) from \((X_0, P_0) \sim \rho_0 \otimes \mathcal{N}(0,I)\) and the HOH chain \((Y_k, P_k')\) from \((Y_0, P_0') \sim \pi \otimes \mathcal{N}(0,I)\) are coupled through the same noise \(\{\xi_i\}\). Introducing the auxiliary process \((\bar Y_k, \bar P'_k) = \mathsf{U}_{\text{\normalfont\scshape baoab},h}^{\xi_{k-1}}\cdots \mathsf{U}_{\text{\normalfont\scshape baoab},h}^{\xi_0}(Y_0, P'_0)\), which starts at \((Y_0, P_0')\) but runs BAOAB with the same noise, the split 18 becomes \[\label{eqn-prop-baoab-weak-bias-iteration-decomposition} |X_k - Y_k|_{2,\ell^\infty} \leq \underbrace{|X_k - \bar Y_k|_{2,\ell^\infty}}_{(a)} + \underbrace{|\bar Y_k - Y_k|_{2,\ell^\infty}}_{(b)}\,,\tag{89}\] the contraction between two BAOAB chains \((a)\) and the BAOAB-vs-HOH discretization error \((b)\).

Contraction term. Using the decomposition \(\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{\xi_{k-1}}\cdots \mathsf{U}_{\text{\normalfont\scshape baoab},h}^{\xi_0} = \mathsf{U}_{\text{\normalfont\scshape ab},h} \mathsf{U}_{\text{\normalfont\scshape abao},h}^{\xi_{k-1}} \cdots \mathsf{U}_{\text{\normalfont\scshape abao},h}^{\xi_{1}} \mathsf{U}_{\text{\normalfont\scshape bao},h}^{\xi_0}\) of Section 2.3 and the weighted position propagators 26 , \[\label{eqn-weighted-matrix-representation-weak-interaction} X_k - \bar Y_k = \begin{bmatrix} I & 0 \end{bmatrix} W^{-1} M_{\text{\normalfont\scshape ab}}^w (H_{k}) M_{\text{\normalfont\scshape abao}}^w(H_{k-1}) \cdots M_{\text{\normalfont\scshape abao}}^w(H_1)M_{\text{\normalfont\scshape bao}}^w(H_0) W\Delta Z_0\,,\tag{90}\] with \(\Delta Z_0 = (X_0-Y_0, P_0 - P'_0)\) and the matrices defined through the intermediate iterates \[\begin{align} &z_1\mathrel{\vcenter{:}}= \mathsf{U}_{\text{\normalfont\scshape bao},h}^{\xi_0}(X_0, P_0)\,, \quad z_1'\mathrel{\vcenter{:}}= \mathsf{U}_{\text{\normalfont\scshape bao},h}^{\xi_0}(Y_0, P_0')\,, \quad z_{k+1} \mathrel{\vcenter{:}}= \mathsf{U}_{\text{\normalfont\scshape ab},h}z_{k}\,, \quad z_{k+1}' \mathrel{\vcenter{:}}= \mathsf{U}_{\text{\normalfont\scshape ab},h} z_{k}'\,,\\ &z_i \mathrel{\vcenter{:}}= \mathsf{U}_{\text{\normalfont\scshape abao},h}^{\xi_{i-1}} \cdots \mathsf{U}_{\text{\normalfont\scshape abao},h}^{\xi_1}z_1\,, \quad z_i' \mathrel{\vcenter{:}}= \mathsf{U}_{\text{\normalfont\scshape abao},h}^{\xi_{i-1}} \cdots \mathsf{U}_{\text{\normalfont\scshape abao},h}^{\xi_1}z_1'\,, \quad 2\leq i \leq k\,, \end{align}\] satisfying \[\begin{align} & W(z_1 - z_1') = M_{\text{\normalfont\scshape bao}}^w(H_0)W\Delta Z_0\,, \quad W(z_{k+1} - z_{k+1}') = M_{\text{\normalfont\scshape ab}}^w(H_k)W(z_{k}-z_{k}')\,,\\ & W(z_i -z_i') = M_{\text{\normalfont\scshape abao}}^{w}(H_{i-1}) W(z_{i-1}- z_{i-1}')\,, \quad 2\leq i \leq k\,. \end{align}\] By Proposition 12 and \(|\begin{bmatrix} I & 0 \end{bmatrix}W^{-1}|_{\ell^\infty} = 1 + \frac{b}{\sqrt{a-b^2}}\), \[(a) \leq 6\Big(1 + \tfrac{b}{\sqrt{a-b^2}}\Big) ( 1 - \tfrac12c(h))^{k-1}|\Delta Z_0|_{2,\ell^\infty_w}\,.\]

Discretization error. Write \((\Delta_q^{l}, \Delta_p^{l}) = (\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{l} - \mathsf{U}_{\text{\normalfont\scshape hoh},h}^{l})(q^*, p^*)\). The matrix representation of \(\mathsf{\Pi}_1\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{k-(i+1)\tilde{l}}(\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{\tilde{l}} - \mathsf{U}_{\text{\normalfont\scshape hoh},h}^{\tilde{l}})(q^*, p^*)\) parallels 90 , with \(\Delta Z_0 = (\Delta_q^{\tilde{l}}, \Delta_p^{\tilde{l}})\) and Hessians along the intermediate iterates. Each term of the telescoping decomposition 19 is thus bounded by \[\label{eqn-baoab-discretization-error-individual-term-with-propagation} \begin{align} &|\mathsf{\Pi}_1\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{k-(i+1)\tilde{l}}(\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{\tilde{l}} - \mathsf{U}_{\text{\normalfont\scshape hoh},h}^{\tilde{l}}) (q^*, p^*)|_{2,\ell^\infty}\\ &\qquad \leq 6\Big(1 + \tfrac{b}{\sqrt{a-b^2}}\Big) ( 1 - \tfrac12c(h))^{k-(i+1)\tilde{l} - 1}|(\Delta_{q}^{\tilde{l}}, \Delta_p^{\tilde{l}})|_{2,\ell^\infty_w}\,, \end{align}\tag{91}\] \[\label{eqn-baoab-last-step-discretization-error} \left|\mathsf{\Pi}_1\left(\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{k-\lfloor k/\tilde{l} \rfloor \tilde{l}} - \mathsf{U}_{\text{\normalfont\scshape hoh},h}^{k - \lfloor k/\tilde{l} \rfloor \tilde{l}}\right)(q^*, p^*)\right|_{2,\ell^\infty} \leq \Big(1 + \tfrac{b}{\sqrt{a-b^2}}\Big)|(\Delta_q^{k - \lfloor k/\tilde{l}\rfloor \tilde{l}}, \Delta_p^{k - \lfloor k/\tilde{l}\rfloor \tilde{l}})|_{2,\ell^\infty_w}\,.\tag{92}\] In 91 , when \(i = \lfloor k/\tilde{l} \rfloor -1\) there may be no internal ABAO steps and the exponent can be negative, so we instead bound the operator by \(6(1+\tfrac{b}{\sqrt{a-b^2}})\), accounting for the \(\mathsf{B}\mathsf{A}\mathsf{O}\) and \(\mathsf{A}\mathsf{B}\) prefactors. Proposition 25 gives for any \(l \leq \tilde{l}\) \[|(\Delta_q^l, \Delta_p^l)|_{2,\ell^\infty_w} \leq \exp\left[(\tilde{l}-1)h\Big(\tfrac32\tfrac{C_5}{\sqrt\beta} + \tfrac{1}{16}\tfrac{C_6}{\beta^{3/2}} + \sqrt\beta\Big)\right] E_{\tilde{l}}\,.\] Substituting these into 89 , with \(X_0, Y_0\) optimally coupled so \(W_{2,\ell^\infty}(\rho_0, \pi) = |X_0-Y_0|_{2,\ell^\infty}\) and \(P_0' = P_0\), \[\label{eqn-term-by-term-k-step-error-UL-weak} \begin{align} &|X_k - Y_k|_{2,\ell^\infty} \leq 6 \Big(1 + \tfrac{b}{\sqrt{a-b^2}}\Big) \left( 1 - \tfrac12c(h)\right)^{k-1}W_{2,\ell^\infty}(\rho_0,\pi) + \Big(1 + \tfrac{b}{\sqrt{a-b^2}}\Big)\\ &\quad \cdot\left( 7 + \sum_{i=0}^{\lfloor k / \tilde{l}\rfloor-2} 6(1-\tfrac12c(h))^{k-(i+1)\tilde{l}-1}\right) \exp\left[(\tilde{l}-1)h\Big(\tfrac32\tfrac{C_5}{\sqrt\beta} + \tfrac{1}{16}\tfrac{C_6}{\beta^{3/2}}+\sqrt\beta\Big)\right] E_{\tilde{l}}\,. \end{align}\tag{93}\] Under the weak interaction assumption \(|H_i^{(O)}|_{\ell^\infty} \leq \alpha/20\), the constants of Proposition 24 satisfy \(C_5 \leq \tfrac{21}{20}\beta\) and \(C_6 \leq (\tfrac{21}{20})^2\beta^2\), so that \(C_5/\beta\) and \(C_6/\beta^2\) are \(O(1)\). Since \((\tilde{l}-1)h\sqrt\beta = O(1)\) by the choice \(\tilde{l} = \lceil 1/(2\sqrt\beta h)\rceil\), the exponential is \(O(1)\), and substituting \(E_{\tilde{l}}\) gives \[\exp\left[(\tilde{l}-1)h\Big(\tfrac32\tfrac{C_5}{\sqrt\beta} + \tfrac{1}{16}\tfrac{C_6}{\beta^{3/2}} + \sqrt\beta\Big)\right] E_{\tilde{l}} = O\Big(\frac{h^2\sqrt\beta}{1-\eta}\sqrt{\log(2d)}\Big)\,.\]

As \(k \to \infty\), the first term vanishes. For the geometric sum, with \(c(h) = \frac{\alpha h^2}{4(1-\eta)}\) and \(\tilde{l} \geq 1/(2\sqrt\beta h)\), the inequality \(1-(1-x)^{\tilde{l}} \geq \tilde{l} x - \tfrac12(\tilde{l} x)^2\) for \(x>0\) gives \[1 - (1- \tfrac12c(h))^{\tilde{l}} \geq \frac{63}{1024}\frac{\alpha h}{\sqrt\beta (1-\eta)}\,,\] so that \(\sum_{i=0}^{\lfloor k/\tilde{l}\rfloor - 2} 6(1-\tfrac12 c(h))^{k-(i+1)\tilde{l} - 1} \leq 6/(1-(1-\tfrac12 c(h))^{\tilde{l}}) = O(\tfrac{\sqrt\beta(1-\eta)}{\alpha h})\). Since \(1 + \tfrac{b}{\sqrt{a-b^2}} = O(1)\), 93 becomes \[W_{2,\ell^\infty}(\pi_h, \pi) \leq O(1)\cdot O\Big(\frac{\sqrt\beta(1-\eta)}{\alpha h}\Big)\cdot O\Big(\frac{h^2\sqrt\beta}{1-\eta}\sqrt{\log(2d)}\Big) = O\Big(\frac{\beta}{\alpha}h\sqrt{\log(2d)}\Big)\,.\] ◻

Proof of the UL case of Theorem 2. As in the HMC sparse case, we induce an \(\ell^\infty\) contraction from the \(\ell^2\) contraction over multiple steps. We partition the \(k\) steps into blocks of \(B\) steps, with \(B\) specified below. We use the three processes \((X_k, P_k), (Y_k, P'_k), (\bar Y_k, \bar P'_k)\) of the proof of the UL case of Theorem 1 with \(P_0 = P_0' = \bar P'_0\), and the decomposition \[\label{eqn-prop-baoab-sparse-bias-iteration-decomposition} |X_k - Y_k|_{2,\ell^\infty} \leq \underbrace{|X_k - \bar Y_k|_{2,\ell^\infty}}_{(a)} + \underbrace{|\bar Y_k - Y_k|_{2,\ell^\infty}}_{(b)}\,.\tag{94}\]

Contraction term. Corollary 1 gives \[\label{eqn-part-a-baoab-sparse} (a) \leq 14\sqrt{d}(1-c(h))^{\frac{k-1}{2}}|X_0 - Y_0|_{2,\ell^\infty}\,,\tag{95}\] which vanishes as \(k \to \infty\).

Discretization error. The telescoping decomposition 19 accumulates a local error every \(\tilde{l} = \lceil \frac{1}{2\sqrt\beta h} \rceil\) steps, each propagated through the later steps. Write \((\Delta_q^l, \Delta_p^l) = (\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{l} - \mathsf{U}_{\text{\normalfont\scshape hoh},h}^{l})(q^*, p^*)\) for an \(l\)-step error, with \((q^*, p^*) \sim \pi \otimes \mathcal{N}(0,I)\), and recall that \(\mathsf{U}_{\text{\normalfont\scshape baoab},h}^l, \mathsf{U}_{\text{\normalfont\scshape hoh},h}^l\) are the \(l\)-fold compositions omitting the coupled noise. For the sparse case the constants of Proposition 24 are \(C_5 = \beta\sqrt{s_1}\) and \(C_6 = \beta^2\sqrt{s_2}\), so Proposition 25 bounds the \(\tilde{l}\)-step error, \[|(\Delta_q^{\tilde{l}}, \Delta_p^{\tilde{l}})|_{2,\ell^\infty_w} \leq \exp\left[(\tilde{l}-1)h\Big(\tfrac32\tfrac{C_5}{\sqrt\beta} + \tfrac{1}{16}\tfrac{C_6}{\beta^{3/2}} + \sqrt\beta\Big)\right] E_{\tilde{l}} = O\Big(\frac{h^2\sqrt\beta}{1-\eta}\sqrt{\log(2d)}\Big)\,.\] Setting the block length \(B \mathrel{\vcenter{:}}= K_1 \tilde{l}\) with \(K_1\) specified later and grouping the trailing steps in 19 into full blocks of \(B\) steps followed by residual steps, \[\label{eqn-multi-step-baoab-discretization-error-block-decomposition-propagation} \begin{align} (b) & \leq \sum_{i=0}^{\lceil\frac{\lfloor k/ \tilde{l} \rfloor}{K_1}\rceil}\sum_{k_1=0}^{K_1-1}\left| \mathsf{\Pi}_1 \mathsf{U}_{\text{\normalfont\scshape baoab},h}^{j_{k_1}}\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{iB}(\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{\tilde{l}} - \mathsf{U}_{\text{\normalfont\scshape hoh},h}^{\tilde{l}})(q^*, p^*)\right|_{2,\ell^\infty}\\ & \quad + \left| \mathsf{\Pi}_1\left(\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{k - \lfloor k/\tilde{l} \rfloor \tilde{l}} - \mathsf{U}_{\text{\normalfont\scshape hoh},h}^{k - \lfloor k / \tilde{l} \rfloor \tilde{l}} \right)(q^*, p^*)\right|_{2,\ell^\infty}\,, \end{align}\tag{96}\] where \(i\) counts the full \(B\)-step blocks, \(k_1 \in \{0, \dots, K_1-1\}\) indexes the residual \(\tilde{l}\)-step groups, and \(j_{k_1} \mathrel{\vcenter{:}}= k - \lfloor k /\tilde{l}\rfloor \tilde{l} + k_1\tilde{l}\) is the corresponding number of trailing steps.

Each summand in the first sum is one \(\tilde{l}\)-step error propagated through \(iB\) full-block steps and \(j_{k_1}\) trailing steps. For two initial states \((x,p), (y, p')\), Proposition 8 and 27 give the weighted matrix-polynomial form of the propagated position difference over the residual steps, \[\label{eqn-baoab-individual-discretization-error-matrix-polynomial} \mathsf{\Pi}_1\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{j_{k_1}} (x,p) - \mathsf{\Pi}_1\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{j_{k_1}} (y,p') = \begin{bmatrix} \tilde{T}_{j_{k_1}}^\eta & - \tfrac{b}{\sqrt{a-b^2}} \tilde{T}_{j_{k_1}}^\eta + \tfrac{1}{\sqrt{a-b^2}}h\tilde{U}_{j_{k_1}-1}^\eta \end{bmatrix} W \begin{bmatrix} x- y \\ p - p' \end{bmatrix}\,,\tag{97}\] the polynomials evaluated at matrices of the form \(I - \tfrac{h^2}{2}H\) with \(H\) the integrated Hessians along the intermediate iterates of the two trajectories. Applying this with \((x,p) = \mathsf{U}_{\text{\normalfont\scshape baoab},h}^{iB}\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{\tilde{l}} (q^*, p^*)\) and \((y, p') = \mathsf{U}_{\text{\normalfont\scshape baoab},h}^{iB}\mathsf{U}_{\text{\normalfont\scshape hoh},h}^{\tilde{l}}(q^*, p^*)\), the \(iB\) full-block steps contract the \(\tilde{l}\)-step error by Corollary 1, \[\label{eqn-baoab-contraction-in-discretization-error} |(x-y, p-p')|_{2,\ell^\infty_w} \leq 7\sqrt{2d}(1-c(h))^{\frac{iB-1}{2}} |(\Delta_q^{\tilde{l}}, \Delta_p^{\tilde{l}})|_{2,\ell^\infty_w}\,,\tag{98}\] while Proposition 14 bounds the \(\ell^\infty\) norm of the matrix-polynomial row in 97 . Together they give \[\begin{align} &\left| \mathsf{\Pi}_1 \mathsf{U}_{\text{\normalfont\scshape baoab},h}^{j_{k_1}}\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{iB}(\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{\tilde{l}} - \mathsf{U}_{\text{\normalfont\scshape hoh},h}^{\tilde{l}})(q^*, p^*)\right|_{2,\ell^\infty}\\ &\qquad \leq 7\sqrt{2d}\Big(14\sqrt{2}+1+\tfrac{2}{\sqrt{3}}\Big)\sqrt{s_{r_{j_{k_1}}}}\, (1-c(h))^{\frac{iB-1}{2}}|(\Delta_q^{\tilde{l}}, \Delta_p^{\tilde{l}})|_{2,\ell^\infty_w}\,. \end{align}\] For \(i=0\) where there are no intermediate full blocks we instead employ \[\left| \mathsf{\Pi}_1 \mathsf{U}_{\text{\normalfont\scshape baoab},h}^{j_{k_1}}(\mathsf{U}_{\text{\normalfont\scshape baoab},h}^{\tilde{l}} - \mathsf{U}_{\text{\normalfont\scshape hoh},h}^{\tilde{l}})(q^*, p^*)\right|_{2,\ell^\infty} \leq \left(14\sqrt{2}+1+\tfrac{2}{\sqrt{3}}\right)\sqrt{s_{r_{j_{k_1}}}}|(\Delta_q^{\tilde{l}}, \Delta_p^{\tilde{l}})|_{2,\ell^\infty_w}\,.\] We now fix \(K_1 = \lceil \frac{16\sqrt\beta(1-\eta)\log(14\sqrt{2d})}{h\alpha}+1\rceil\), so that the per-block factor contracts, \(7\sqrt{2d}(1-c(h))^{\frac{B-1}{2}} \leq \tfrac12\), making the contraction over \(i\) blocks at most \(2^{-i}\). Substituting the \(\tilde{l}\)-step error bound and these per-term estimates into 96 , together with the last-term bound 92 , and summing the geometric factor \(2^{-i}\) over \(i\), we obtain \[\label{eqn-part-b-baoab-sparse} (b) \leq \left(1 + 2\left( 14\sqrt{2}+1+\tfrac{2}{\sqrt{3}}\right)\sum_{k_1=0}^{K_1-1} \sqrt{s_{r_{j_{k_1}}}} \right)O\Big(\frac{h^2\sqrt\beta}{1-\eta}\sqrt{\log(2d)}\Big)\,.\tag{99}\] By an integral comparison as in 88 , with \(s_{r_i}\) increasing in \(i\), \[\begin{align} \sum_{k_1=0}^{K_1-1} \sqrt{s_{r_{j_{k_1}}}} \leq \sum_{k_1=1}^{K_1}\sqrt{s_{r_{k_1\tilde{l}}}} &\leq \sqrt{C}\frac{\left(\tfrac{\tilde{l}e^2\beta (K_1+1)}{1-\eta}h^2+\log\sqrt{d}+2\right)^{\frac{n}{2} + 1}}{(\frac{n}{2} + 1) \frac{\tilde{l} e^2\beta}{1-\eta}h^2}\\ &= \frac{1-\eta}{\sqrt\beta h}\left(O(\tfrac{\beta}{\alpha}\log(2d))\right)^{\frac{n}{2} +1}\,. \end{align}\] Combining with 95 and 99 in 94 and taking \(k \to \infty\), \[W_{2,\ell^\infty}(\pi_h, \pi) = h \sqrt{\log(2d)} \left(O(\tfrac{\beta}{\alpha}\log(2d))\right)^{\frac{n}{2} +1}\,.\] ◻