B.1 Precise assumptions

Let

(10) \begin{equation} N_k = \{ x \in \partial U_k | \exists ! i \text{ such that } \exists u \in W_i \text{ such that } \phi(u) = x\} \end{equation}

be the set of points of $\partial U_k$ for which the normal $n(x) = n(u)$ is well-defined. Since $\dim_H((\overline{U_k})^{\mathrm{o}}\setminus U_k) \leq d_k - 2$ , for all points for which the normal exists, the boundary separates $U_k$ and $\mathbb{R}^{d_k}\setminus U_k$ , and does not correspond to an ‘internal’ boundary of $U_k$ that is removed. By convention, we assume that n(x) is the outer normal.

B.2. Proof of Theorem 2

We abuse notation and write $\partial_v f$ for the derivative at $t=0$ of $f(x+tv,v)$ with respect to t for a fixed v, which corresponds to the term $\Xi f$ .

B.2.2. Integration of the infinitesimal generator over E.

Proposition 9. Let $f\in \mathcal{D}(A)$ and let $k \in K$ . For all $v\in \mathcal{V}^k$ we have

\begin{equation*} \begin{split} \int_{U_k} \pi_k(x) \partial_v f(k,x,v)\;\textrm{d} x=-\int_{U_k} f(k,x,v) \partial_v \pi_k(x) \;\textrm{d} x\;\textrm{d} v \\[5pt] + \int_{\partial U_k\cap N_k}\pi_k(x) f(k,x,v)\,|\langle n(x), v\rangle|\;\textrm{d}\sigma(x), \end{split} \end{equation*}

where $\sigma$ is the Lebesgue measure of the boundary (seen as a Riemannian manifold).

Proof. We would like to use an integration-by-parts result on the integral in question, which is precisely detailed in Appendix C.

So now let $v \in \mathcal{V}^k$ . Assumptions 3 and 4 imply that $U_k$ satisfies Assumption 6 of Appendix C. Furthermore, ; thus for all $x\in U_k$ , $t\rightarrow f(k,x+tv,v)$ is absolutely continuous. Finally, Lemma 5 implies that $\partial_v f(k,x,v) \in L_1(\mu)$ ; hence we can use Proposition 10 of Appendix C on $U_k$ with the function $z\mapsto f(z) \pi(z)$ . In this context, for $x \in \partial U_k$ , if $\langle n(x),v \rangle > 0$ ,

\begin{align*} \pi_k(x^-) f(k,x^-,v) &= \lim_{t\uparrow 0}\pi_k(x+tv) f(k,x+tv,v) = \pi_k(x) f(k,x,v), \\[5pt] \pi_k\big(x^+\big) f\big(k,x^+,v\big) &= 0; \end{align*}

otherwise, if $\langle n(x),v \rangle > 0$ ,

\begin{align*} \pi_k(x^-) f(k,x^-,v) &= 0, \\[5pt] \pi_k(x^+) f(k,x^+,v) &= \lim_{t\downarrow 0}\pi_k(x+tv) f(k,x+tv,v) = \pi_k(x) f(k,x,v). \end{align*}

This yields

\begin{align*} \int_{U_k} \pi_k(x) \partial_v f(k,x,v)\;\textrm{d} x &=-\int_{U_k} f(k,x,v) \partial_v \pi_k(x) \;\textrm{d} x\;\textrm{d} v\\[5pt] &+ \int_{\partial U_k\setminus C_v}\pi_k(x) f(k,x,v)\, \langle n(x), v \rangle \;\textrm{d}\sigma(x), \end{align*}

where we have removed the absolute value around $\langle n(x), v \rangle$ to account for the sign difference of $\pi_k(x^-) f(k,x^-,v) - \pi_k(x^+) f(k,x^+,v)$ . Furthermore, $\partial U_k\setminus C_v \subset \partial U_k \cap N_k$ , and these two sets differ by a set of zero measure. Hence

\begin{equation*} \begin{split} \int_{U_k}\partial_v \pi_k(x) f(k,x,v)\;\textrm{d} x=-\int_{U_k} f(k,x,v) \partial_v \pi_k(x) \;\textrm{d} x\;\textrm{d} v\\[5pt] + \int_{\partial U_k\cap N_k}\pi_k(x) f(k,x,v)\, \langle n(x), v\rangle\;\textrm{d}\sigma(x), \end{split} \end{equation*}

which concludes the proof.

Lemma 6. We have

\begin{align*} \int_{E_k} \partial_v f(k,x,v) &\mu(k,x,v) \;\textrm{d} x \;\textrm{d} v = -\int_{U_k\times \mathcal V^k} f(k,x,v) \nabla_x \mu \cdot v \;\textrm{d} x \;\textrm{d} v\\[5pt] &+ \int_{(\partial U_k \cap N_k)\times \mathcal V^k} f(k,x,v) \mu(k,x,v) \langle v, n(x) \rangle \;\textrm{d}\sigma(x) \;\textrm{d} v.\end{align*}

Proof. Since $\partial_v f$ is in $L_1(\mu)$ , $\partial_v f(k,x,v) \mu(k,x,v)$ is integrable. Hence by Fubini’s theorem,

\begin{align*} \int_{U_k\times\mathcal{V}^k} \partial_v f(k,x,v) \mu(k,x,v) \;\textrm{d} x\;\textrm{d} v = \int_{\mathcal{V}^k}\int_{U_k} \partial_v f(k,x,v) \mu(k,x,v) \;\textrm{d} x\;\textrm{d} v.\end{align*}

Using Proposition 9,

\begin{align*} \int_{U_k\times\mathcal{V}^k} \partial_v f(k,x,v) \mu(k,x,v) &\;\textrm{d} x\;\textrm{d} v = - \int_{\mathcal{V}^k}\int_{U_k} f(k,x,v) (\nabla \pi_k(x) \cdot v) p_k(v) \;\textrm{d} x\;\textrm{d} v \\[5pt] &+ \int_{\mathcal V^k}\int_{\partial U_k \cap N_k} f(k,x,v) \mu(k,x,v) \langle v, n(x) \rangle \;\textrm{d}\sigma(x) \;\textrm{d} v.\end{align*}

Using (iii) from Assumption 2, $f(k,x,v) (\nabla \pi_k(x) \cdot v) p_k(v)$ is integrable, and $f(k,x,v) \mu(k,x,v)$ is bounded. Hence we can use Fubini a second time on both terms to get the result.

C.1. Integration over open domain

We extend f to $\mathbb R^{n}\setminus \partial U$ with $f(y)=0$ for every $y\notin \overline U$ , so that we can suppose that $U=\mathbb R^n\setminus \partial U$ .

Theorem 4. Let U be an open set of $\mathbb R^n$ satisfying Assumption 6, for some fixed $v \in \mathbb R^n$ with $\| v\|=1$ . Let

\begin{equation*} \mathcal N_v=\left(\bigcup_{i=1}^k(\phi_i(\partial W_i)\cup M_i^0)\right)\bigcup\left(\bigcup_{1\leq i<j\leq k} M_i\cap M_j)\right) \end{equation*}

be the set on which normals are ill-defined. Then for any f satisfying Assumption 7, the following hold:

1. 1. We have $\dim_H \mathcal N_v\leq 2$ and $\dim_H p_v(\mathcal N_v)\leq n-2$ .

2. 2. For each $y\in\partial U\setminus \mathcal N_v$ , the limits

   \begin{align*} \lim_{t\downarrow 0}f(y+tv)=f(y^+) \text{ and } \lim_{t\uparrow 0}f(y+tv)=f(y^-) \end{align*}
   exist.

3. 3. The normal n(y) is well-defined at each point $y\in\partial U\setminus \mathcal N_v$ , and

   \begin{align*} \int_U\partial_v f(y)\;\textrm{d} y=\int_{\partial U\setminus \mathcal N_v}(f(y^-)-f(y^+))\,|n(y)\cdot v|\;\textrm{d}\sigma(y), \end{align*}
   where $\sigma$ is the Lebesgue measure on $\partial U$ .

We can use the theorem with a product $f=g\pi$ where

• • $g\;:\;U\rightarrow\mathbb R$ is measurable, bounded, and absolutely continuous on each set $U\cap(y+\mathbb R v)$ , $y\in\mathbb R^n$ , and $\partial_t g\in L^1(U)$ ; and

• • $\pi\;:\;U\rightarrow\mathbb R$ is in $C^1(U)\cap L^1(U)$ and is bounded with bounded derivatives, $\lim_{\|y\|\rightarrow\infty}\pi(y)=0$ , and the derivative $\partial_v\pi\in L^1(U)$ .

C.2. Proof of Theorem 4

For the first point, let

\begin{align*} \mathcal N_v&=\left(\bigcup_{i=1}^k(\phi_i(\partial W_i)\cup M_i^0)\right)\bigcup\left(\bigcup_{1\leq i<j\leq k} M_i\cap M_j)\right),\\[5pt] \mathcal N&=p_v^{-1}(p_v(\mathcal N_v)) \text{ and } V=H\setminus p_v(\mathcal N_v).\end{align*}

By Assumptions 6.1, 6.4, and 6.5, $\dim_H\!(p_v(\mathcal N_v))\leq n-2$ ; therefore $p_v(\mathcal N_v)$ has zero H-Lebesgue measure.

For the second point, the fact that $f(y+\!)$ and $f(y-\!)$ exist is a direct consequence of Assumption 7(i).

Finally we consider the third point. For all $z\in H$ , define

\begin{align*} E(z)=\{t\in \mathbb R\;:\;z+tv\in\partial U\}.\end{align*}

By Assumption 6.6, the set E(z) has 2k elements at most for all $z\in V$ . Set

\begin{align*}U'=U\setminus\mathcal N.\end{align*}

Since $\dim_H(p_v(\mathcal N_v))\leq n-2$ , it follows that $\dim_H \mathcal N\leq n-1$ . Therefore,

\begin{align*}\int_U\partial_vf(y)\;\textrm{d} y=\int_{U'}\partial_vf(y)\;\textrm{d} y.\end{align*}

By Fubini’s theorem,

\begin{align*}\int_{U'}\partial_v f(y)\;\textrm{d} y&=\int_V\left(\int_{\mathbb R\setminus E(z)}\partial_tf(z+tv)\;\textrm{d} t\right)\;\textrm{d} z.\end{align*}

By Assumption 7, for almost all $z\in V$ and for each connected component (a, b) of $\mathbb R\setminus E(z)$ ,

\begin{align*}\int_{(a,b)}\partial_v f(z+tv)\;\textrm{d} t=f((z+tb)^-)-f((z+ta)^+).\end{align*}

Taking into account that $\lim_{t\rightarrow\pm\infty}f(z+tv)=0$ , we obtain

\begin{align*} \int_{U}\partial_tf(y)\;\textrm{d} y=\int_V\sum_{t\in E(z)}(f((z+tv)^-)-f((z+tv)^+))\;\textrm{d} z.\end{align*}

Now we want to see that the latter integral is equal to

\begin{align*}\int_{\partial U\setminus \mathcal N_v}(f(y^-)-f(y^+))\,|n(y)\cdot v|\;\textrm{d}\sigma(y).\end{align*}

Set

\begin{align*} I&=\{1,\dots,k\}\times\{+,-\},\\[5pt] J(z)&=\{(i,s)\in I\;:\;\exists t\in E(z), z+tv\in M_i^s\} \text{ for } z\in H,\\[5pt] V_J&=\{z\in V \;:\; J(z)=J\} \text{ for } J\subset I, \text{ and }\\[5pt] V_i^s&=p_v(M_i^s) \text{ for } (i,s)\in I.\end{align*}

By Assumption 6.6, for each $(i,s)\in I$ , the map $p_v\circ \phi_i\;:\;W_i^s\rightarrow V_i^s$ is a bijection, so that we can define the map $F_i^s\;:\;V_i^s\rightarrow M_i^s$ by $F_i^s(z)=\phi_i((p_v\circ \phi_i)^{-1}(z))$ . Using the definition of the set $V=H\setminus\mathcal N$ , we see that for each $z\in V$ and each $t\in E(z)$ , there exists $(i,s)\in I$ unique such that $z+tv=F_i^s(z)\in M_i^s$ . Furthermore, for each $(i,s)\in I$ , $V_i^s\cap V=\cup_{J\ni (i,s)}V_J$ . It follows that

\begin{align*} \int_U\partial_t f(y)\;\textrm{d} y&=\int_V\sum_{t\in E(z)}(f((z+tv)^-)-f((z+tv)^+))\;\textrm{d} z\\[5pt] &=\sum_{J\subset I}\int_{V_J} \sum_{t\in E(z)}(f((z+tv)^-)-f((z+tv)^+))\;\textrm{d} z\\[5pt] &=\sum_{J\subset I}\int_{V_J} \sum_{(i,s)\in J}(f(F_i^s(z)^-)-f(F_i^s(z)^+))\;\textrm{d} z\\[5pt] &=\sum_{(i,s)\in I}\sum_{J\ni (i,s)} \int_{V_J}(f(F_i^s(z)^-)-f(F_i^s(z)^+))\;\textrm{d} z\\[5pt] &=\sum_{(i,s)\in I} \int_{V_i^s\cap V}(f(F_i^s(z)^-)-f(F_i^s(z)^+))\;\textrm{d} z.\end{align*}

Since the differential of each $\phi_i$ is always one-to-one and since $D\phi_i(x)(u).u$ is never orthogonal to v for $x\in W_i^s $ and $u\neq 0$ , the local inverse function theorem implies that the maps $F_i^s$ are $C^1$ . Furthermore, the image of $F_i^s$ is $M_i^s$ and the normal to $M_i^s$ at $y=F_i^s(z)$ is $n(y)=\pm n_i((p_v\circ \phi)^{-1}(z))$ . Therefore,

\begin{align*}\int_{V_i^s\cap V}(f(F_i^s(z)^-)-f(F_i^s(z)^+))\;\textrm{d} z=\int_{M_i^s\setminus\mathcal N_v}(f(y^-)-f(y^+))|n(y)|\;\textrm{d}\sigma(y).\end{align*}

Finally, since $\partial U\setminus \mathcal N_v=\cup_{(i,s)\in I}(M_i^s\setminus \mathcal N_v)$ ,

\begin{align*} \sum_{(i,s)\in I} \int_{V_i^s\cap V}(f(F_i^s(z)^-)-f(F_i^s(z)^+))\;\textrm{d} z=\int_{\partial U\setminus \mathcal N_v}(f(y^-)-f(y^+))|n(y)|\;\textrm{d}\sigma(y).\end{align*}

D.1. Bouncy particle sampler: proof of Proposition 5

For the specified $Q_{\mathrm{BPS}}$ , we first derive the form of the associated probability kernel on velocities, $Q'_{\!\!x}$ , remembering that $Q'_{\!\!x}$ is obtained by flipping the velocity and then applying the transition defined by $Q_{\mathrm{BPS}}$ . This gives

\begin{align*}Q'_{\!\!x}(v^{\prime}|v)=\left\{\begin{array}{c@{\quad}l}1, &v^{\prime}= -v \mbox{ and } v\in\mathcal{V}_x^-,\\[5pt] \pi_{k_1(x)}(x)/\pi_{k_2(x)}(x), & v^{\prime}= -v \mbox{ and } v\in\mathcal{V}_x^+,\\[5pt] 1-\pi_{k_1(x)}(x)/\pi_{k_2(x)}(x), & v^{\prime}= -v + 2 \langle n,v\rangle n \mbox{ and } v\in\mathcal{V}_x^+.\end{array}\right.\end{align*}

The transition kernel $Q'_{\!\!x}$ allows for two possible transitions, either $v^{\prime}=- v$ or $v^{\prime}=v - 2 \langle n,v\rangle n$ . These transitions have the following properties:

1. (i) In the first case, if $v \in\mathcal{V}_x^+$ then $v^{\prime} \in\mathcal{V}_x^-$ , and vice versa. In the second case, if $v\in \mathcal{V}_x^+$ then $v^{\prime}\in \mathcal{V}_x^+$ .

2. (ii) For either transition, $\langle v^{\prime},v^{\prime}\rangle=\langle v,v \rangle$ , and $|\langle n,v^{\prime}\rangle|=|\langle n,v\rangle|$ . Furthermore, by the spherical symmetry of p(v) for the bouncy particle sampler, the first of these means that $p(v) = p(v^{\prime})$ .

We need to show that $l_x(dv^{\prime})=\int Q'_{\!\!x}(dv^{\prime}|v)l_x(v) \;\textrm{d} v$ , where

\begin{align*} l_x(v) =\left\{ \begin{array}{c@{\quad}l} |\langle n,v \rangle| p(v) \pi_{k_2(x)}(x) & \forall v\in \mathcal{V}_x^+,\\[5pt] |\langle n,v \rangle| p(v) \pi_{k_1(x)}(x) & \forall v\in \mathcal{V}_x^-. \end{array} \right.\end{align*}

We will show this holds first for $v^{\prime}\in \mathcal{V}_x^+$ and then for $v^{\prime}\in \mathcal{V}_x^-$ .

If $v^{\prime}\in \mathcal{V}_x^+$ then there are two possible transitions, from $v=-v^{\prime} \in \mathcal{V}_x^-$ and from $v\in \mathcal{V}_x^+$ where $v=-v'+2\langle v',n \rangle n$ . Let $v^*=-v'+2\langle v',n \rangle n$ . Since the Jacobian of both transformations is 1,

\begin{eqnarray*}{{\int Q'_{\!\!x}(dv^{\prime}|v)l_x(v) \;\textrm{d} v}= 1\cdot\left( |\langle n,-v^{\prime} \rangle| p(\!-\!v^{\prime}) \pi_{k_1(x)}(x)\right) {\;\textrm{d} v^{\prime}}}\\[5pt] & & +\left(1-\frac{\pi_{k_1(x)}(x)}{\pi_{k_2(x)}(x)}\right)\left(|\langle n,v^* \rangle| p(v^*) \pi_{k_2(x)}(x)\right) {\;\textrm{d} v^{\prime}},\\[5pt] &=& |\langle n,v^{\prime} \rangle| p(v^{\prime}) \left(\pi_{k_1(x)}(x) + \left(1-\frac{\pi_{k_1(x)}(x)}{\pi_{k_2(x)}(x)}\right) \pi_{k_2(x)}(x) \right) {\;\textrm{d} v^{\prime}},\end{eqnarray*}

where the last equality comes from applying Property (ii) of the transition. The last expression simplifies to $l_x(\textrm{d} v^{\prime})$ as required.

For $v^{\prime}\in \mathcal{V}_x^-$ we have only one transition and thus

\begin{align*}{\int Q'_{\!\!x}(dv^{\prime}|v)l_x(v) dv}= \left(\frac{\pi_{k_1(x)}(x)}{\pi_{k_2(x)}(x)} \right) \left( |\langle n,-v^{\prime} \rangle| p(\!-\!v^{\prime}) \pi_{k_2(x)}(x)\right) {\;\textrm{d} v^{\prime}},\end{align*}

which, again using Property (ii), is $l_x(v^{\prime})$ .

D.2. Coordinate sampler: proof of Proposition 6

We follow a similar argument to that of the previous section. First we write down the form of $Q'_{\!\!x}$ derived from $Q_{\mathrm{CS}}$ :

\begin{align*}Q'_{\!\!x}(v^{\prime}|v)=\left\{\begin{array}{c@{\quad}l}1, &v^{\prime}= -v \mbox{ and } v\in\mathcal{V}_x^-,\\[5pt] \pi_{k_1(x)}(x)/\pi_{k_2(x)}(x), & v^{\prime}= -v \mbox{ and } v\in\mathcal{V}_x^+,\\[5pt] (1-\pi_{k_1(x)}(x)/\pi_{k_2(x)}(x))\frac{\langle v^{\prime},n\rangle}{K}, & v^{\prime} \in\mathcal{V}_x^+ \mbox{ and } v\in\mathcal{V}_x^+,\end{array}\right.\end{align*}

where $K=\sum_{v\in\mathcal{V}_x^+} |\langle n,v \rangle|$ .

For the coordinate sampler, $p(v)=1/(2d)$ for each of the possible values for v. Now we need to show $l_x(v^{\prime})=\sum Q'_{\!\!x}(v^{\prime}|v)l_x(v)$ . We will consider the cases $v^{\prime}\in \mathcal{V}_x^+$ and $v^{\prime}\in\mathcal{V}_x^-$ separately. For the latter case, the argument is the same as for the bouncy particle sampler. Thus we just present the case for $v^{\prime}\in \mathcal{V}_x^+$ .

We have

\begin{eqnarray*}{\sum Q'_{\!\!x}(v^{\prime}|v)l_x(v)= 1\cdot\left( |\langle n,-v^{\prime} \rangle| p(\!-\!v^{\prime}) \pi_{k_1(x)}(x)\right)}\\[5pt] & & +\;\langle v^{\prime},n \rangle \left(1-\frac{\pi_{k_1(x)}(x)}{\pi_{k_2(x)}(x)}\right)\sum_{v\in\mathcal{V}_x^+} \frac{|\langle n,v \rangle|}{K} p(v) \pi_{k_2(x)}(x)\\[5pt] &=& \;p(v^{\prime}) \langle v^{\prime},n \rangle\left( \pi_{k_1(x)}(x) + \left(1-\frac{\pi_{k_1(x)}(x)}{\pi_{k_2(x)}(x)}\right)\pi_{k_2(x)}(x)\sum_{v\in\mathcal{V}_x^+} \frac{|\langle n,v \rangle|}{K}\right) \\[5pt] &= & \;p(v^{\prime}) \langle v^{\prime},n \rangle \left( \pi_{k_1(x)}(x) + \left(1-\frac{\pi_{k_1(x)}(x)}{\pi_{k_2(x)}(x)}\right)\pi_{k_2(x)}(x)\right).\end{eqnarray*}

The second equality comes from the fact that p(v) is constant for all $v\in\mathcal{V}$ . The third equality comes from the definition of K.

D.3. Zigzag sampler: proof of Proposition 7

As discussed, the kernel for the velocity does not depend on k. Without loss of generality, we can set $k=1$ for implementing the algorithm that defines $Q_{\mathrm{ZZ}}$ . We need to show that $Q'_{\!\!x}$ keeps $l_x$ invariant. We will prove this by showing that the following stronger detailed balance condition holds:

\begin{align*}l_x(v)Q'_{\!\!x}(v^{\prime}|v)=l_x(v^{\prime})Q'_{\!\!x}(v|v^{\prime}), \quad\forall v,v^{\prime}\in \mathcal{V}.\end{align*}

As $Q'_{\!\!x}(v^{\prime}|v)=Q_{\mathrm{ZZ}}(v^{\prime}|-v)$ , writing the detailed balance condition for pairs $-v$ and v ^′, we have that it suffices to show

\begin{align*}l_x(\!-\!v)Q_{\mathrm{ZZ}}(v^{\prime}|v)=l_x(v^{\prime})Q_{\mathrm{ZZ}}(\!-\!v|-v^{\prime}), \quad\forall v,v^{\prime}\in \mathcal{V}.\end{align*}

By a slight abuse of notation, let $k(v)=k_1(x)$ if $v\in\mathcal{V}^-_x$ and $k(v)=k_2(x)$ if $v\in\mathcal{V}^+_x$ . Then we can write $l_x(v)=|\langle n,v \rangle|p(v)\pi_{k(v)}(x)$ . Thus, using the fact that p(v) defines a uniform distribution on $\mathcal{V}$ , we have that the detailed balance condition simplifies to

\begin{align*}|\langle n,v \rangle| \pi_{k(\!-\!v)}(x)Q_{\mathrm{ZZ}}(v^{\prime}|v)=|\langle n,v^{\prime} \rangle| \pi_{k(v^{\prime})}(x)Q_{\mathrm{ZZ}}(\!-\!v|-v^{\prime}), \quad\forall v,v^{\prime}\in \mathcal{V}.\end{align*}

This can be viewed as matching the probability that we have a velocity v and transition to v ^′ with one where we flip the velocities, starting at $-v^{\prime}$ and transitioning to $-v$ .

We show that the detailed balance condition holds separately for different combinations of whether $v\in\mathcal{V}_x^+$ or $v\in\mathcal{V}_x^-$ and whether $v^{\prime}\in\mathcal{V}_x^+$ or $v^{\prime}\in\mathcal{V}_x^-$ .

First assume $v\in\mathcal{V}_x^+$ and $v^{\prime}\in\mathcal{V}_x^-$ . This corresponds to a trajectory that is moving from the lower- to the higher-density region, but that reflects off the boundary and stays in the lower-density region. It is straightforward to see that the events that change the velocity only increase $\langle n,v \rangle$ , the speed at which the trajectory moves through the boundary region to the higher-density region. Thus a transition from $\mathcal{V}_x^+$ to $\mathcal{V}_x^-$ is impossible, and $Q_{\mathrm{ZZ}}(v^{\prime}|v)=0$ . Similarly, $-v^{\prime}\in\mathcal{V}_x^+$ and $-v\in\mathcal{V}_x^-$ so $Q_{\mathrm{ZZ}}(\!-\!v|-v^{\prime})=0$ . Hence the detailed balance conditions trivially hold in this case.

Next assume $v\in\mathcal{V}_x^-$ and $v^{\prime}\in\mathcal{V}_x^+$ . This corresponds to a trajectory that is moving from the higher- to the lower-density region, but that reflects off the boundary and stays in the higher-density region. In this case $k(\!-\!v)=k(v^{\prime})$ and thus the detailed balance condition becomes

(11) \begin{equation} |\langle n, v\rangle| Q_{\mathrm{ZZ}}(v^{\prime}|v) =|\langle n, v^{\prime} \rangle| Q_{\mathrm{ZZ}}(\!-\!v|-v^{\prime}), \quad \forall v\in\mathcal{V}_x^-,\;v^{\prime}\in\mathcal{V}_x^+.\end{equation}

To prove that the detailed balance condition holds, we will first obtain an expression for $Q_{\mathrm{ZZ}}(v^{\prime}|v)$ , and then introduce a coupling between a transition for v to v ^′ and one from $-v^{\prime}$ to $-v$ to link it to a similar expression for $Q_{\mathrm{ZZ}}(\!-\!v|-v^{\prime})$ .

The randomness in the algorithm that defines $Q_{\mathrm{ZZ}}$ only comes through the randomness of the event times simulated in Step (a) of Section 6.2.3. Remember that $\tau_i$ is the time at which component i of the velocity will switch, if the trajectory is still within the boundary region. Each $\tau_i$ is (conditionally) independent of the others, and has an exponential distribution with rate $\max\{0, -n_iv_i\}$ , where $n_i$ is the component of the ith coordinate of the unit normal n. If $n_iv_i\geq0$ , then $\tau_i=\infty$ .

It is helpful to introduce three sets of components:

• • Let $\mathcal{S}_1$ be the set of components i such that $v^{\prime}_i=v_i$ and $n_iv_i<0$ .

• • Let $\mathcal{S}_2$ be the set of components i such that $v^{\prime}_i=-v_i$ .

• • Let $\mathcal{S}_3$ be the set of components i such that $v^{\prime}_i=v_i$ and $n_iv_i\geq0$ .

So $\mathcal{S}_1$ is the set of components of the velocity v that are moving the particle towards the low-density region, and are unchanged by the transition to v ^′; $\mathcal{S}_2$ is the set of components that flip during the transition from v to v ^′; and $\mathcal{S}_3$ is the set of components of the velocity v that are moving the particle towards the high-density region, and are unchanged by the transition to v ^′.

Only components i of the velocity for which $n_iv_i<0$ can change during the transition from v to v ^′. This means that if there exists $i\in\mathcal{S}_2$ such that $n_iv_i \geq 0$ , then the transition from v to v ^′ is impossible. By the same argument, the transition from $-v^{\prime}$ to $-v$ is impossible. Thus in this case $Q_{\mathrm{ZZ}}(v^{\prime}|v)=Q_{\mathrm{ZZ}}(\!-\!v|-v^{\prime})=0$ and detailed balance trivially holds. So in the following we will assume that $n_iv_i<0$ for $i \in\mathcal{S}_2$ .

By a similar argument we have that the set $\mathcal{S}_1$ is the set of indices of the velocity that could have changed during the transition from v to v ^′, but did not, whereas $\mathcal{S}_3$ is the set of indices of the velocity that could never have changed during the transition.

To ease notation, let $m=|\mathcal{S}_2|$ , the number of indices in the set $\mathcal{S}_2$ , and note that $m\geq 1$ as $v\neq v^{\prime}$ . Without loss of generality we can relabel the coordinates so that $\mathcal{S}_2=\{1,\ldots,m\}$ , and we will use $\tau_{1:m}$ to denote the vector of event times for the coordinates in $\mathcal{S}_2$ .

We now introduce a function of time, t, that depends on $\tau_{1:n}$ . This is

\begin{align*}h(t;\;\tau_{1:n})=\sum_{i=1}^{m} n_iv_i(t - {2}\max\{0,t-\tau_i\})+\sum_{i \in \mathcal{S}_1} n_iv_i t+ \sum_{i \in \mathcal{S}_3} n_iv_i t.\end{align*}

This can be viewed as the net distance travelled by the trajectory up to time t in the direction of the normal n, given that only velocity coordinates in $\mathcal{S}_2$ can change, and these change at times $\tau_{1:m}$ . This function is important as it determines when the trajectory leaves the boundary region, and determines the termination of the simulation algorithm in Step (b). As $v\in\mathcal{V}_x^-$ and $v^{\prime}\in\mathcal{V}_x^+$ , and the changes in velocity in the direction of n are monotone as we flip components, we have that $h(t;\;\tau_{1:n})$ is strictly decreasing at $t=0$ , is strictly increasing for large enough t, and is unimodal. As $h(0;\tau_{1:m})=0$ , this means that there is a unique $t^*(\tau_{1:m})>0$ such that $h(t^*(\tau_{1:m});\tau_{1:m})=0$ . This is the exit time from the boundary region calculated in Step (b) of the algorithm.

We can now define the set $\mathcal{T}$ of values of $\tau_{1:m}$ that are consistent with a transition from v to v ^′. The conditions are that all components of the velocity must flip before $t^*$ , and that the trajectory must not pass through the boundary region—see the other stopping criteria in Step (b) of the algorithm. This gives us that

\begin{align*}\mathcal{T}=\left\{\tau_{1:m} \;:\; \tau_i \leq t^*(\tau_{1:m}),\;i=1,\ldots,m; \;\;\min_{0<t<t^*(\tau_{1:m})} h(t;\;\tau_{1:n}) >-C\right\}.\end{align*}

The probability of a transition from v to v ^′ is thus the probability that $\tau_{1:m}\in\mathcal{T}$ times the probability that $\tau_i>t^*(\tau_{1:m})$ for $i\in\mathcal{S}_1$ . As each $\tau_i$ , $i \in \mathcal{S}_1$ or $i\in\mathcal{S}_2$ , has an independent exponential distribution with rate $-n_iv_i$ ,

\begin{align*}Q_{\mathrm{ZZ}}(v^{\prime}|v) =\int_\mathcal{T} \left(\prod_{i\in\mathcal{S}_1} \exp\{n_iv_it^*(\tau_{1:m})\} \right)\left(\prod_{i=1}^m (\!-\!n_iv_i)\exp\{n_iv_i\tau_i\}\;\textrm{d} \tau_{1:m}\right).\end{align*}

Now consider the reverse transition, from $-v^{\prime}$ to $-v$ . Under our existing definitions of $\mathcal{S}_1$ , $\mathcal{S}_2$ , and $\mathcal{S}_3$ , we have that $\mathcal{S}_2$ is still the set of indices that the flip for the transition from $-v^{\prime}$ to $-v$ , but now $\mathcal{S}_1$ is the set of components of the velocity that could never have flipped, while $\mathcal{S}_3$ is the set of components that could have flipped but did not.

We can define the same quantities for the reverse transition from $-v^{\prime}$ to $-v$ . We will use tildes to denote quantities that relate to this transition. So $\tilde{\tau}_{1:m}$ will be the vector of flip times for components in $i\in\mathcal{S}_2$ . We have

\begin{align*}\tilde{h}(t;\;\tilde{\tau}_{1:n})=\sum_{i=1}^{m} n_iv_i(t - {2}\max\{0,t-\tilde{\tau}_i\})-\sum_{i \in \mathcal{S}_1} n_iv_i t- \sum_{i \in \mathcal{S}_3} n_iv_i t,\end{align*}

using the fact that $-v^{\prime}_i=v_i$ for $i \in \mathcal{S}_2$ and $-v^{\prime}_i=-v_i$ otherwise. By the same argument as above, there is a unique $\tilde{t}^*(\tilde{\tau}_{1:m})>0$ such that $\tilde{h}(\tilde{t}^*(\tilde{\tau}_{1:m});\tilde{\tau}_{1:m})=0$ . The set of possible values of $\tilde{\tau}_{1:m}$ that are consistent with the transition from $-v^{\prime}$ to v is

\begin{align*}\tilde{\mathcal{T}}=\left\{\tilde{\tau}_{1:m}\; :\; \tilde{\tau}_i \leq \tilde{t}^*(\tilde{\tau}_{1:m}),\;i=1,\ldots,m; \;\;\min_{0<t<\tilde{t}^*(\tilde{\tau}_{1:m})} \tilde{h}(t;\;\tilde{\tau}_{1:n}) > -C\right\}.\end{align*}

Finally we can write down the transition probability as before, remembering that the rate of flipping for components $i\in\mathcal{S}_2$ is $-n_iv_i$ as before, but for $i\in\mathcal{S}_3$ it is $n_iv_i$ . Thus

(12) \begin{eqnarray}{Q_{\mathrm{ZZ}}(\!-\!v|-v^{\prime}) =} \nonumber\\[5pt] & &\left(\prod_{i\in\mathcal{S}_3} \exp\{-n_iv_i\tilde{t}^*(\tilde{\tau}_{1:m})\} \right)\left(\int_{\tilde{\mathcal{T}}}\prod_{i=1}^m (\!-\!n_iv_i)\exp\{n_iv_i\tilde{\tau}_i\}\;\textrm{d} \tilde{\tau}_{1:m}\right). \end{eqnarray}

To relate the two transition probabilities, we introduce a coupling between $\tau_{1:m}$ and $\tilde{\tau}_{1:m}$ , so $\tilde{\tau}_{1:m}=g(\tau_{1:m})$ , where

\begin{align*}\tilde{\tau}_i=g(\tau_{1:m})_i=t^*(\tau_{1:m})-\tau_i.\end{align*}

This coupling is a natural one. If we consider the path through the boundary region given by $\tau_{1:m}$ that transitions from v to v ^′, we can reverse that path to get a path that transitions from $-v^{\prime}$ to $-v$ . For the forward path a flip of component i at time $\tau_i$ occurs at a time $t^*(\tau_{1:m})-\tau_i$ prior to the end of the path. Thus for the reverse path the flip would occur at time $t^*(\tau_{1:m})-\tau_i$ .

It is straightforward to show that if $\tilde{\tau}_{1:m}=g(\tau_{1:m})$ then $h(t;\;\tau_{1:m})=\tilde{h}(t^*(\tau_{1:m})-t;\;\tilde{\tau}_{1:m})$ . This result is intuitive; it is saying the distance of the forward trajectory within the boundary region at time t is equal to the distance of the backward trajectory within the boundary region at time $t^*(\tau_{1:m})-t$ . This immediately implies that $t^*(\tau_{1:m})=\tilde{t}^*(\tilde{\tau}_{1:m})$ : the exit times for the forward and backward trajectories are the same. Furthermore, if we consider the second constraint on $\tau_{1:m}$ in the definition of $\mathcal{T}$ , then we have

\begin{align*}\min_{0<t<t^*(\tau_{1:m})} h(t;\;\tau_{1:n}) =\min_{0<t<\tilde{t}^*(\tilde{\tau}_{1:m})} \tilde{h}(t;\;\tilde{\tau}_{1:n}),\end{align*}

for $\tilde{\tau}_{1:m}=g(\tau_{1:m})$ . Combining this with the fact that $\tau_i\leq t^*(\tau_{1:m})$ , we have $\tilde{\tau}_{1:m}\leq \tilde{t}^*(\tilde{\tau}_{1:m})$ . We have that the function g maps $\tau_{1:m}\in\mathcal{T}$ to $\tilde{\tau}_{1:m}\in\tilde{\mathcal{T}}$ . Furthermore, the function g is invertible, and by similar arguments we have that $g^{-1}$ maps $\tilde{\tau}_{1:m}\in\tilde{\mathcal{T}}$ to $\tau_{1:m}\in\mathcal{T}$ . Hence g is a bijection from $\mathcal{T}$ to $\tilde{\mathcal{T}}$ .

The function g defines a linear map between $\tau_{1:m}$ and $\tilde{\tau}_{1:m}$ . For $\tau_{1:m}\in\mathcal{T}$ we have that, by definition of $t^*(\tau_{1:m})$ ,

(13) \begin{equation} \sum_{i=1}^m n_iv_i(2\tau_i-t^*(\tau_{1:m}))+\sum_{i=m+1}^dn_iv_it^*(\tau_{1:m})=0.\end{equation}

This gives that

\begin{align*}t^*(\tau_{1:m})=\sum_{i=1}^m \left( \frac{-2n_iv_i}{K} \right)\tau_i, \mbox{ where } K= -\sum_{i=1}^m n_iv_i+\sum_{i=m+1}^d n_iv_i.\end{align*}

Furthermore, using that v ^′ is equal to v except that $v_i$ is flipped for $i=1,\ldots,m$ , we have $K=\langle v^{\prime},n\rangle $ .

Let $b_{1:m}$ be the $1\times m$ vector whose ith entry is $b_i=2v_in_i/K$ . If we let $1_m$ denote the $1\times m$ vector of ones, and $I_m$ the $m\times m$ identity matrix, then we have

\begin{align*}\tilde{\tau}_{1:m} = g(\tau_{1:m}) = (b_{1:m}1_{1:m}^\top-I_m)\tau_{1:m}=A\tau_{1:m},\end{align*}

where the $m \times m$ matrix $A=(b_{1:m}1_{1:m}^\top-I_m)$ . In the following argument we will make the change of variables $\tilde{\tau}_{1:m}=g(\tau_{1:m})=A\tau_{1:m}$ , and we will need the determinant of the Jacobian of this transformation. Using the matrix determinant lemma, this is given by

\begin{align*}|\mbox{det}(A)|=|(1-1_{1:m}^\top b_{1:m})||\mbox{det}(\!-\!I_m)|=\left|1-\sum_{i=1}^m b_i \right|.\end{align*}

This simplifies to

\begin{align*}\left|1-\sum_{i=1}^m b_i \right|=\left|1+\sum_{i=1}^m \frac{2v_in_i}{\langle v^{\prime},n\rangle}\right|=\left|\frac{\langle v^{\prime},n\rangle+\sum_{i=1}^m 2v_in_i}{\langle v^{\prime},n\rangle}\right|=\left|\frac{\langle v,n\rangle}{\langle v^{\prime},n\rangle}\right|.\end{align*}

So now, taking the definition of $Q_{\mathrm{ZZ}}(v^{\prime}|v)$ and applying the change of variables $\tilde{\tau}_{1:m}=g(\tau_{1:m})$ , we get

\begin{eqnarray*}{Q_{\mathrm{ZZ}}(v^{\prime}|v) = \int_\mathcal{T} \prod_{i\in\mathcal{S}_1} \exp\{n_iv_it^*(\tau_{1:m})\}\prod_{i=1}^m \left((\!-\!n_iv_i)\exp\{n_iv_i\tau_i\}\right)\;\textrm{d} \tau_{1:m}}\\[5pt] &=& \;\int_{\tilde{\mathcal{T}}}\exp\!\left\{\sum_{i\in\mathcal{S}_1} n_iv_i \tilde{t}^*(\tilde{\tau}_{1:m})\right\}\left|\frac{\langle v,n\rangle}{\langle v^{\prime},n\rangle}\right|\prod_{i=1}^m \left((\!-\!n_iv_i)\exp\{n_iv_i(\tilde{t}^*(\tilde{\tau}_{1:m})-\tilde{\tau}_i)\}\right)\;\textrm{d} \tilde{\tau}_{1:m} \\[5pt] &=&\;\left|\frac{\langle v,n\rangle}{\langle v^{\prime},n\rangle}\right|\int_{\tilde{\mathcal{T}}}\left(\prod_{i=1}^m (\!-\!n_iv_i)\right)\exp\!\left\{\sum_{i\in\mathcal{S}_1} n_iv_i \tilde{t}^*(\tilde{\tau}_{1:m})+ \sum_{i=1}^m n_iv_i(\tilde{t}^*(\tilde{\tau}_{1:m})-\tilde{\tau}_i)\right\}\;\textrm{d} \tilde{\tau}_{1:m} \\[5pt] &=&\;\left|\frac{\langle v,n\rangle}{\langle v^{\prime},n\rangle}\right|\int_{\tilde{\mathcal{T}}}\left(\prod_{i=1}^m (\!-\!n_iv_i)\right)\exp\!\left\{\sum_{i\in\mathcal{S}_3} -n_iv_i \tilde{t}^*(\tilde{\tau}_{1:m})+\sum_{i=1}^m n_iv_i \tilde{\tau}_i \right\}\;\textrm{d} \tilde{\tau}_{1:m},\end{eqnarray*}

where the final equality comes from the definition of $\tilde{t}^*(\tilde{\tau}_{1:m})=t^*(\tau_{1:m})$ , using (13) after substituting in $\tau_i=\tilde{t}^*(\tilde{\tau}_{1:m})-\tilde{\tau}_i$ .

By comparing the final expression with (12), we get that

\begin{align*}Q_{\mathrm{ZZ}}(v^{\prime}|v) = \left|\frac{\langle v,n\rangle}{\langle v^{\prime},n\rangle}\right| Q_{\mathrm{ZZ}}(\!-\!v|v^{\prime}),\end{align*}

which satisfies (11) as required.

The final combination involves $v,v^{\prime}\in\mathcal{V}_x^+$ and $-v^{\prime},-v\in\mathcal{V}_x^-$ , or vice versa. The detailed balance condition in this case becomes

\begin{align*}\pi_{k_1(x)}(x)|\langle n, v\rangle| Q_{\mathrm{ZZ}}(v^{\prime}|v) =\pi_{k_2(x)}(x)|\langle n, v^{\prime} \rangle| Q_{\mathrm{ZZ}}(\!-\!v|-v^{\prime}), \quad \forall v,v^{\prime}\in\mathcal{V}_x^-.\end{align*}

We can show this using an argument similar to the one above, with the same coupling of paths from v to v ^′ with paths from v ^′ to v. The main differences are as follows. First, the definition of $\mathcal{T}$ is simplified to

\begin{align*}\left\{\tau_{1:m} \;:\; \tau_i \leq t^*(\tau_{1:m}),\;i=1,\ldots,m\right\},\end{align*}

as, by monotonicity of the changes in velocity, we do not need to check whether the other exit condition in Step (b) holds. Second, the definition of $t^*(\tau_{1:m})$ changes, with it being the value of t for which $h(t;\;\tau_{1:m})=C$ . For $\tau_{1:m}\in\mathcal{T}$ , this becomes

\begin{align*}\sum_{i=1}^m n_iv_i(2\tau_i-t^*(\tau_{1:m}))+\sum_{i=m+1}^d n_iv_it^*(\tau_{1:m})=C\end{align*}

because of the different exit condition in Step (b). We have similar changes to the definitions of $\tilde{\mathcal{T}}$ and $\tilde{t}^*(\tilde{\tau}_{1:m})$ .

However, we can define $Q_{\mathrm{ZZ}}(v^{\prime}|v)$ and $Q_{\mathrm{ZZ}}(\!-\!v|-v^{\prime})$ in a similar way. Furthermore we can use the same linear transformation g, which is still a bijection between $\mathcal{T}$ and $\tilde{\mathcal{T}}$ . Whilst the definition $t^*$ has changed, this only introduces an additive constant into the linear transformation defined by g, and thus the Jacobian of the transformation is unchanged. Following the argument above, we thus get to the same expression for $Q_{\mathrm{ZZ}}(v^{\prime}|v)$ after making the change of variables:

\begin{eqnarray*}{Q_{\mathrm{ZZ}}(v^{\prime}|v)=\left|\frac{\langle v,n\rangle}{\langle v^{\prime},n\rangle}\right|\int_{\tilde{\mathcal{T}}}\left(\prod_{i=1}^m (\!-\!n_iv_i)\right) }\\[5pt] &\times&\exp\!\left\{\sum_{i\in\mathcal{S}_1} n_iv_i \tilde{t}^*(\tilde{\tau}_{1:m})+ \sum_{i=1}^m n_iv_i(\tilde{t}^*(\tilde{\tau}_{1:m})-\tilde{\tau}_i)\right\}\;\textrm{d} \tilde{\tau}_{1:m}.\end{eqnarray*}

Now substituting in our new definition of $\tilde{t}^*(\tilde{\tau}_{1:m})=t^*(\tau_{1:m})$ we get

\begin{align*}Q_{\mathrm{ZZ}}(v^{\prime}|v)=\left|\frac{\langle v,n\rangle}{\langle v^{\prime},n\rangle}\right| Q_{\mathrm{ZZ}} (\!-\!v|v^{\prime}) \exp\{C\},\end{align*}

where the additional factor of $\exp\{C\}$ is due to the different definition of $t^*$ . As $C=\log(\pi_{k_2(x)}(x)/\pi_{k_1(x)}(x) )$ this rearranges to

\begin{align*}\pi_{k_1(x)}(x)|\langle v^{\prime},n\rangle|Q_{\mathrm{ZZ}}(v^{\prime}|v) =\pi_{k_2(x)}(x)|\langle v,n\rangle|Q_{\mathrm{ZZ}}(\!-\!v|v^{\prime}),\end{align*}

as required.

D.4. Bouncy particle sampler with reflect and refract: proof of Proposition 8

We first write down the form of $Q_x'$ derived from $Q_{\text{RR}}$ . To do this it is helpful to define functions that define the new velocity at a reflection or refraction. Remember that $\Delta U$ is the change in log density as we move from the region with lower density, $k_1$ , to the region with higher density, $k_2$ , and is strictly positive. We will define the reflection of a velocity v by $L(v)=v - 2 \langle v,n\rangle n$ . For refractions we need to distinguish between transitions from $k_1$ to $k_2$ and vice versa. The refraction for a transition from $k_1$ to $k_2$ is $R_1(v) = \sqrt{\langle v,n\rangle^2 + 2 \Delta U}n + v - \langle v,n\rangle n$ , and for $k_2$ to $k_1$ it is $R_2(v)= \sqrt{\langle v,n\rangle^2 - 2 \Delta U}n + v - \langle v,n\rangle n$ , with $R_2$ being defined only for v such that $\langle v,n\rangle^2 > 2 \Delta U $ .

Using these functions, we can write the transition kernel $Q_x'$ as

\begin{align*}Q_x'(v^{\prime}|v)= \left\{ \begin{array}{c@{\quad}l}1 & \mbox{ for } v\in\mathcal{V}_x^- \mbox{ and } v^{\prime}=R_1(\!-\!v),\\[5pt] 1 & \mbox{ for } v\in\mathcal{V}_x^+,\; \langle v,n\rangle^2 > 2 \Delta U,\mbox{ and } v^{\prime}=R_2(\!-\!v),\\[5pt] 1 & \mbox{ for } v\in\mathcal{V}_x^+,\; \langle v,n\rangle^2 \leq 2 \Delta U, \mbox{ and } v^{\prime}=L(\!-\!v),\\[5pt] 0 & \mbox{otherwise}.\end{array}\right.\end{align*}

We need to show that $l_x(dv^{\prime}) = \int Q'_{\!\!x}(dv^{\prime}|v) l_x(v) dv$ with

\begin{align*} l_x(v) =\left\{ \begin{array}{c@{\quad}l} |\langle n,v \rangle| p(v) \pi_{k_2(x)}(x) & \forall v\in \mathcal{V}_x^+,\\[5pt] |\langle n,v \rangle| p(v) \pi_{k_1(x)}(x) & \forall v\in \mathcal{V}_x^-. \end{array} \right.\end{align*}

Since, for each v ^′, there is only a finite number of v such that $Q'_{\!\!x}(v^{\prime}|v) > 0$ , the former is equivalent to $l_x(v^{\prime})=\sum Q'_{\!\!x}(v^{\prime}|v)l_x(v) J_{v,v^{\prime}}$ , where $J_{v,v^{\prime}}$ is the Jacobian of the deterministic transformation from v to v ^′.

If $v^{\prime}\in\mathcal{V}_x^+$ and $\langle v^{\prime},n\rangle^2 \leq 2 \Delta U$ , then the only v for which $Q'_{\!\!x}(v^{\prime}|v)$ is non-zero is $v=L(\!-\!v^{\prime})$ . The result trivially holds as this transition has unit Jacobian, $\|v\|=\|v^{\prime}\|$ , so $p(v)=p(v^{\prime})$ ; $|\langle n,v \rangle|= |\langle n,v^{\prime} \rangle|$ ; $v^{\prime}\in\mathcal{V}_x^+$ and $Q'_{\!\!x}(v^{\prime}|v)=1$ .

If $v^{\prime}\in\mathcal{V}_x^+$ and $\langle v,n\rangle^2 > 2 \Delta U$ , then the only v for which $Q'_{\!\!x}(v^{\prime}|v)$ is non-zero is for v such that $v^{\prime}=R_1(v)$ , which is equivalent to $v=R_2(\!-\!v^{\prime})\in\mathcal{V}_x^-$ . So as $p(v)=C_1 \exp\{-\frac{1}{2}\|v\|^2\}$ for some constant $C_1$ ,

(14) \begin{align}\sum Q'_{\!\!x}(v^{\prime}|v)l_x(v)&= l_x(R_2(\!-\!v^{\prime})) \nonumber \\[5pt] &= |\langle n,R_2(\!-\!v^{\prime})\rangle| \cdot p(R_2(\!-\!v^{\prime})) \cdot \pi_{k_1}(x) \nonumber \\[5pt] &= |\langle n,R_2(\!-\!v^{\prime})\rangle| \left(C_1 \exp\!\left\{-\frac{1}{2}\|R_2(\!-\!v^{\prime})\|^2\right\} \right) \pi_{k_1}(x). \end{align}

Now

\begin{align*}\|R_2(\!-\!v^{\prime})\|^2 &=\left\|\sqrt{\langle v^{\prime},n\rangle^2 - 2 \Delta U}n - v^{\prime} + \langle v^{\prime},n\rangle n \right\|^2 =\langle v^{\prime},n\rangle^2 - 2 \Delta U + \|v^{\prime} - \langle v^{\prime},n\rangle n\|^2 \\[5pt] &=\|v^{\prime}\|^2- 2 \Delta U,\end{align*}

where the third equality comes from n being orthogonal to $v^{\prime} - \langle v^{\prime},n\rangle n$ . This gives

\begin{align*}\left(C_1 \exp\!\left\{-\frac{1}{2} (\|v^{\prime}\|^2 - 2\Delta U)\right\} \right) \pi_{k_1}(x) =p(v^{\prime}) \exp\{\Delta U\} \pi_{k_1}(x).\end{align*}

The definition of $\Delta U$ gives that $\exp\{\Delta U\} \pi_{k_1}(x)=\pi_{k_2}(x)$ .

Now $v^{\prime}=\sqrt{\langle v,n\rangle^2 + 2 \Delta U}n - v + \langle v,n\rangle n$ . Up to a sign change, this only changes the velocity in the direction of n. Thus the Jacobian of the transformation is given by

\begin{align*}J_{v,v^{\prime}}^{-1}=\left|\frac{d \langle v^{\prime},n \rangle}{d \langle v,n \rangle}\right| =\left|\frac{d }{d\langle v,n \rangle} \sqrt{\langle v,n\rangle^2 + 2 \Delta U} \right|= \left| \frac{\langle v,n\rangle}{\sqrt{\langle v,n\rangle^2 + 2 \Delta U}}\right| = \frac{|\langle v,n \rangle|}{|\langle v^{\prime},n \rangle|}.\end{align*}

Thus, substituting in these expressions, we have

\begin{align*}\sum Q'_{\!\!x}(v^{\prime}|v)l_x(v) J_{v,v^{\prime}} &= |\langle n,R_2(\!-\!v^{\prime})\rangle|p(v^{\prime}) \pi_{k_2}(x) \frac{|\langle v^{\prime},n \rangle|}{|\langle v,n \rangle|} \\[5pt] &= |\langle n,v\rangle| p(v^{\prime}) \pi_{k_2}(x) \frac{|\langle v^{\prime},n \rangle|}{|\langle v,n \rangle|} \\[5pt] &= |\langle n,v^{\prime}\rangle| p(v^{\prime}) \pi_{k_2}(x) \\[5pt] &= l_x(v^{\prime}),\end{align*}

as required.

An equivalent calculation holds for $v^{\prime}\in\mathcal{V}_x^{-}$ , which is omitted for brevity.