# Dark Energy predictions from GREA: Background and linear perturbation theory

Juan García-Bellido<EMAIL_ADDRESS>Instituto de Física Teórica
UAM-CSIC, Universidad Autonóma de Madrid, Cantoblanco 28049 Madrid, Spain

###### Abstract

General Relativistic Entropic Acceleration (GREA) theory provides a covariant
formalism for out-of-equilibrium phenomena in GR, extending the Einstein
equations with an entropic force that behaves like bulk viscosity with a
negative effective pressure. In particular, the growth of entropy associated
with the homogeneous causal horizon can explain the present acceleration of
the Universe, without introducing a cosmological constant. The dynamics of the
accelerated universe is characterized by a single parameter $\alpha$, the
ratio of the causal horizon to the curvature scale, which provides a unique
history of the Universe distinguishable from that of $\Lambda$CDM. In
particular, we explain the coincidence problem and the Hubble tension by
shifting the coasting point to higher redshifts. All background observables
are correlated among themselves due to their common dependence on $\alpha$.
This scenario gives a specific evolution for the effective equation of state
parameter, $w(a)$. Furthermore, we study the linear growth of matter
perturbations in the context of a homogeneous expanding background driven by
the entropy of the causal horizon. We find that the rate of growth of matter
fluctuations in GREA slows down due to the accelerated expansion and
alleviates the $\sigma_{8}$ tension of $\Lambda$CDM. We compute the growth
function of matter fluctuations, the redshift space distortions in the galaxy
correlation function, as well as the redshift evolution of the baryon acoustic
oscillation scale, and find that the integrated Sachs-Wolfe effect is
significantly larger than in $\Lambda$CDM. It is interesting to note that many
of the tensions and anomalies of the standard model of cosmology are
alleviated by the inclusion of this transient period of acceleration of the
Universe based on known fundamental physics. In the near future we will be
able to constrain this theory with present data from deep galaxy surveys.

††preprint: IFT-UAM/CSIC-24-69

##  I I. Introduction

General Relativistic Entropic Acceleration (GREA) theory Espinosa-Portalés and
García-Bellido (2021) gives a covariant formalism for out of equilibrium
dynamics in the context of general relativity. A consequence of GREA is the
explicit breaking of the time reversal invariance when there is entropy
production. This drives an entropic force that behaves effectively like bulk
viscosity with a negative effective pressure. The natural consequence of this
formalism is that the source for space-time curvature is the Helmholtz free
energy, $F=U-TS$, and not just matter and radiation García-Bellido (2024,
2023).

When the entropy production is associated with the growth of horizons one can
describe their thermodynamical effects as a GHY boundary term Gibbons and
Hawking (1977); York (1972). In the cosmological context, this theory predicts
that the present acceleration of the universe could arise from the growth of
entropy associated with cosmic García-Bellido and Espinosa-Portalés (2021) and
black hole horizons García-Bellido (2024), without the need to introduce a
cosmological constant, providing an alternative dynamics to that of the
Standard Cosmological Model ($\Lambda$CDM).

Already before the recent evidence of significant deviations from the standard
model of cosmology through state-of-the-art observations by DES-Y5 SNe Abbott
et al. (2024) and DESI-Y1 BAO Adame et al. (2024), GREA predictions gave a
better description of cosmological observations than $\Lambda$CDM Arjona et
al. (2022). In particular, GREA explains the coincidence problem and can
alleviate the Hubble tension by shifting the coasting point (when the universe
transitioned from matter domination to acceleration) to higher redshifts and
extending the period of acceleration.

In this paper, we perform a detailed calculation of the background and linear
matter perturbations with the assumption that the entropy growth responsible
for the accelerated expansion comes from the homogeneous cosmological horizon.
We show that, in this case, the dynamics can be described in terms of just one
parameter, $\alpha$, the ratio of the spatial curvature to the horizon
distance, which gives the size of the cosmological horizon and the rate of
growth of the horizon entropy. A given value of $\alpha$ defines a unique
history of the universe, with a dynamics that differs from $\Lambda$CDM and
gives an explicit prediction for all background observables, like the values
of the matter and dark energy content of the universe, the rate of expansion
today, the coasting point, the value of the equation of state of dark energy
and its derivative, the angular diameter and luminosity distances, etc.

In order to compare with the plethora of cosmological data, we need not just
the background evolution but also the linear growth of matter fluctuations. We
derive a second order differential equation in GREA for the amplitude of
matter perturbations in a homogeneous background, and compare with the
corresponding one in $\Lambda$CDM. The solutions match those during the matter
era, but start deviating around redshift one, when the entropic forces begin
to take over and start accelerating the universe.

In Section II we describe the extended Einstein field equations in the
presence of out-of-equilibrium phenomena associated with the growth of
horizons. We recall from Ref. Espinosa-Portalés and García-Bellido (2021) that
the first law of thermodynamics can be incorporated into Einstein’s equations
as an effective bulk viscosity, which drives an effective negative pressure as
a consequence of the second law of thermodynamics. In Section III we explore
the background phenomenology of GREA in the context of entropic forces arising
from the cosmological horizon. In Section IV we derive the linear perturbation
theory of GREA and compare cosmological observables with the corresponding
ones in $\Lambda$CDM. Finally, we conclude in Section V.

##  II II. The Einstein field equations in GREA

We give here a short summary of the covariant formalism developed in Ref.
Espinosa-Portalés and García-Bellido (2021) that describes the extended
Einstein field equations in the presence of out-of-equilibrium phenomenon
associated with the growth of horizons. The gravitational action of the
Standard Model of Cosmology ($\Lambda$CDM) is given by

${\cal S}=\frac{1}{2\kappa}\int_{\cal
M}\\!d^{4}x\sqrt{-g}\Big{[}R+2\kappa{\cal
L}_{m}\Big{]}-\frac{1}{\kappa}\int_{\cal M}\\!d^{4}x\sqrt{-g}\,\Lambda\,,$

where ${\cal L}_{m}$ is the matter Lagrangian, $\kappa=8\pi G$ and $\Lambda$
is the cosmological constant. In the context of GREA we substitute the bulk
term of the cosmological constant for the GHY boundary term of the horizon
García-Bellido and Espinosa-Portalés (2021),

${\cal S}=\frac{1}{2\kappa}\int_{\cal
M}\\!d^{4}x\sqrt{-g}\Big{[}R+2\kappa{\cal
L}_{m}\Big{]}+\frac{1}{\kappa}\int_{\cal H}\\!d^{3}x\sqrt{h}\,K\,,$

where $K$ is the extrinsic curvature and $h$ the 3D metric of the cosmological
horizon, ${\cal H}=\partial{\cal M}$.

The main difference of GREA with respect to $\Lambda$CDM is the breaking of
time-reversal invariance at the level of the action by the entropic force term
that arises due to the growth of the cosmic horizon. This introduces a
cosmological arrow of time. Most of the evolution of the universe has been
quasi-adiabatic since these entropic forces acted only during inflation,
reheating after inflation and at the present time, when the universe is old
enough and the cosmic horizon associated with radiation and matter has grown
sufficiently large for its entropy growth to dominate over the attraction of
matter and induce the present acceleration. Eventually, even this period of
acceleration will end as the entropic term gets diluted and we will end in an
empty Minkowsky space-time García-Bellido and Espinosa-Portalés (2021).

Therefore, the Einstein field equations, extended to out-of-equilibrium
phenomena, can be written as Espinosa-Portalés and García-Bellido (2021),

$G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}R\,g_{\mu\nu}=\kappa\,(T_{\mu\nu}-f_{\mu\nu})\equiv\kappa\,{\cal
T}_{\mu\nu}\,,$ (1)

where $f_{\mu\nu}$ arises from the first law of thermodynamics, and introduces
an effective negative pressure associated with the growth of entropy,
$p_{S}=-TdS/dV<0$, according to the second law of thermodynamics. This extra
component to the Einstein equations can be interpreted as an effective bulk
viscosity term of a real (non-ideal) fluid Espinosa-Portalés and García-
Bellido (2021); Gagnon and Lesgourgues (2011), with
$\Theta=D_{\lambda}u^{\lambda}$ the trace of the congruence of geodesics,

$f_{\mu\nu}=\zeta\,\Theta\,(g_{\mu\nu}+u_{\mu}u_{\nu})=\zeta\,\Theta\,h_{\mu\nu}\,,$
(2)

such that the covariantly-conserved energy-momentum tensor has the form of a
perfect fluid tensor,

$\displaystyle{\cal T}^{\mu\nu}$ $\displaystyle=$ $\displaystyle
p\,g^{\mu\nu}+(\rho+p)u^{\mu}u^{\nu}-\zeta\,\Theta\,h^{\mu\nu}$ (3)
$\displaystyle=$
$\displaystyle\tilde{p}\,g^{\mu\nu}+(\rho+\tilde{p})u^{\mu}u^{\nu}\,,$

with $\tilde{p}=p+p_{S}$, and the bulk viscosity coefficient $\zeta$ can be
written as Espinosa-Portalés and García-Bellido (2021)

$\zeta=\frac{T}{\Theta}\frac{dS}{dV}\,.$ (4)

In the case of an expanding universe, $\Theta=\frac{d}{dt}\ln V=3H$ and the
coefficient becomes $\zeta=T\dot{S}/(9H^{2}a^{3})$, see García-Bellido and
Espinosa-Portalés (2021), with $S$ the entropy per comoving volume of the
Universe. Entropy production therefore implies $\zeta>0$.

Note that the matter contribution to the Hamiltonian constraint arises from
variations of the matter action with respect to the lapse function (and
afterwards setting $N(t)=c$),

${\cal S}_{m}=\int\\!d^{4}x\sqrt{-g}\,{\cal L}_{m}=-\int
dt\,N(t)\,V_{c}\,a^{3}\rho\,,$ (5)

where $V_{c}$ is the comoving volume,

$\displaystyle
V_{c}=\int\\!\frac{4\pi\,r^{2}\,dr}{\sqrt{1-k\,r^{2}}}\,=\,\frac{2\pi}{(-k)^{3/2}}\times$
(6)
$\displaystyle\left[\ln\left(\sqrt{1-k\,r^{2}}-\sqrt{-k\,r^{2}}\right)+\sqrt{-k\,r^{2}(1-k\,r^{2})}\right]\,,$

for a given radial coordinate $r$. The conservation of the total energy
momentum tensor is derived from the first law of thermodynamics,

$\displaystyle dU+pdV-TdS=0\hskip 8.53581pt\Rightarrow$ (7) $\displaystyle
V_{c}\left[d(\rho\,a^{3})+p\,d(a^{3})-T\,dS_{c}\right]=0\hskip
8.53581pt\Rightarrow$ (8)
$\displaystyle\dot{\rho}+3H(\rho+p)=\frac{T\dot{S}_{c}}{a^{3}}\,,$ (9)

where $S_{c}$ is the entropy per comoving volume. For adiabatic expansion, or
constant $S_{c}$, we recover the usual conservation equation Espinosa-Portalés
and García-Bellido (2021).

###  II.1 The causal cosmological horizon

Let us describe here the entropic forces induced by the causal cosmological
horizon of a FLRW universe García-Bellido and Espinosa-Portalés (2021). We
start by considering an arbitrary comoving 2-sphere around the origin of
coordinates. Then the trace of its extrinsic curvature is given by (we set
$c=1$ everywhere)

$\sqrt{h}K=-2N(t)\,r\,a\sqrt{1-kr^{2}}\sin\theta$ (10)

where $r=\sinh(\eta\sqrt{-k})/\sqrt{-k}$ along the lightcone, and $\eta$ is
the conformal time. The boundary term for the causal cosmological horizon,
$d_{H}=a\,\eta$, can be written as García-Bellido and Espinosa-Portalés (2021)

$\displaystyle S_{\rm GHY}$ $\displaystyle=$ $\displaystyle-\frac{1}{2G}\int
dt\,N(t)\,\frac{a}{\sqrt{-k}}\sinh(2\eta\sqrt{-k})$ $\displaystyle=$
$\displaystyle-\int dt\,N(t)\,T_{H}S_{H}\ \equiv-\int
dt\,N(t)\,V_{c}\,a^{3}\rho_{H}\,,$

where $T_{H}$ is the temperature and $S_{H}$ the entropy associated with the
causal cosmological horizon García-Bellido and Espinosa-Portalés (2021)

$k_{\rm B}T_{H}=\frac{\hbar}{2\pi}\
\frac{a\,\sinh(2\eta\sqrt{-k})}{d_{H}^{2}\sqrt{-k}}\,,\hskip
14.22636ptS_{H}=\frac{k_{\rm B}\pi}{\hbar}\frac{d_{H}^{2}}{G}\,.$ (12)

The fact that we can naturally assign a temperature and an entropy to a
hypersurface is a signal of the existence of an underlying quantum description
of gravity and thermodynamics Hawking (1975). This is made explicit by the
appearance of $\hbar$ in both quantities. Their product, however, does not
depend on $\hbar$ and leads to a classical emergent phenomenon, the
acceleration of the universe. Note also that variation w.r.t. the lapse
function will give an extra contribution to the Hamiltonian constraint coming
from the boundary term Gibbons and Hawking (1977), whose origin is related to
the quantum degrees of freedom of the horizon, a phenomenon which has
connections with the Holographic principle ’t Hooft (1993); Susskind (1995).

Figure 1: The horizon distance as a function of the scale factor in GREA
(color lines), compared with $\Lambda$CDM (black dot-dashed line) and matter
dominated (dotted blue line). It is clear that the ratio of the curvature
scale to the horizon distance is very approximately one, for all realizations
of GREA.

Figure 2: The upper panel shows the evolution of the inverse comoving horizon
with the corresponding coasting points (dashed lines) for each value of
$\alpha$ (in color). The lower panel shows the evolution of the effective DE
equation of state $w(a)$.

##  III III. GREA phenomenology

We can then solve the dynamical Friedmann equations for the accelerated
universe in the context of an entropic force arising from the causal horizon.
We will assume here a concrete scenario of open inflation García-Bellido and
Espinosa-Portalés (2021) in which the causal horizon coincides with the
boundary of a bubble wall separating our flat FLRW universe from an open empty
space, therefore matching conditions require that both the scale factor and
its derivative should be continuous at that hypersurface, i.e. $3H_{\rm
in}^{2}=\kappa\,\rho(a_{\rm in})$ and $H_{\rm out}^{2}=-k/a_{\rm in}^{2}$,
which determines $\sqrt{-k}=a_{0}H(a_{0})$ today. We write the Hamiltonian
constraint in conformal time (where primes denote derivatives w.r.t.
$\tau=\sqrt{-k}\,\eta=a_{0}H_{0}\eta$) as 111Note that there is a difference
w.r.t. Ref. García-Bellido and Espinosa-Portalés (2021), and we can now set
$\Omega_{\rm K}=0$ inside the causal patch.

$\displaystyle\left(\frac{a^{\prime}}{a_{0}}\right)^{2}=\Omega_{\rm
M}\left(\frac{a}{a_{0}}\right)+\frac{4\pi}{3}\left(\frac{a}{a_{0}}\right)^{2}\frac{\sinh(2\tau)}{(-k)^{3/2}V_{c}}\,,$
(13)

where the comoving volume (6) is now

$(-k)^{3/2}V_{c}=\pi\left[\sinh(2\sqrt{-k}\,\eta_{0})-2\sqrt{-k}\,\eta_{0}\right]\,.$
(14)

We can then fix $\eta_{0}$ and determine the ratio between the causal horizon
today and the curvature scale ($c=1$),

$\alpha\,H(a_{0})\,d_{H}(a_{0})\equiv\frac{a_{0}\,\eta_{0}}{a_{0}/\sqrt{-k}}=\sqrt{-k}\,\eta_{0}\,.$
(15)

This ratio can be used to parametrize the different GREA scenarios and compare
with $\Lambda$CDM, see Fig. 1, where

$H_{0}\,d_{H}(a)=\frac{2}{\sqrt{\Omega}_{\rm
M}}\,a^{3/2}\cdot{}_{2}F_{1}\left[\frac{1}{2},\,\frac{1}{6},\,\frac{7}{6},\,\frac{\Omega_{\rm
M}-1}{\Omega_{\rm M}}\,a^{3}\right]\,,$ (16)

is the expression for $d_{H}(a)$ in $\Lambda$CDM.

Solving the Hamiltonian constraint equation (13) with initial conditions set
by the Cosmic Microwave Background, deep in the matter era, where
$a_{i}(\tau)=\tau^{2}\,\Omega_{\rm M}/4$, with $\Omega_{\rm M}=0.31$,
$\Omega_{\rm K}=0$ and $H_{0}=67.8$ km/s/Mpc Aghanim et al. (2020), we find an
effective dark energy contribution $\Omega_{\Lambda}\simeq 0.70$ today, for a
narrow range of values of $\alpha\simeq 1$, see Table 1.

$\alpha$ | $\Omega_{\rm M}$ | $\Omega_{\Lambda}$ | $h_{0}$ | $w_{0}$ | $w_{a}$ | $z_{c}$
---|---|---|---|---|---|---
1.139 |  0.368 |  0.632 |  0.786 |  -0.908 |  -0.315 |  0.931
1.119 |  0.336 |  0.664 |  0.749 |  -0.937 |  -0.331 |  0.840
1.099 |  0.307 |  0.693 |  0.714 |  -0.966 |  -0.345 |  0.757
1.081 |  0.279 |  0.721 |  0.682 |  -0.996 |  -0.354 |  0.678
1.062 |  0.254 |  0.746 |  0.651 |  -1.028 |  -0.365 |  0.603
1.043 |  0.231 |  0.769 |  0.622 |  -1.060 |  -0.367 |  0.533
|  0.310 |  0.690 |  0.678 |  -1.000 |  -0.000 |  0.645
Table 1: The GREA values of the DM and DE content, the rate of expansion
today, $H(0)=100\,h_{0}$ km/s/Mpc, the effective EOS parameters
($w_{0},\,w_{a}$) and the redshift $z_{c}$ of the corresponding coasting
point, for different values of the parameter $\alpha$. From top to bottom, the
values of $\alpha$ correspond to the (red, orange, yellow, green, blue,
purple) colors in the rest of the figures. We compare those values with the
corresponding ones for $\Lambda$CDM (bottom row).

It is important to emphasize that GREA has exactly the same number of free
parameters as $\Lambda$CDM, i.e. ($\alpha,\,\Omega_{\rm M},\,H_{0}$) versus
($\Omega_{\Lambda},\,\Omega_{\rm M},\,H_{0}$), and nevertheless has very
different dynamics, both at the present times and in the far future. Moreover,
while in $\Lambda$CDM there is a coincidence (“why now”) problem, i.e. why is
$\Omega_{\Lambda}\simeq\Omega_{\rm M}$ precisely today?, in the case of GREA
the growth of the cosmic horizon during the matter era is responsible for the
transition to a dominance of the entropic repulsive forces associated with the
acceleration of the universe. The moment at which this happens depends on the
actual value of the spatial boundary today in units of the causal horizon, as
can be seen in Fig. 2a, where we show the coasting points for the different
realizations of GREA. In this scenario, the value of $\alpha$ depends on how
far is our worldline from the bubble wall, as pictured in Fig. 1 of
Ref.García-Bellido and Espinosa-Portalés (2021). If the bubble wall would have
been further away, the onset of acceleration would have happened later for the
same value of $\Omega_{\rm M}$ inside our universe. However, in an open
universe the largest volume is near the bubble wall, so it is expected that
most observers will lie far from the center of the inflated bubble and near
its edges. Therefore the finetuning associated with the coincidence problem is
strongly alleviated. Note also that different GREA scenarios populate the
parameter space of the cosmological background in a way that can be easily
distinguished from $\Lambda$CDM, see Fig. 3, with multiple panels, showing
that of every 2D parameter space there is always a value of $\alpha$ that can
accommodate a universe like that of $\Lambda$CDM, except for the value of
$w_{a}$, which is always large and negative, in agreement with recent
observations.

Figure 3: Different GREA scenarios in the plane of $\Omega_{\rm M}$, the
coasting point $z_{c}$, the present rate of expansion $H(0)$ (in units of 100
km/s/Mpc) and the EOS parameters $(w_{0},\,w_{a})$. We also plot the
$\Lambda$CDM values (black dashed lines).

Furthermore, in GREA the value of the rate of expansion today is not fixed at
the CMB, like in the case of $\Lambda$CDM. Here, depending on $\alpha$, the
present rate of expansion, $H(0)$, could be larger or smaller than the
$\Lambda$CDM value derived from the CMB, see Fig. 4a and Table I. In fact, for
certain values of the spatial curvature one can resolve in a natural way the
so-called Hubble tension Di Valentino et al. (2021).

Finally, while in $\Lambda$CDM, or its variant $w_{0}w_{a}\Lambda$CDM, the
value of the equation of state of Dark Energy is given by
$w(a)=w_{0}+w_{a}(1-a)$, as an expansion around the present value, in GREA we
have the whole function $w(a)$ to compare with observations, as can be seen in
Fig. 2b. The evolution of the scale factor in GREA allows one to compute the
effective equation of state of dark energy,

$\displaystyle 1+w(a)$ $\displaystyle=$ $\displaystyle\frac{-d\ln}{d\ln
a^{3}}\left[\frac{H^{2}(a)}{H_{0}^{2}}-\Omega_{\rm
M}\\!\left(\frac{a_{0}}{a}\right)^{3}\right]$ (17) $\displaystyle=$
$\displaystyle\frac{-d\ln}{d\ln
a^{3}}\left[\frac{\sinh(2\tau)}{a^{2}(\tau)}\right]\,.$

In the matter dominated era, $\tau\propto\sqrt{a}$, and therefore $1+w(a)\to
1/2$ as $a\to a_{i}$, see Fig. 2b. In the far future, a combination of
curvature and decaying entropy production gives a dependence $1+w(a)\to 2/3$
as $a\to\infty$. The rate of expansion therefore will decay as curvature in
the far future, diluting the universe, and ending in flat Minkowski space-
time, see Fig.4b. Note that this is very different from the standard
$\Lambda$CDM scenario with the asymptote to empty de Sitter.

Figure 4: The upper panel shows the ratio of the evolution of the rate of
expansion in GREA (color lines) versus $\Lambda$CDM (black dot-dashed line).
The lower panel shows the actual evolution in GREA and $\Lambda$CDM, with the
past (matter) and future (curvature) asymptotics (dotted blue lines).

In fact, in GREA theory there is no freedom to chose ($\Omega_{\rm
M},\,H(0),\,w_{0},\,w_{a}$) independently, as it happens in
$w_{0}w_{a}\Lambda$CDM. The actual values of $(w_{0},\,w_{a})$ are precisely
given by the theory, just like the present rate of expansion,
$H(0)=100\,h_{0}$ km/s/Mpc, once you fix the ratio $\alpha$, see Fig. 3 for
some examples. As a consequence, if the cosmological data favours values that
are far away from $\Lambda$CDM $(w_{0}=-1,\,w_{a}=0)$, then we can test the
whole theory by comparing those values with the corresponding predictions for
both $H(0)$ and ($\Omega_{\rm M},\,\Omega_{\Lambda}$), as can be seen in Table
1.

We note in Fig. 5 that the GREA scenario is not in tension with the recent
DESI-BAO Adame et al. (2024) \+ DES-SNIa Abbott et al. (2024) data, although a
proper analysis, including the background and linear perturbations of GREA,
has to be performed with the full data set (CMB, SNIa, LSS, BAO, RSD, ISW,
$H_{0}$, etc.). In fact, the amount of information that we can extract from
present observations is larger than what can be inferred from their projection
on to the plane ($w_{0},\,w_{a}$), and future data will be even more
constraining at different redshifts.

Figure 5: The plot shows the parameter space $(w_{0},\,w_{a})$. Cosmological
data seemed to prefer GREA from $\Lambda$CDM at the 2$\sigma$ level, already
in 2021, see Ref. Arjona et al. (2022). Nowadays, the SN-Ia from DES Y5, and
the BAO from DESI Y1, give compelling evidence that $\Lambda$CDM is excluded
at close to 4$\sigma$, while GREA is within the 2$\sigma$ contours. Figure
adapted from Ref. Adame et al. (2024).

##  IV IV. Linear perturbation theory in GREA

So far the description of GREA has been done at the level of the background
evolution, while we have ignored the small matter fluctuations that give rise
to galaxies and the large scale structure (LSS) of the universe. In order to
understand the effect of GREA on the growth of matter perturbation we need to
develop a linear perturbation theory in the accelerating background of GREA,
which will drive a rate of growth that may be distinguishable in principle
from that of $\Lambda$CDM. Furthermore, in order to compare with observations,
this linear growth has to be properly characterized as a function of the free
parameters of the theory.

The equations driving the evolution of matter fluctuations can be derived from
the covariant conservation of the total energy-momentum tensor (3). The
background metric in conformal time is that of a curved FLRW with a Newtonian
potential $\Phi$ and curvature fluctuation $\Psi$,

$ds^{2}=a^{2}(\eta)\left[-(1+2\Phi)\,d\eta^{2}+(1-2\Psi)\gamma_{ij}dx^{i}dx^{j}\right]\,.$
(18)

In the absence of shear and vorticity (i.e. for a spatially-symmetric matter
energy-momentum tensor), the Newtonian potential and curvature fluctuation
satisfy $\Phi=\Psi$. Then the Bardeen equation Bardeen (1980) for the gauge
invariant curvature fluctuation in flat space ($K=0$) with a cosmological
constant and entropic pressure can be written, in terms of a general equation
of state $w=\sum p_{i}/\sum\rho_{i}$, as

$\Phi^{\prime\prime}(\eta)+3{\cal H}\Phi^{\prime}(\eta)-3w{\cal
H}^{2}\Phi=c_{s}^{2}\,\nabla^{2}\Phi\,,$ (19)

where the Poisson equation

$\nabla^{2}\Phi=4\pi G\rho\,a^{2}\,\delta$ (20)

relates the curvature fluctuation $\Phi$ to the gauge-invariant density
contrast, $\delta$. Expanding Bardeen’s equation (19) for a matter fluid with
negligible speed of sound, $c_{s}^{2}\simeq 0$, we find

$\Phi^{\prime\prime}(\eta)+3{\cal H}\Phi^{\prime}(\eta)+(\Lambda a^{2}-8\pi
G\,a^{2}p_{S})\Phi(\eta)=0\,,$ (21)

with $p_{S}=-\rho_{H}$ is the entropic pressure (13). We can then rewrite the
Bardeen equation, using the Poisson equation, $\Phi\propto\delta/a$, and the
identity

${\cal H}^{\prime}+2{\cal H}^{2}=\frac{3}{2}H_{0}^{2}\frac{\Omega_{\rm
M}}{a}+\Lambda a^{2}-8\pi G\,a^{2}p_{S},$

to finally write

$\delta^{\prime\prime}(\tau)+\frac{a^{\prime}}{a}\,\delta^{\prime}(\tau)-\frac{3}{2}\,\frac{\Omega_{\rm
M}}{a(\tau)}\,\delta(\tau)=0\,,$ (22)

where $\tau=a_{0}H_{0}\eta$, and the time evolution is determined by eq.(13).
We can now integrate with initial condition deep in the matter era
$\delta_{i}(\tau)\propto a_{i}(\tau)\propto\tau^{2}$. Once the solution is
found for the density contrast, we can derive the time evolution of the
Newtonian potential $\Phi$ via the Poisson equation for matter fluctuations
(20).

Figure 6: We show the damped evolution of the density contrast (upper panel)
and the Newtonian potential (lower panel), as a function of the scale factor
for $\Lambda$CDM (black dot-dashed line) and GREA (colors correspond to those
of Fig. 2). The gray line corresponds to the growth in an open matter
dominated cosmology. It is clear that both $\Lambda$CDM and GREA perturbations
deviate significantly from a matter-only scenario after $a\sim 0.2$ or $z\sim
4$.

For flat $\Lambda$CDM, the solution to eq. (22) can be written explicitly in
terms of Hypergeometric functions,

$\displaystyle\delta(a)$ $\displaystyle=$
$\displaystyle\delta_{0}\,a\cdot{}_{2}F_{1}\left[1,\,\frac{1}{3},\,\frac{11}{6},\,\frac{\Omega_{\rm
M}-1}{\Omega_{\rm M}}\,a^{3}\right]\,,$ (23) $\displaystyle\tau(a)$
$\displaystyle=$ $\displaystyle\frac{2\sqrt{a}}{\sqrt{\Omega}_{\rm
M}}\,\cdot{}_{2}F_{1}\left[\frac{1}{2},\,\frac{1}{6},\,\frac{7}{6},\,\frac{\Omega_{\rm
M}-1}{\Omega_{\rm M}}\,a^{3}\right]\,,$ (24) $\displaystyle\Phi(a)$
$\displaystyle=$
$\displaystyle\Phi_{0}\,\cdot{}_{2}F_{1}\left[1,\,\frac{1}{3},\,\frac{11}{6},\,\frac{\Omega_{\rm
M}-1}{\Omega_{\rm M}}\,a^{3}\right]\,.$ (25)

In the case of GREA, we must solve first the background evolution (13), and
then the linear perturbation eq. (22). We compare in Fig. 6 the time evolution
of both the density contrast $\delta(a)$ and the Newtonian potential $\Phi(a)$
for $\Lambda$CDM and GREA for the values of $\alpha$ in Table 1.

Figure 7: We show the angular diameter distance for $\Lambda$CDM (black dot-
dashed line) and GREA (colors correspond to those of Fig. 2).

###  IV.1 The growth function and index

In Large Scale Structure, the redshift-space distortions in the two-point
correlation function $\xi(r)$ depends on a particular combination of the
density contrast that appears in the velocity perturbation known as the growth
function Amendola et al. (2018), $f(a)=d\ln\delta(a)/d\ln a$, where
$\delta(a)$ is the linear matter perturbation. In the case of $\Lambda$CDM the
growth function has a compact expression Bueno Belloso et al. (2011),

$f(a)=\Omega_{\rm M}^{1/2}(a)\frac{P_{-1/6}^{-5/6}\Big{[}\Omega_{\rm
M}^{-1/2}(a)\Big{]}}{P_{1/6}^{-5/6}\Big{[}\Omega_{\rm M}^{-1/2}(a)\Big{]}}\
\equiv\ \Omega_{\rm M}^{\gamma}(a)\,,$ (26)

where $\Omega_{\rm M}(a)=\Omega_{\rm M}/a^{3}H^{2}(a)$ and $P_{n}^{m}(z)$ are
the associated Legendre polynomials. The growth index $\gamma$ is defined as
the power of $\Omega_{\rm M}(a)$ in the growth function,

$\gamma(a)=\frac{1}{2}+\frac{1}{\ln\Omega_{\rm
M}(a)}\ln\left[\frac{P_{-1/6}^{-5/6}\Big{[}\Omega_{\rm
M}^{-1/2}(a)\Big{]}}{P_{1/6}^{-5/6}\Big{[}\Omega_{\rm
M}^{-1/2}(a)\Big{]}}\right]\,.$ (27)

We have shown in Fig. 8 the growth index for $\Lambda$CDM and GREA for
different values of $\Omega_{\rm K}$ as a function of redshift. Future surveys
should be able to distinguish between these two scenarios.

Figure 8: We show the growth function $f(z)\,\sigma_{8}(z)$ and growth index
$\gamma(z)$ for $\Lambda$CDM (black dot-dashed line) and GREA (colors
correspond to those of Fig. 2).

###  IV.2 The BAO scale and the $S_{8}$ tension

One of the most important physical quantities that one can measure with deep
galaxy surveys is the BAO scale, corresponding to the size of the sonic
horizon as a function of redshift. The angular BAO scale depends on both the
sonic horizon at the baryon drag epoch, $\chi_{s}(z_{\rm drag})=110$ Mpc/h,
and the angular diameter distance, while the radial BAO scale depends also on
the rate of expansion,

$\displaystyle\theta_{\rm BAO}(z)$ $\displaystyle=$
$\displaystyle\frac{\chi_{s}(z_{\rm drag})}{d_{A}(z)}\,,$ (28) $\displaystyle
r_{\rm BAO}(z)$ $\displaystyle=$ $\displaystyle\frac{H(z)}{c}\,\chi_{s}(z_{\rm
drag})\,.$ (29)

The angular diameter distance for flat $\Lambda$CDM is given explicitly in
terms of $\tau(a)$ of eq. (24)

$H_{0}\,d_{A}(z)=\frac{\tau(0)-\tau(z)}{1+z}\,.$ (30)

We show in Fig. 7 the angular diameter distance as a function of redshift.

The luminosity distance needed to make connection with the SN-Ia data is
simply related to the angular diameter distance via the Etherington relation,

$d_{L}(z)=(1+z)^{2}\,d_{A}(z)\,.$

Another cosmological observable which is derived from weak lensing is the
magnitude of the fluctuations today on a scale of 8 Mpc/h. For $\Lambda$CDM
and GREA we have different predictions and one can write the ratio as a
function of $\alpha$, as shown in Fig. 9. For those values of $\alpha$ that
resolve the $H_{0}$ tension, predicting today $H(0)=73\pm 1$ km/s/Mpc Di
Valentino et al. (2021), the $\sigma_{8}$ value predicted by GREA is around 4%
lower than that of $\Lambda$CDM, resolving also the $S_{8}$ tension Di
Valentino et al. (2021).

We should be cautious here, since in order to claim a resolution of the so-
called $H_{0}$ and $S_{8}$ tensions within GREA requires a full CMB, weak
lensing and Cepheid-TRGB analysis. In fact, recent claims of resolution of
each of these tensions have been discussed recently and it could be that they
are associated with unknown systematics yet to be understood.

Figure 9: We show the ratio of the amplitude of fluctuations (red) at 8 Mpc/h
and the ISW effect (blue) between GREA and $\Lambda$CDM, as a function of
$\alpha$. There is a value for which both ratios are one, for $\alpha=1.092$.

###  IV.3 The ISW effect

An important observable that can be obtained from the CMB maps is the
Integrated Sachs-Wolfe (ISW) effect, which arises because of the decrease in
the Newtonian potential $\Phi$ due to the presence of an accelerating
expansion, as photons travel since decoupling and feel the redshift induced by
the expansion of the Universe. This corresponds to a temperature decrement
given by

$\frac{\delta
T}{T}=2\int_{0}^{\eta_{0}}d\eta\,\Phi^{\prime}(\eta_{0}-\eta)\,.$ (31)

The equation for $\Phi(\eta)$ can be derived in flat $\Lambda$CDM as
$\Phi^{\prime\prime}(\tau)+3(a^{\prime}/a)\,\Phi^{\prime}(\tau)+\Lambda\,a^{2}\Phi(\tau)=0$,
see eq. (21).

For GREA the Newtonian potential evolves according to the density contrast,
see Fig. 6b. One can integrate along the line of sight from the last
scattering surface, see eq. (31), and compare with the predicted value in
$\Lambda$CDM. We show the ratio in Fig. 9. While the effect of GREA on
$\sigma_{8}$ can be at most a few percent, in the case of the ISW effect, the
correction can reach 20 or 30%.

Note that we have studied here only the homogenous ISW effect, due to the
global expansion of the universe. There is a much more pronounced ISW effect
towards deep voids in the cosmic web, where photons traversing large voids get
an extra decrement in the ISW temperature contrast due to the growth of the
Newtonian potential as the void becomes emptier. Such an effect is known to be
significantly larger and cannot be accounted for by the evolution of the LSS
around large voids in $\Lambda$CDM Kovács et al. (2019). However, GREA can
account for this extra acceleration due to the production of entropy
associated with the formation of the cosmic web and the growth of voids as a
consequence of gravitational collapse of large structures like filaments,
sheets and superclusters García-Bellido (2024). This effect is not computed
here and may be relevant for the comparison with observations in the future.

##  V V. Conclusions

General Relativistic Entropic Acceleration Theory is a covariant framework
that can explain the present acceleration of the universe from the entropic
force associated with the growth of the cosmological horizon, without the need
to introduce a fundamental cosmological constant, i.e. $\Lambda=0$ always. The
cosmic entropic acceleration associated with the homogeneous expansion of the
cosmological causal horizon can be described in terms of just one parameter
$\alpha$, the ratio between the horizon distance and the spatial curvature.

The entropic acceleration due to the causal horizon starts to be important
only recently, thanks to the $\sinh(2\tau)$ dependence in Eq. (13),
alleviating the coincidence problem. It is the fact that the universe is old
and big that GREA started to dominate the expansion of the universe only
recently. In the past, such a term was exponentially negligible, for ratios
$\alpha\simeq 1$. One only has to live near the edge of the bubble, where most
of the volume is. For larger horizon distances the acceleration would have
started later for the same amount of matter.

The recent evidence for deviations from the $\Lambda$CDM paradigm seems to
agree remarkably well with the predictions of GREA at the homogeneous
background level. In particular, it resolves the Hubble tension by shifting
the coasting point $z_{c}$ to higher redshifts and extending the period of
accelerated expansion with respect to $\Lambda$CDM. Nevertheless, this
entropic acceleration doesn’t last for ever and eventually the GREA entropic
term is diluted just like matter, so that the universe ends in a flat
Minkowski space-time, instead of de Sitter.

We find that some background observables in GREA, like $\Omega_{\rm M}$ and
the effective equation of state parameter today $w_{0}$, remain very close to
those of $\Lambda$CDM, for values of $\alpha\simeq 1$, see Fig. 2. In fact, we
find a dark energy parameter today, $\Omega_{\Lambda}$, consistent with
observations, which gives an effective $\Lambda$ that is 120 orders of
magnitude below the Planck scale, without finetuning. Meanwhile, for the same
values of $\alpha$, the EOS slope $w_{a}$ is found to be significantly large
and negative, and closer to the observed values today, without introducing a
new parameter to explain the data.

We have also studied the linear growth of matter perturbations in the context
of the homogeneous expanding background driven by the growth of entropy from
the cosmological horizon. We find that the rate of growth of matter
fluctuations in GREA slows down due to the accelerated expansion and
alleviates the $\sigma_{8}$ tension. We further study the redshift-space
distortions induced by the growth function of matter fluctuations, which
enters the two-point correlation function of galaxies measured in LSS surveys.
We find that the corresponding ISW effect is significantly larger than in
$\Lambda$CDM, but not enough to explain the observed anomalous ISW effect.

The main effect of GREA on the formation of structure is the accelerated
expansion due to the rapid growth of entropy of the cosmological horizon
during the matter era. However, there are also local accelerations due to the
growth of entropy associated with the growth of entropy associated with the
gravitational collapse of large scales structures in the cosmic web, which is
responsible for emptying the large voids and can induce an additional
ananomalous ISW effect. These inhomogeneous local accelerations have not been
discussed here and will be the subject of a future publication.

In summary, it is encouraging that most of the tensions and anomalies of the
standard model of cosmology are resolved, or at least alleviated, by the
inclusion of this transient period of acceleration of the universe based on
known fundamental physics. We leave for a future publication the detailed
comparison of the predictions of GREA theory with the present cosmological
data from deep and wide galaxy surveys like DES, DESI and Euclid.

##  VI Acknowledgements

The author thanks Julien Lesgourgues for constructive criticisms on a
preliminary draft that has helped improve the paper. He also acknowledges
support from the Spanish Research Project PID2021-123012NB-C43 [MICINN-FEDER],
and the Centro de Excelencia Severo Ochoa Program CEX2020-001007-S at IFT.

## References

  * Espinosa-Portalés and García-Bellido (2021) L. Espinosa-Portalés and J. García-Bellido, Phys. Dark Univ. 34, 100893 (2021), eprint 2106.16012.
  * García-Bellido (2024) J. García-Bellido, Phys. Dark Univ. 44, 101491 (2024), eprint 2306.10593.
  * García-Bellido (2023) J. García-Bellido (2023), eprint 2302.08537.
  * Gibbons and Hawking (1977) G. W. Gibbons and S. W. Hawking, Phys. Rev. D 15, 2752 (1977).
  * York (1972) J. W. York, Jr., Phys. Rev. Lett. 28, 1082 (1972).
  * García-Bellido and Espinosa-Portalés (2021) J. García-Bellido and L. Espinosa-Portalés, Phys. Dark Univ. 34, 100892 (2021), eprint 2106.16014.
  * Abbott et al. (2024) T. M. C. Abbott et al. (DES) (2024), eprint 2401.02929.
  * Adame et al. (2024) A. G. Adame et al. (DESI) (2024), eprint 2404.03002.
  * Arjona et al. (2022) R. Arjona, L. Espinosa-Portalés, J. García-Bellido, and S. Nesseris, Phys. Dark Univ. 36, 101029 (2022), eprint 2111.13083.
  * Gagnon and Lesgourgues (2011) J.-S. Gagnon and J. Lesgourgues, JCAP 09, 026 (2011), eprint 1107.1503.
  * Hawking (1975) S. W. Hawking, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * ’t Hooft (1993) G. ’t Hooft, Conf. Proc. C 930308, 284 (1993), eprint gr-qc/9310026.
  * Susskind (1995) L. Susskind, J. Math. Phys. 36, 6377 (1995), eprint hep-th/9409089.
  * Aghanim et al. (2020) N. Aghanim et al. (Planck), Astron. Astrophys. 641, A6 (2020), eprint 1807.06209.
  * Di Valentino et al. (2021) E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, Class. Quant. Grav. 38, 153001 (2021), eprint 2103.01183.
  * Bardeen (1980) J. M. Bardeen, Phys. Rev. D 22, 1882 (1980).
  * Amendola et al. (2018) L. Amendola et al., Living Rev. Rel. 21, 2 (2018), eprint 1606.00180.
  * Bueno Belloso et al. (2011) A. Bueno Belloso, J. García-Bellido, and D. Sapone, JCAP 10, 010 (2011), eprint 1105.4825.
  * Kovács et al. (2019) A. Kovács et al. (DES), Mon. Not. Roy. Astron. Soc. 484, 5267 (2019), eprint 1811.07812.


# Sampling from a Gaussian distribution conditioned on the level set of a
piecewise affine, continuous function

Jesse Windle
<EMAIL_ADDRESS>
###### Abstract

We consider how to use Hamiltonian Monte Carlo to sample from a distribution
whose log-density is piecewise quadratic, conditioned on the sample lying on
the level set of a piecewise affine, continuous function.

##  1 Introduction

Our motivation comes from modeling monocot (e.g. maize) root growth in a
piecewise linear fashion and inferring the path of a root given a point along
its trajectory. We defer details of this application to another paper. Our
focus here is a question that arose when studying that problem.

Suppose that a priori $X\sim N(\mu,I_{n})$. Subsequently, we learn that
$\ell(X)=0$ where $\ell$ is a piecewise affine, continuous function. How does
one sample $(X|\ell(X)=0)$?

The answer to this question extends to the more general problem, which can be
roughly phrased as: how does one sample $(X|\ell(X)=0)$ when, marginally, $X$
comes from a density proportional to

$\exp[-\frac{1}{2}x^{\prime}M_{i}x+r_{i}^{\prime}x+k_{i}\Big{]},\text{ if
}x\in R_{i},\;i=1,\ldots K$

where $R_{i}\subseteq\mathbb{R}^{n}$ and $R=\cup_{i=1}^{K}R_{i}$ is a
connected set. Notice that, in effect, this provides a way to sample any
distribution, approximately, on discretized volumes, conditioned upon an
observation lying on surface embedded in $\mathbb{R}^{n}$, which has been
discretized into linear pieces.

We will use Hamiltonian Monte Carlo (HMC) to do this. To be clear: our
application of interest is the discretized version, not the approximation of
something continuous, since we can use HMC to simulate from a smooth density.
In that case, one must discretize an ODE. We have, in effect, made a trade for
what must be discretized, sampling exactly from the above discretization by
being able to solve the dynamics of motion exactly.

Part of the aim of this paper is to encourage a physics-based intuition of
particle dynamics for HMC, specifically as it relates to how one deals with
truncations and step changes in densities. To that end, we have synthesized
the work of Pakman and Paninski (2014), who deal with truncations, and Mohasel
Afshar and Domke (2015), who deal with steps, and then place that within the
context of a particle moving on a manifold within a higher dimensional space.

Accompanying this paper is a Python package and example Jupyter notebooks
(Windle (2023)) to show how to use the software. In many modeling scenarios,
the excellent HMC software Stan (Stan Development Team (2023a)) will work,
even with truncations or for some non-smooth densities (see e.g. Stan
Development Team (2023b), section 3.7), and we encourage you to use that
software. However, to the best of our knowledge and as of this writing, Stan
may be unable to accomodate some use cases, like arbitrary truncations by
hyperplane in multiple dimensions or when the piecewise affine, continuous
function used for conditioning is complex.

##  2 Background

We will be using HMC to generate samples from the distribution above.
Hamiltonian Monte Carlo is inspired by Hamiltonian mechanics. (Kibble and
Berkshire (2004) is a well-known introductory text to classical mechanics. I
found the HMC notes by Vishnoi (2021) to be a useful, quick introduction and
borrow from them here. Betancourt (2017) is a more extensive introduction.
Neal (2011) is another useful review for some of the key concepts.)

In physical systems that are governed by a time-independent potential, the
total energy can be written as

$E=V(r)+\frac{1}{2m}\|m\dot{r}\|^{2},$

which is constant. It is often the case that the positions and momenta of the
system can be written in only a few coordinates, called generalized positions,
$x$ and momenta, $s$ so that $E=H(x,s)$. One can show using variational
methods, that the evolution of a particle in $xs$-phase space is governed by

$\displaystyle\dot{x}$ $\displaystyle=\nabla_{s}H$ $\displaystyle\dot{s}$
$\displaystyle=-\nabla_{x}H.$

Let $\varphi_{\tau}$ be the map
$(x(\tau),s(\tau))=\varphi_{\tau}(x_{0},s_{0})$ that represents the state of
the system after time $\tau$ and starting at $(x_{0},s_{0})$. If one defines a
distribution in phase space so that its density is proportional to
$\exp\\{-H(x,s)\\}$, and if one samples random variables $(X,S)$ from this
distribution, then the map of the evolution $\varphi_{\tau}(X,S)$ will
preserve the distribution, where $\tau$ is any time step. In other words, the
distribution is stationary under the mapping $\varphi_{\tau}$.

For the purposes of random variate generation, one can consider
$H(x,s)=V(x)+T(s)$ where $V(x)$ is the negative logarithm of a distribution of
interest and $T(s)$ is the negative logarithm of a proposal distribution. We
can construct a Markov chain as follows. Let $X_{0}$ be an arbitrary starting
position. Let $X_{n}=\varphi_{\tau}(X_{n-1},S_{n-1})$ where $S_{n}$ is a
random variable sampled from the distribution derived from $\exp(-T(s))$. When
$T$ and $V$ satisfy regularity conditions and $\tau$ is an acceptable step
size, $(X_{n})_{n=1}^{\infty}$ is a Markov chain that converges to the
distribution of interest.

###  2.1 Example: Gaussian distribution

Let us consider sampling from a Gaussian random variable using HMC. This may
seem trivial, since we already known how to sample from a Gaussian
distribution without HMC, but it will be instructive as well as useful for
later. Suppose

$\displaystyle V(x)=\frac{1}{2}x^{\prime}Mx-r^{\prime}x$ $\displaystyle
T(s)=\frac{1}{2}s^{\prime}M^{-1}s.$

The Hamiltonian for the system is then

$H(x,s)=\frac{1}{2}x^{\prime}Mx-r^{\prime}x+\frac{1}{2}s^{\prime}M^{-1}s$

and the evolution in this phase space is given by

$\displaystyle\dot{x}$ $\displaystyle=\frac{dH}{ds}=M^{-1}s$
$\displaystyle\dot{s}$ $\displaystyle=-\frac{dH}{dx}=-Mx+r.$

Note that $s$ in this case _is not_ the actual momenta, but rather a linear
combination of the particle’s momentum. Together this implies that the motion
of $x$ is governed by

$\ddot{x}+x=\mu$

where $\mu=M^{-1}r$, which has the solution

$x_{i}(t)=\mu_{i}+a_{i}\sin(t)+b_{i}\cos(t)$ (1)

with $b_{i}=x_{i}(0)-\mu_{i}$ and
$a_{i}=\dot{x}_{i}(0)=-(M^{-1})_{i\cdot}s_{0}$. We can sample $s$ from the
correct marginal distribution, $s\sim N(0,M)$, to generate the Markov chain,
or by the transformation $\dot{x}=-M^{-1}$, we see that we can also sample
$\dot{x}(0)$ directly as $N(\mu,M^{-1})$ — yes, we are saying we can sample
$N(0,M^{-1})$ to generate a Markov chain whose stationary distribution is
$N(\mu,M^{-1})$. If we chose $\tau=\pi/2$, we can even generate independent
and identical samples. This may seem rather tautological at the moment, but it
will be very useful when dealing with truncated Gaussian random variables or
piecewise quadratic likelihoods.

###  2.2 Example: piecewise (or truncated) Gaussian

To see why where HMC becomes useful, let us now consider a potential defined
piecewise, so that

$V(x)=\begin{cases}\frac{1}{2}x^{\prime}M_{1}x-r_{1}^{\prime}x+k_{1},\;\text{
if }f^{\prime}x+g>0,\text{ region 1}\\\
\frac{1}{2}x^{\prime}M_{2}x-r_{2}^{\prime}x+k_{2},\;\text{ if
}f^{\prime}x+g<0,\text{ region 2}\end{cases}.$ (2)

Previously, the constant term in the polynomial did not matter, but now it
does since that will weight each piece differently. Also, we could let e.g.
$k_{2}=\infty$ to effectively truncate the distribution.

We can now evolve a particle by Hamiltonian dynamics, but we have to be
mindful of the boundary. The key question is: when does the particle run into
the boundary? Following Pakman and Paninski (2014), let

$\displaystyle K(t)$ $\displaystyle=f^{\prime}x(t)+g$
$\displaystyle=f^{\prime}a\sin(t)+f^{\prime}b\cos(t)+(f^{\prime}\mu+g)$
$\displaystyle=u\cos(t+\phi)+h$

where

$u=\sqrt{(f^{\prime}a)^{2}+(f^{\prime}b)^{2}},$
$\tan\phi=\frac{f^{\prime}a}{f^{\prime}b}$

and

$h=f^{\prime}\mu+g.$

The particle hits a boundary at the first time that $K(t)=0$, which can only
happen if $u>|h|$. (As a technical matter, when computing $\phi$, one should
use the atan2 function, which returns the angle in $[-\pi,\pi]$ so that
$f^{\prime}a$ and $f^{\prime}b$ have the correct signs.)

We want to find the time $t>0$ in which a particle hits a boundary, as it is
going in a direction that would exit the region. To that end, if
$u_{j}<|h_{j}|$, let $\tau=\infty$, otherwise let

$\tau=\min\Big{\\{}t^{(n)}=\arccos(-h/u)-\phi+\pi
n:n\in\mathbb{Z},\;t^{(n)}>0,\;K^{\prime}(\tau^{(n)})<0\Big{\\}}.$

To clarify: first, we have to find positive $\tau$ hence the adjustment by
multiples of $\pi$. Second, we want to make sure that the particle is exiting
the boundary, hence the condition on the derivative. If $\tau\geq\tau_{max}$,
$\tau_{max}$ the max step size, then we take a step $\tau_{max}$. If
$\tau<\tau_{max}$, then we need to adjust the momentum of the particle as it
encounters the boundary and then let it continue on its way.

Suppose the particle starts in region 1 and at time $\tau<\tau_{max}$ it
encounters the boundary. The intuition of how to approach the problem is as
follows. Suppose $V(x)$ is piecewise constant so that only $k_{1}$ and $k_{2}$
are non-zero. We can think of this potential as a plain (region 1) and a
plateau (region 2), both of which are totally flat with a constant force of
gravity perpendicular to both surfaces. We can relax the problem to make the
vertical face of the plateau slightly tilted, and even smooth. The key is that
the shape of the relaxation only changes in the direction orthogonal to the
boundary. Under that assumption, the value that clearly matters is the
momentum in this direction. If there is enough momentum to overcome the
potential energy due to gravity, a particle will ascend the plateau. If it
does not, then it will return to region 1\. As we make the relaxation steeper
and steeper, the time in spends on this portion of the curve decreases, so
that in the limit, it will either ascend the plateau and lose momentum in the
direction orthogonal to the boundary, or it will reflect.

The direction orthogonal to the boundary, pointing in, is $u=f/\|f\|$, hence
the velocity (and momentum since $m=1$) is
$v_{\perp}=u^{\prime}\dot{x}(\tau)$. The associated kinetic energy is
$E_{\perp}=\frac{1}{2}v_{\perp}^{2}$. Let

$\Delta V=\lim_{z\rightarrow x(\tau),f^{\prime}z+g<0}V(z)-\lim_{z\rightarrow
x(\tau),f^{\prime}z+g>0}V(z).$

If $E_{\perp}>\Delta V$, then the particle will move into the other region
with a subsequent reduction in momentum in the $u$ direction of

$v_{new}=\sqrt{2(E_{\perp}-\Delta V)},$

which is the key point of Mohasel Afshar and Domke (2015), so that

$\dot{x}(\tau^{+})=(I-uu^{\prime})\dot{x}(\tau^{-})+v_{new}(-u),$

(the minus in front of $u$ ensures we are pointing into the next region),
otherwise all of the momentum in question will be reflected;

$\dot{x}(\tau^{+})=(I-uu^{\prime})\dot{x}(\tau^{-})-v_{\perp}u,$

which we can summarize as

$\dot{x}(\tau^{+})=\dot{x}(\tau^{-})+u\begin{cases}-2v_{\perp},\;v_{\perp}^{2}<2\Delta
V\\\ -v_{\perp}-\sqrt{v_{\perp}^{2}-2\Delta V},v_{\perp}^{2}\geq 2\Delta
V.\end{cases}$

Note, when $k_{2}=\infty$, the boundary is effectively a truncation of the
distribution, and you are guaranteed to reflect.

Given this new velocity, we can re-initialize the dynamics, but now in the
appropriate region with initial conditions $(x(\tau),\dot{x}(\tau^{+}))$ and
max running time $\tau_{max}-\tau$.

##  3 Parameterization of distribution

First, let us define the distribution and its parameters. In the previous
section, we considered a truncated Gaussian distribution. Here, we consider a
Gaussian, subject to a specific type of piecewise linear contraints. We can
think of this as starting with the marginal distribution

$\displaystyle X\sim\exp(-V(x))$ $\displaystyle
V(x)=\begin{cases}\frac{1}{2}x^{\prime}M_{j}x-{r_{j}}^{\prime}x+k_{j},\;x\in
R_{j}\\\ R_{j}=\\{x:{F_{j}}x+g_{j}\geq 0\\}\\\ j=1,\ldots,J\end{cases},$

and then wanting to sample $(X|\ell(X)=0)$ where $\ell$ is a continuous
function that can be written as

$\ell(x)={A_{j}}^{\prime}x+y_{j},\text{ if }x\in R_{j}.$

That provides a general sense of what we want to do, but we need to restrict
the problem further in order to avoid degenerate cases and provide a succinct
parameterization. To begin, let

${f_{i}}^{\prime}x+g_{i}\geq 0,\;i=1,\ldots,n$

be a collection of hyperplanes.

Let the n-tuple $s_{j}\in\\{-1,0,1\\}^{n}$ define the $j$th region, $R_{j}$ by
imposing the constraints

$s_{ji}\Big{[}{f_{i}}^{\prime}x+g_{i}\Big{]}\geq 0,\;i=1,\ldots,n.$

When $s_{ji}=0$, then the constraint is not active; if it is non-zero it is
active and its sign determines which half of the half-space to use. We also
want to keep track of the intended transitions. To that end, let the magnitude
of $s_{ji}$ indicate the index of the adjoining region and when $|s_{ji}|=j$,
we will assume it represents a hard boundary. Let $L$ is the $J\times m$
matrix of these $m$-tuples so that $s_{ji}=L_{ji}$ ($L$ for lookup). If a
particle can move from $j$ to $j^{*}$ and visa versa, then there should be
some constraint $i$ such that $j^{*}=|L_{ji}|\neq 0$ and $j=|L_{j^{*}i}|\neq
0$.

Thus, to recap, the parameterization includes:

  * •

The potential $V$: that is $(M_{j},r_{j},k_{j})$, where $M_{j}$ is a $d\times
d$ matrix, $r_{j}$ is an $\mathbb{R}^{d}$ vector, and $k_{j}$ a scalar, for
$j=1,\ldots,J$.

  * •

$f_{i},i=1,\ldots,m$, vectors in $\mathbb{R}^{n}$ used to construct the
hyperplane boundaries, or $F$ a $m\times n$ dimensional matrix whose rows are
the $f_{i}$.

  * •

$g_{i},i=1,\ldots,m$, the scalars that control the offsets of the hyperplanes,
or $g$ an $m$-dimensional vector.

  * •

$s_{j},j=1,\dots J$, the m-tuples defining the regions, or $L$ a $J\times m$
dimensional matrix whose rows are the $m$-tuples.

  * •

The continuous, piecewise affine function $\ell$: that is $A_{j}$ and $y_{j}$,
where $A_{j}$ is a $n\times d$ and $y_{j}$ is a $\mathbb{R}^{d}$ vector,
$j=1,\ldots,J$.

###  3.1 Restrictions on parameters

Since we will be using HMC, as a technical matter, it will be helpful to avoid
defining regions so that a face of one region abuts the faces of two adjacent
regions. This will prevent us from having to search through all the regions to
determine which region the particle is moving into. We can avoid this by
imposing: for each $j$, for each $i$, there is at most one $j^{\prime}$ such
that $0\neq s_{ij}=-s_{j^{\prime}i}$. We also want $\ell$ to be continuous.
The user should define the regions so that the interior of $R$ is connected.

##  4 Dynamics

Let us return to the potential (2):

$V(x)=\begin{cases}V_{1}(x):=\frac{1}{2}x^{\prime}M_{1}x-r_{1}^{\prime}x+k_{1},\;\text{
if }f^{\prime}x+g>0\;(R_{1})\\\
V_{2}(x):=\frac{1}{2}x^{\prime}M_{2}x-r_{2}^{\prime}x+k_{2},\;\text{ if
}f^{\prime}x+g<0\;(R_{2})\end{cases}.$

and consider the negative log density of the marginal distribution of
interest. Our aim is to sample from conditional upon $\ell(x)=0$ where

$\ell(x)=\begin{cases}A_{1}^{\prime}x+y_{1},\;f^{\prime}x+g>0\;(R_{1})\\\
A_{2}^{\prime}x+y_{2},\;f^{\prime}x+g<0\;(R_{2})\end{cases}$

and $\ell$ is continuous. There are two issues to address: (i) what equations
should drive the dynamics of the particle, and (ii) what happens when a
particle hits the boundary?

###  4.1 Dynamics within a region

Suppose a particle’s dynamics are governed by the unconstrained force that has
the form

$F(x)=-Mx+r,$

which corresponds to a Gaussian potential with precision $M$ and mean
$\mu=M^{-1}r$, but that we then constrain the motion of a particle to the
plane $A^{\prime}x=y$. How do we adjust the force and subsequent potential to
accurately represent the particle’s dynamics, and how can we connect it to the
unconstrained potential?

Let $Q_{1}$ be the unitary column space of $A$ and $Q_{2}$ be the unitary null
space. We need to make two observations. First, the force within the
constrained space is

$-Q_{2}Q_{2}^{\prime}[Mx-r].$

Second, if we decompose $x$ as $Q_{1}Q_{1}^{\prime}x+Q_{2}Q_{2}^{\prime}x$,
then the portion $Q_{1}Q_{1}^{\prime}x$ is constant. Specifically, letting
$A=Q_{1}R$, we see that $Q_{1}^{\prime}x={R^{\prime}}^{-1}y=:z_{1}^{*}$, so
that $Q_{1}Q_{1}^{\prime}x=Q_{1}z_{1}^{*}$.

Hence the force of interest is

$\displaystyle F(x)$
$\displaystyle=-Q_{2}Q_{2}^{\prime}MQ_{2}Q_{2}^{\prime}x-Q_{2}Q_{2}^{\prime}MQ_{1}Q_{1}^{\prime}x+Q_{2}Q_{2}^{\prime}r$
$\displaystyle=-Q_{2}(Q_{2}^{\prime}MQ_{2})Q_{2}x+Q_{2}Q_{2}^{\prime}[r-MQ_{1}z_{1}^{*}].$

Letting $\tilde{r}=[r-MQ_{1}z_{1}^{*}]=M[\mu-Q_{1}z_{1}^{*}]$, the prior
expression implies the potential

$V(x)=\frac{1}{2}x^{\prime}Q_{2}(Q_{2}^{\prime}MQ_{2})Q_{2}x-(Q_{2}Q_{2}^{\prime}\tilde{r})^{\prime}x.$

To write down the Hamiltonian for the purposes of generating random variates,
we need to chose a kinetic energy term that makes generating the dynamics
easy.

From section 2, we know that $s^{\prime}M^{-1}s$ is a good choice, where $M$
is the potential. This suggests we should choose a pseudo inverse of
$Q_{2}(Q_{2}^{\prime}MQ_{2})^{-1}Q_{2}^{\prime}$ when constructing the kinetic
energy term, which we also know is a good choice by solving the problem in
lower dimensional coordinates. To that end, we want

$H(x,s)=\frac{1}{2}x^{\prime}(Q_{2}(Q_{2}^{\prime}MQ_{2})Q_{2}^{\prime})x-(Q_{2}Q_{2}^{\prime}\tilde{r})^{\prime}x+\frac{1}{2}s^{\prime}Q_{2}(Q_{2}^{\prime}MQ_{2})^{-1}Q_{2}s,$

which leads to the dynamics

$\ddot{x}+Q_{2}Q_{2}^{\prime}x=Q_{2}(Q_{2}^{\prime}MQ_{2})^{-1}Q_{2}^{\prime}\tilde{r}.$

This has homogeneous and particular solutions

$x_{h}(t)=a\sin(t)+b\cos(t)$

and

$x_{p}(t)=Q_{2}(Q_{2}^{\prime}MQ_{2})^{-1}Q_{2}^{\prime}\tilde{r}+Q_{1}z_{1}^{*}.$

Then the solution is $x=x_{h}+x_{p}$ where

$b=x(0)-x_{p}$

and

$a=\dot{x}(0),$

with the requirement that $Q_{1}^{\prime}\dot{x}(0)=0$. To generate a sample
from the marginal of $s$, we can sample $s\sim
N(0,Q_{2}(Q_{2}^{\prime}MQ_{2})Q_{2}^{\prime})$, which is like endowing
$\dot{x}(0)$ with initial velocity from $\dot{x}(0)\sim
N(0,Q_{2}(Q_{2}^{\prime}MQ_{2})^{-1}Q_{2}^{\prime})$.

###  4.2 Change in dynamics between regions

To begin, let us consider a particle that is traveling free of any force, but
whose movement is restricted to $R_{1}$ and $R_{2}$. Understanding the
velocity of the particle, and how we can decompose it is key. Keep in mind
that, if the particle is moving from region 1 to region 2, then the particle’s
velocity must be in the null space of $A_{1}^{\prime}$, and when it crosses
the boundary, it will be in the null space of $A_{2}^{\prime}$. It stands to
reason that the components of the velocity that are shared between each region
will not change.

Let us assume we can find a semi-orthonormal matrix $U_{0}$ and unit vectors
$u_{i},i=1,2$ so that $u_{i}$ is orthogonal to $U_{0}$ and the span of $u_{i}$
and $U_{0}$ is the null space of $A_{i}^{\prime}$, $i=1,2$. Further, assume
there is another semi-orthonormal matrix $U_{c}$ such that
$U_{c}^{\prime}U_{0}=0$ and $u_{i}^{\prime}U_{c}=0$, $i=1,2$. The
justification is given in Proposition 2 and Corollary 3 of Appendix A. To
orient $u_{1}$ and $u_{2}$, assume they are pointing into each region so that
$f^{\prime}u_{1}>0$ and $f^{\prime}u_{2}<0$.

How does a particle move from $R_{1}$ into $R_{2}$? Our approach is to set up
a hard boundary so that a particle traveling in $R_{1}$ has an inelastic
collision that bounces it directly into $R_{2}$. As before let
$\dot{x}(\tau^{-})$ be the particle’s velocity when it encounters the
boundary. There will be an instantaneous change in velocity to
$\dot{x}(\tau^{+})$ where

$x(\tau^{+})=(I-u_{\perp}u_{\perp}^{\prime})x(\tau^{-})-u_{\perp}u_{\perp}^{\prime}x(\tau^{-}).$
(3)

where

$u_{\perp}=\frac{u_{2}+u_{1}}{\|u_{2}+u_{1}\|}.$

It will be helpful to define

$u_{\|}=\frac{u_{2}-u_{1}}{\|u_{2}-u_{1}\|}.$

as well. Geometrically speaking, $u_{\perp}$ is the vector that bisects
$u_{2}$ and $u_{1}$ and $u_{\|}$ is the vector the bisects $u_{2}$ and
$-u_{1}$. Because $u_{1}$ and $u_{2}$ are orthogonal to $U_{0}$ and $U_{c}$,
we have that the concatenation of $U=[u_{\perp}\;u_{\|}\;U_{0}\;U_{c}]$ (or
$U=[u_{\perp}\;U_{0}\;U_{c}]$ when $u_{1}=-u_{2}$) is an orthonormal matrix
and $U^{\prime}U=I$, which tells us we can also decompose (3) as

$\dot{x}(\tau^{+})=-u_{\perp}u_{\perp}^{\prime}\dot{x}(\tau^{-})+u_{\|}u_{\|}^{\prime}\dot{x}(\tau^{-})+U_{0}U_{0}^{\prime}\dot{x}(\tau^{-})+\underbrace{U_{c}U_{c}^{\prime}\dot{x}(\tau^{-})}_{=0}.$

The last term is eliminated because $\dot{x}(\tau^{-})$ is in the null space
of $A_{1}^{\prime}$.

To see that the particle does indeed take its trajectory in $u_{2}$, we can
decompose $\dot{x}(t^{-})$ as $-\alpha u_{1}+w$ where $w\in U_{0}$ and
$\alpha>0$, since the particle is leaving the region. The subsequent
trajectory is

$\displaystyle(u_{2}+u_{1})\frac{(u_{2}+u_{1})^{\prime}}{\|u_{2}+u_{1}\|^{2}}\alpha
u_{1}-(u_{2}-u_{1})\frac{(u_{2}-u_{1})^{\prime}}{\|u_{2}-u_{1}\|^{2}}\alpha
u_{1}+w$
$\displaystyle=\alpha(u_{2}+u_{1})\frac{(u_{2}^{\prime}u_{1}+1)}{2(1+u_{2}^{\prime}u_{1})}-\alpha(u_{2}-u_{1})\frac{u_{2}^{\prime}u_{1}-1}{2(1-u_{2}^{\prime}u_{1})}+w$
$\displaystyle=\alpha u_{2}+w.$

In other words, we can project the velocity of the particle into $-u_{1}$ (the
direction of exit) and then transfer all of that velocity to $u_{2}$ and keep
the remaining components of velocity the same. Since the particle is traveling
in $R_{1}$ to begin, with, this is equivalent to removing the velocity due to
$u_{1}$ and then redirecting it to $u_{2}$. That is

$\displaystyle\dot{x}(\tau^{+})$
$\displaystyle=(I-u_{1}u_{1}^{\prime})\dot{x}(\tau^{-})+u_{2}(-u_{1})^{\prime}\dot{x}(\tau^{-}).$

When the potential is continuous at $x(\tau)$, then the total energy will not
change, since momentum is conserved. In order to tackle potentials that are
defined piecewise, we must consider how the momentum of the particle changes.

Instead of considering a flat plain and plateau, consider a plain that is
tilted in the direction of $u_{1}$ and a plateau that is tilted in the
direction of $u_{2}$. The sheer face of the plateau is predominantly in the
direction of $u_{\perp}$, so that it is nearly vertical. Now imagine that we
relax the sheer face so that it is steep and that there is a force of gravity
in the direction of $-u_{\perp}$ and only in the portion of the relaxation. In
order for a particle to ascend to the plateau, we must have the momentum in
$u_{1}$, $v_{1}=u_{1}^{\prime}\dot{x}(\tau^{-})$, and its corresponding energy
exceed the potential energy we must climb:

$E_{u_{1}}:=\frac{1}{2}(\dot{x}(\tau^{-})^{\prime}u_{1})^{2}>\Delta
V=:V_{2}(x(\tau^{-}))-V_{1}(x(\tau^{-})).$

The drop in energy is then

$E_{new}=E_{u_{1}}-\Delta V$

and the subsequent momentum in the $u_{2}$ direction is

$\sqrt{2E_{new}}u_{2}.$

In terms of $v_{1}$, we have

$\sqrt{2E_{new}}=|v_{1}|\sqrt{1-\frac{2\Delta
V}{v_{1}^{2}}}=-v_{1}\sqrt{1-\frac{2\Delta V}{v_{1}^{2}}},$

which clarifies that we lose a fractional portion of the speed in the $u_{1}$
direction as it relates to the change in potential energy. If there is not
enough energy to ascend the potential, then the momentum in the direction of
$u_{1}$ is effectively reversed. That is

$x(\tau^{+})=(I-u_{1}u_{1}^{\prime})\dot{x}(\tau^{-})+\begin{cases}u_{2}\sqrt{2E_{new}},\;\text{
if }E_{new}\geq 0\\\ -u_{1}(u_{1}^{\prime}\dot{x}(\tau^{-})),\;\text{ if
}E_{new}<0\end{cases}.$

######  Aside 1.

Given the above equation, there is a useful computational trick to note.
Assume that the particle is traveling in region 1. We can make an artificial
boundary — that has no effect — by assuming that $\Delta V=0$ and
$u_{2}=-u_{1}$ where $u_{1}$ is a vector in the null space of
$A_{1}^{\prime}$. Similarly, we can force there to be a hard boundary by
taking $\Delta V=0$ and $u_{2}=u_{1}$ where $u_{1}$ is a vector in the null
space of $A_{1}^{\prime}$. These interpretations are useful when vectorizing
the dynamics for many particles, since algorithm above can then be used for
all potential scenarios: no boundary, a hard boundary, or a transition
boundary.

##  5 Examples

###  5.1 One-norm

Here we consider sampling a Gaussian, given that the sample must live on a
level set of the one-norm. Specifically, we want to sample

$X\sim\Big{(}N(0,\Sigma)\;|\>\|X\|_{1}=1\Big{)}.$

Figure 1 shows these samples for $\Sigma=I_{3}$ and $\Sigma$ diagonal with
diagonal elements $(0.1,10.,10.)$.

Figure 1: Sampling $X\sim N(0,\Sigma)$ given $\|X\|=1$. The left image is just
the constraint. The middle image shows a sample for $\Sigma=I_{3}$. The right
image shows a sample for $\Sigma$ diagonal with diagonal elements
$(10,0.1,0.1)$.

###  5.2 Polygonal top

Extending the one norm example, here we have constructed a regular polygon in
the xy-plane and then connected its vertices to the z-axis above and below the
xy-plane to create a top. Figure 2 shows the distribution of $N(0,D)$
restricted to a six-sided top for various covariance structures $D$. In this
case

Figure 2: Sampling $X\sim N(0,D)$ given that $X$ lies on a 6-sided top for
various covariance structures. The colors represent different regions in our
partition of $\ell$.

###  5.3 Positive part

Here we consider sampling on the subspace that arises when using the
$(\cdot)^{+}$ function. Our motivation comes from modeling root growth in a
piecewise linear fashion. We consider a root that has $K=3$ kinks in it. The
depth of the root after traveling a distance of $r>0$ in the x-direction is

$\ell(\Delta x)=\sum_{i=1}^{K+1}\Delta m_{i}(r-\sum_{j=1}^{i}\Delta
x_{i})^{+}$

where $\Delta m_{i},\i=1,\ldots,K+1$ is taken as a parameter here that
determines how the slope of the root changes. The variable we will be sampling
is $\Delta x_{i}\geq 0,i=1,\ldots,3$. The constraint $\ell(\Delta x)=d$ can be
described by the three regions shown in Figure 3.

Figure 3: The regions used to define $\ell(\Delta x)=d$.

Figure 4 shows samples of

$\Delta X\sim\Big{(}N(\mathbb{1}_{3},I_{3})\;\Big{|}\;\ell(\Delta
X)=d\Big{)}.$

Figure 4: Samples of $\Delta X\sim(N(\mathbb{1}_{3},I_{3})\;|\;\ell(\Delta
X)=d)$

##  6 Conclusion

Here we have shown how to start with a distribution that is piecewise
quadratic in the logarithm and sample from that distribution, conditional upon
the sample lying on the the level set of a piecewise affine, continuous
function using Hamiltonian Monte Carlo.

This relies heavily on knowing the exact dynamics of a particle under a
quadratic potential, which can then easily accomodate truncations and step
changes in densities. To that end, we have synthesized the work of Pakman and
Paninski (2014), who deal with truncations, and Mohasel Afshar and Domke
(2015), who deal with steps, and then placed that within the context of a
particle moving on a manifold within a higher dimensional space.

We have written a Python package implementing these methods (Windle (2023)).
However, for many use cases, slightly relaxing the requirement that the sample
come from a lower dimensional subspace enables the use of the HMC software
Stan, even with univariate truncations and for some non-smooth densities. We
encourage that route, if possible. If one needs to consider complex
truncations or piecewise affine conditioning conditions, then our software may
be useful.

## References

  * Betancourt [2017] M. Betancourt.  A conceptual introduction to Hamiltonian Monte Carlo.  _arXiv preprint arXiv:1701.02434_ , 2017.
  * Kibble and Berkshire [2004] T. Kibble and F. H. Berkshire.  _Classical mechanics_.  World Scientific Publishing Company, 2004.  ISBN 1911298283.
  * Mohasel Afshar and Domke [2015] H. Mohasel Afshar and J. Domke.  Reflection, refraction, and hamiltonian monte carlo.  _Advances in neural information processing systems_ , 28, 2015.  URL https://proceedings.neurips.cc/paper/2015/file/8303a79b1e19a194f1875981be5bdb6f-Paper.pdf.
  * Neal [2011] R. M. Neal.  MCMC using Hamiltonian dynamics.  _Handbook of markov chain monte carlo_ , 2(11):2, 2011.
  * Pakman and Paninski [2014] A. Pakman and L. Paninski.  Exact hamiltonian monte carlo for truncated multivariate gaussians.  _Journal of Computational and Graphical Statistics_ , 23(2):518–542, 2014.  ISSN 1061-8600.  URL https://arxiv.org/abs/1208.4118.
  * Stan Development Team [2023a] Stan Development Team.  Stan Modeling Language Users Guide and Reference Manual, 2.31, 2023a.
  * Stan Development Team [2023b] Stan Development Team.  Stan Functions Reference.  Technical report, 2023b.  URL https://mc-stan.org/docs/2_31/functions-reference.
  * Vishnoi [2021] N. K. Vishnoi.  An introduction to Hamiltonian Monte Carlo method for sampling.  _arXiv preprint arXiv:2108.12107_ , 2021.
  * Windle [2023] J. Windle.  CTGauss, 2023.  URL https://github.com/jwindle/ctgauss.

##  Appendix A Intersection of subspaces

Let

$\ell(x)=\begin{cases}A_{1}^{\prime}x+y_{1},\;f^{\prime}x+g>0\;(R_{1})\\\
A_{2}^{\prime}x+y_{2},\;f^{\prime}x+g<0\;(R_{2})\end{cases}$ (4)

where $f$ is an $n$-dimensional vector, and $A_{1}$ and $A_{2}$ are $n\times
d$ matrices. We want to show that the null spaces of $A_{1}^{\prime}$ and
$A_{2}^{\prime}$ share a common basis $U_{0}$ and are spanned by $U_{0}$ and
single vectors $u_{1}$ and $u_{2}$, respectively.

We can find a common basis easily — just solve $[A_{1}\;A_{2}]^{\prime}v=0.$
But that does not guarantee that there are single vectors $u_{1}$ and $u_{2}$
that suffice for spanning the respective null spaces along with the common
elements. The continuity of $\ell$ is critical for that purpose.

######  Proposition 2.

Suppose $\ell$ is continuous and $f$ is not contained in either the column
space of $A_{1}$ or $A_{2}$ and $A_{1}$ and $A_{2}$ are full rank. Then there
exists a semi-orthonormal matrix $U_{0}$ and unit vectors $u_{1}$ and $u_{2}$,
such that $u_{i}$ and the columns of $U_{0}$ span the null space of
$A_{i}^{\prime}$ and $u_{i}^{\prime}U_{0}=0$, $i=1,2$.

###### Proof.

Make the requisite assumptions and construct $C_{1}=[A_{1}\;f]$. Since $A_{1}$
is full rank and $f$ is not in the column space of $A_{1}$, we know $C_{1}$ is
full rank. Let $Q_{1}R_{1}$ be the QR decomposition of $C$ and let $Q_{0}$ be
the orthogonal complement to $Q_{1}$ so $Q=[Q_{1}\;Q_{0}]$ is unitary. Instead
of working in $x$-space, consider the transformation $x=Qz$ which is
$Q_{1}z_{1}+Q_{0}z_{2}$ after decomposing $z$ into two components. The
contraint $\ell(x)=0$ for region 1 becomes

$A_{1}^{\prime}Q_{1}z_{1}+y_{1}=0$

where $z_{2}$ is free. Similarly, at the boundary $f^{\prime}x+g=0$, we have
$f^{\prime}Q_{1}z_{1}+g=0$ and $z_{2}$ is free. Hence, $z_{1}$ must satisfy

$\displaystyle[A_{1}\;f]^{\prime}Q_{1}z_{1}+b_{1}$ $\displaystyle=0$
$\displaystyle R_{1}^{\prime}z_{1}+b_{1}$ $\displaystyle=0.$

where $b_{1}^{\prime}=[y_{1}^{\prime}\;g]$. The continuity at the boundary
requires that

$A_{2}^{\prime}Qz+y_{2}=A_{2}^{\prime}Q_{1}z_{1}+A_{2}^{\prime}Q_{0}z_{2}+y_{2}=A_{2}^{\prime}Q_{0}z_{2}+(A_{2}^{\prime}Q_{1}z_{1}+y_{2})=0$

for all $z_{2}$. Since this must hold for all $z_{2}\in\mathbb{R}^{n-d-1}$, we
must have both

$\displaystyle A_{2}^{\prime}Q_{0}=0,$ (5) $\displaystyle
A_{2}^{\prime}Q_{1}z_{1}+y_{2}=0.$ (6)

Equation (6) is an algebraic relationship that must be satisfied. Equation (5)
is more critical and it says that $A_{2}^{\prime}$ must have a null space that
includes the span of $Q_{0}$. (Alternatively, $A_{2}$ must have columns that
reside in the span of $Q_{1}$, that is, $A_{2}$ can be constructed using the
columns of $A_{1}$ and $f$.)

Lastly, how do we construct $u_{1}$, $u_{2}$ and $U_{0}$? Equation (5) ensures
that $A_{2}$ is orthogonal to $Q_{0}$, hence $U_{0}=Q_{0}$ characterizes the
shared basis of the null spaces of $A_{1}^{\prime}$ and $A_{2}^{\prime}$. The
QR decomposition of $[A_{1}\;f]$ motivates the rest. In particular, note that
the $(d+1)$th column of $Q_{1}$, call it $u_{1}$ is in the null space of
$A_{1}^{\prime}$. Similarly, let $[A_{2}\;f]=Q_{2}R_{2}$ be the analogous QR
decomposition for $A_{2}$. Then $u_{2}$, the $(d+1)$th column of $Q_{2}$, will
be in the null space of $A_{2}$.

∎

######  Corollary 3.

Given the same suppositions as in Proposition 2, there are unit vectors
$u_{1}$ and $u_{2}$ and semi-orthonormal matrices $U_{0}$ and $U_{c}$ such
that

  1. 1.

The span of $U_{0}$, $u_{1}$, $u_{2}$, and $U_{c}$ is all of $\mathbb{R}^{n}$.

  2. 2.

$U_{0}^{\prime}U_{c}=0$ and $u_{i}^{\prime}U_{0}=u_{i}^{\prime}U_{c}=0$,
$i=1,2$.

  3. 3.

The span of $u_{i}$ and the columns of $U_{0}$ span the null space of
$A_{i}^{\prime}$ and $u_{i}^{\prime}U_{0}=0$, i=1, 2.

  4. 4.

Further, if $A_{1}$ and $A_{2}$ have the same column span, then

    1. (a)

$u_{1}=\pm u_{2}$;

    2. (b)

$[U_{0}\;u_{1}\;U_{c}]$ is an orthonormal matrix.

  5. 5.

if $A_{1}$ and $A_{2}$ do not have the same column span, then

    1. (a)

$U_{c}$ is $n\times(d-1)$;

    2. (b)

$[U_{0}\;u_{b}\;u_{\perp}\;U_{c}]$ — where
$u_{b}=(u_{1}-u_{2})/\|u_{1}-u_{2}\|$ and
$u_{\perp}=(u_{1}+u_{2})/\|u_{1}+u_{2}\|$ — is an orthonormal matrix.

###### Proof.

By Proposition (2), we have computed $U_{0}$, $u_{1}$ and $u_{2}$. Let $U_{c}$
be the matrix representation of an orthonormal basis of the null space of
$[U_{0}\;u_{1}\;u_{2}]^{\prime}$. This decomposition immediately yields (1) -
(5). ∎

##  Appendix B Distribution of a truncated Gaussian conditioned on a
hyperplane

Here we recall some basic facts about a Gaussian distribution and truncated
Gaussian distribution, conditional on the sample lying in some plane.

Suppose that $x\sim N(\mu,\Sigma)$ and we want to find the distribution of
$(x|A^{\prime}x=y)$, when $A$ is full rank. To find the conditional
distribution, we will transform the problem to a different vector of random
variables.

Let $Q=[Q_{1}\;Q_{2}]$ be a unitary matrix such that $Q_{1}$ spans the column
space of $A$ and $Q_{2}$ spans the null space. We want to work with $x=Qz$.
Whether it be by completing the square or using the properties of Gaussian
random variables (and it will matter that we can complete the square later),
the distribution of $z$ is

$\displaystyle z$ $\displaystyle\sim N(\tilde{\mu},\tilde{\Sigma})$
$\displaystyle\tilde{\mu}$ $\displaystyle=Q^{\prime}\mu$
$\displaystyle\tilde{\Sigma}$ $\displaystyle=Q^{\prime}\Sigma Q.$

And, again, by either completing the square or using the properties of
Gaussian random variables, we can decompose the distribution of $z$ as
$p(z_{2}|z_{1})p(z_{1})$ where $\Omega=\tilde{\Sigma}^{-1}$,

$\displaystyle(z_{2}|z_{1})$ $\displaystyle\sim N(\tilde{m}(z_{1}),\tilde{V})$
$\displaystyle\tilde{m}(z_{1})$
$\displaystyle=\tilde{\mu}_{2}+A(z_{1}-\tilde{\mu}_{1})$ $\displaystyle A$
$\displaystyle=\tilde{\Sigma}_{21}\tilde{\Sigma}_{11}^{-1}=\Omega_{22}^{-1}\Omega_{21}$
$\displaystyle\tilde{V}$
$\displaystyle=\tilde{\Sigma}_{22}-A\tilde{\Sigma}_{11}A^{\prime}=\Omega_{22}^{-1}.$

We want to sample $(z|A^{\prime}Qz=y)$, which is to say
$(z|A^{\prime}Q_{1}z_{1}=y)$. In this case, $z_{1}$ is known with certainty
and is $z_{1}=(A^{\prime}Q_{1})^{-1}y$ and we can sample $z_{2}$ from the
conditional distribution above. We can return to the original random vector
$x$ by $x=Qz$ so that

$\displaystyle(x|A^{\prime}x=y)$ $\displaystyle\sim N(m(z_{1}^{*}),V)$
$\displaystyle m(z_{1}^{*})$
$\displaystyle=Q_{1}z_{1}^{*}+Q_{2}(\tilde{\mu}_{2}+A(z_{1}^{*}-\tilde{\mu}_{1}))$
$\displaystyle V$ $\displaystyle=Q_{2}\tilde{V}Q_{2}^{\prime}$ $\displaystyle
z_{1}^{*}$ $\displaystyle=(A^{\prime}Q_{1})^{-1}y.$

##  Appendix C Pseudocode

Data: $M$, $r$, $A$, $y$, mean=False

Result: $x_{p}$, $Q$, $S$

$n\leftarrow\texttt{nrow}(M)$ ;

$d\leftarrow\texttt{ncol}(A)$ ;

$Q,R\leftarrow\texttt{QR}(A,\texttt{mode='complete'})$ ;

$Q_{1}\leftarrow Q[{:},0{:}d]$ ;

$Q_{2}\leftarrow Q[{:},d{:}]$ ;

$R_{1}\leftarrow R[0{:}d,{:}]$ ;

$\Omega_{22}\leftarrow Q_{2}^{\prime}MQ_{2}$ ;

/* $\Omega_{22}=UU^{\prime}\implies Q_{2}\Omega_{22}^{-1}Q_{2}=SS^{\prime}$
where $S=Q_{2}{U^{\prime}}^{-1}$ */

$U\leftarrow\texttt{cholesky($\Omega_{22}$, lower=False)}$ ;

$S^{\prime}\leftarrow\texttt{solve\\_triangular($U$, $Q_{2}^{\prime}$,
trans=$1$, lower=False)}$ ;

/* $U^{\prime}S^{\prime}=Q_{2}^{\prime}$ */

$z_{1}\leftarrow\texttt{solve\\_triangular($R_{1}$, $-y$, trans=$1$,
lower=False)}$ ;

/* $R^{\prime}z_{1}+y=0$ */

$x_{1}\leftarrow Q_{1}z_{1}$ ;

if  _mean_ then

$\tilde{r}\leftarrow M(r-x_{1})$ ;

else

$\tilde{r}\leftarrow r-Mx_{1}$ ;

end if

$x_{p}\leftarrow S(S^{\prime}\tilde{r})+x_{1}$ ;

Algorithm 1 ODEParam

Data: $\phi$, $\mu$, $A$, $y$

Result: $x_{p}$, $Q$, $S$

$n\leftarrow\texttt{ncol}(\mu)$ ;

$d\leftarrow\texttt{nrow}(A)$ ;

$Q,R\leftarrow\texttt{QR}(A,\texttt{mode='complete'})$ ;

$Q_{1}\leftarrow Q[:,0:d]$ ;

$Q_{2}\leftarrow Q[:,d:]$ ;

$S\leftarrow Q_{2}$ ;

$z_{1}\leftarrow\texttt{solve\\_triangular($R$, $-y$, trans=$1$,
lower=False)}$ ;

/* $R^{\prime}z_{1}+y=0$ */

$x_{1}\leftarrow Q_{1}z_{1}$ ;

$\tilde{r}\leftarrow\phi(\mu-x_{1})$ ;

$x_{p}\leftarrow S(S^{\prime}\tilde{r})+x_{1}$ ;

Algorithm 2 IsotropicODEParam: using $M=\phi I_{n}$ and mean $r=\mu$ in Alg

Data: $j$, cache, $M$, $r$, $A$, $y$, mean=False

Result: $x_{p}$, $Q$, $S$

if  _cache $[j]$ is Null_ then

cache[j] = ODEParam($M[j]$, $r[j]$, $A[j]$, $y[j]$, mean)

end if

$(x_{p},Q,S)\leftarrow cache[j]$;

Algorithm 3 GetODEParamForRegion

Data: $j$, $F$, $g$, $L$

Result: $F_{j},g_{j},L_{j}$

$active\leftarrow L[j,{:}]\neq 0$ ;

$signs\leftarrow\texttt{sign}(L[j,active])$ ;

$F_{j}\leftarrow F[active,{:}]\odot signs$ ;

$g_{j}\leftarrow g[active]\odot signs$ ;

$L_{j}\leftarrow L[j,active]$ ;

Algorithm 4 GetBoundariesForRegion

Data: $j$, $M$, $r$, $k$

$V\leftarrow x^{\prime}M[j]x+x^{\prime}r[j]+k[j]$ ;

Algorithm 5 GetPotentialForRegion

Data: $f$, $g$, $A_{1}$, $A_{2}$, $y_{1}$, $y_{2}$, $tol$

Result: $ok$

$B_{1}\leftarrow[A_{1}\;f]$;

$d\leftarrow\texttt{ncol}(B)$;

$(Q,R)\leftarrow\texttt{QR}(B_{1},\texttt{mode='complete'})$;

$Q_{1}\leftarrow Q[{:},0{:}d]$ ;

$Q_{0}\leftarrow Q[{:},d{:}]$ ;

$R_{1}=R[0{:}d,{:}]$;

$z_{1}=\texttt{solve\\_triangular}(R_{1},-y_{1},\texttt{trans=1},\texttt{lower=False})$
;

/* $R_{1}^{\prime}z_{1}+y_{1}=0$ */

$e_{1}\leftarrow\|A_{2}^{\prime}Q_{1}z_{1}+y_{2}\|$ ;

/* $d\times 1$ */

$e_{2}\leftarrow\|A_{2}^{\prime}Q_{0}\|$ ;

/* $d\times(n-d)$ */

$ok\leftarrow(e_{1}<tol)$ and $(e_{2}<tol)$ ;

Algorithm 6 ContinuityCheck

Data: $x_{p}$, $S$, $x_{0}$, $\dot{x}_{0}=\texttt{Null}$

Result: $x_{p}$, $S$

$b\leftarrow x_{0}-x_{p}$ ;

if  _$\dot{x}_{0}$_ then

$a\leftarrow\dot{x}_{0}$

else

$\varepsilon\sim N(0,I_{n-d})$ ;

$a\leftarrow S\varepsilon$ ;

end if

Algorithm 7 ODECoef

Data: $t_{max}$, $x_{p}$, $a$, $b$, $F$, $g$, $L$, $j$

Result: $\tau^{*}$, $j^{*}$, $f$, $x$, $\dot{x}$

$Fa\leftarrow F*a$ ;

$Fb\leftarrow F*b$ ;

$\phi\leftarrow\texttt{arctan2}(-Fa,Fb)$ ;

$h\leftarrow F\;\mu+g$ ;

$u\leftarrow\sqrt{(Fa)^{2}+(Fb)^{2}}$ ;

$\tau\leftarrow[2*t_{max}+2\pi]*J$ ;

/* We have to be careful with $\tau$ */

$\tau\leftarrow\texttt{arccos}(-h/u,\texttt{out}=\tau,\texttt{where}=u\geq|h|)$
;

$\tau\leftarrow\texttt{where}(\tau<0,\tau+\pi,\tau)$ ;

$\tau\leftarrow\texttt{where}(\tau>\pi,\tau-\pi,\tau)$ ;

$dK\leftarrow Fa\odot\cos(\tau)-Fb\odot\sin(\tau)$ ;

$\tau\leftarrow\texttt{where}(dK>0,\tau+\pi,\tau)$ ;

$\tau^{*}=t_{max}$, $j^{*}=j$, $f=a$ ;

/* Default: does not hit boundary */

if  _$\texttt{any}(\tau\leq t_{max})$_ then

;

/* We need a custom argmin to deal with hitting a corner */

$k^{*}\leftarrow\texttt{argmin}(\tau)$ ;

$j^{*}\leftarrow|L[j,k^{*}]|$ ;

/* Need to subtract 1 if using C-indexing */

$\tau^{*}\leftarrow\tau[k^{*}]$ ;

$f=F[k^{*}]$ ;

end if

$x=x_{p}+a\sin(\tau^{*})+b\cos(\tau^{*})$ ;

$\dot{x}=a\cos(\tau^{*})-b\sin(\tau^{*})$ ;

Algorithm 8 EvolveToBoundary

Data: $x$, $\dot{x}$, $u_{1}$

Result: $\dot{x}_{new}$

$v_{1}\leftarrow u_{1}^{\prime}\dot{x}$ ;

$\dot{x}_{new}\leftarrow\dot{x}-2v_{1}u_{1}$ ;

Algorithm 9 WallDynamics

Data: $x$, $\dot{x}$, $j_{1}$, $j_{2}$, $u_{1}$, $u_{2}$, $V_{1}$, $V_{2}$,
check=False

Result: $\dot{x}_{new}$, $j_{new}$

$v_{1}\leftarrow u_{1}^{\prime}\dot{x}$ ;

$E_{v_{1}}\leftarrow\frac{1}{2}v_{1}^{2}$ ;

$\Delta V\leftarrow V_{2}-V_{1}$;

$\dot{x}_{new}\leftarrow\dot{x}-v_{1}u_{1}$ ;

if  _$E_{v_{1}} <\Delta V$_ then

$j_{new}=j_{1}$ ;

$\dot{x}_{new}\leftarrow\dot{x}_{new}-v_{1}u_{1}$ ;

else

$j_{new}=j_{2}$ ;

$\dot{x}_{new}\leftarrow\dot{x}_{new}+\sqrt{2(E_{u_{1}}-\Delta V)}u_{2}$ ;

end if

Algorithm 10 BoundaryDynamics

Data: $f$, $d$, $Q$

Result: $x$

$Q_{1}=Q[{:},0{:}d]$ ;

$x\leftarrow f-Q_{1}Q_{1}^{\prime}f$;

Algorithm 11 resid

Data: $t_{max}$, $j$, $x_{0}$, $\dot{x}_{0}$, $M$, $r$, $A$, $y$, $L$, $F$,
$g$, mean

Result: $x$, $\dot{x}$, $\tau^{*}$, $j^{*}$

$(x_{p},Q,S)\leftarrow$GetODEParam($j$, cache, $M$, $r$, $A$, $y$, mean=False)
;

$a,b\leftarrow$ ODECoef($x_{p}$, $S$, $x_{0}$, $\dot{x}_{0}$) ;

$F_{j},g_{j},L_{j}\leftarrow$ GetBoundariesForRegion ;

$\tau^{*},j^{*},f,x,\dot{x}\leftarrow$ EvolveToBoundary($\tau$, $x_{p}$, $a$,
$b$, $F_{j}$, $g_{j}$, $L_{j}$, $j$) ;

if  _$\tau^{*} >0$_ then

$u_{1}\leftarrow\texttt{resid}(f,d,Q)$ ;

/* $u_{1}=Q[{:},-1]$ from QR of $[A[j]\;f]$ */

if  _$\texttt{hard\\_boundary}\leftarrow j^{*}=j$_ then

$\dot{x}\leftarrow$ WallDynamics($x$, $\dot{x}$, $u_{1}$) ;

else

$(x_{p},Q,S)\leftarrow$GetODEParam($j^{*}$, cache, $M$, $r$, $A$, $y$,
mean=False) ;

$V_{1}\leftarrow V_{j}(x)$ ;

$V_{2}\leftarrow V_{j^{*}}(x)$ ;

$u_{2}\leftarrow\texttt{resid}(f,d,Q)$ ;

/* $u_{1}=Q[{:},-1]$ from QR of $[A[j^{*}]\;f]$ */

$\dot{x}\leftarrow$ BoundaryDynamics($\dot{x}$, $u_{1}$, $u_{2}$, $V_{1}$,
$V_{2}$) ;

end if

end if

Algorithm 12 Evolve1

Data: $t_{max}$, $j$, $x_{0}$, $\dot{x}_{0}$, $M$, $r$, $A$, $y$, $L$, $F$,
$g$, mean

Result: $x$, $\dot{x}$, $\tau^{*}$, $j^{*}$

$(x_{p},S,Q)\leftarrow$GetODEParam($j$, cache, $M$, $r$, $A$, $y$, mean=False)
;

$a,b\leftarrow$ ODECoef($x_{p}$, $S$, $x_{0}$, $\dot{x}_{0}$) ;

$F_{j},g_{j},L_{j}\leftarrow$ GetBoundariesForRegion ;

$\tau^{*},j^{*},f,x,\dot{x}\leftarrow$ EvolveToBoundary($\tau$, $x_{p}$, $a$,
$b$, $F_{j}$, $g_{j}$, $L_{j}$, $j$) ;

$Q_{1}\leftarrow Q[{:},1{:}d]$ ;

$u_{1}\leftarrow f-Q_{1}Q_{1}^{\prime}f$ ;

$(x_{p},S,Q)\leftarrow$GetODEParam($j^{*}$, cache, $M$, $r$, $A$, $y$,
mean=False) ;

$Q_{2}\leftarrow Q[{:},1{:}d]$ ;

$pm\leftarrow 1$ if $(\tau^{*}<t_{max})$ and $(j=j^{*})$ else -1 ;

$u_{2}\leftarrow pm(f-Q_{2}Q_{2}^{\prime}f)$ ;

$V_{1}\leftarrow V_{j}(x)$ ;

$V_{2}\leftarrow V_{j^{*}}(x)$ ;

$\dot{x}\leftarrow$ BoundaryDynamics($\dot{x}$, $u_{1}$, $u_{2}$, $V_{1}$,
$V_{2}$) ;

Algorithm 13 Evolve2

Data: $N$, $t_{max}$, $j$, $x_{0}$, $\dot{x}_{0}$, $M$, $r$, $A$, $y$, $L$,
$F$, $g$, mean

Result: $x$, $\dot{x}$, $t_{remain}$

$n\leftarrow\texttt{len}(x_{0})$ ;

$X\leftarrow\texttt{array}((N,n))$ ;

$\dot{X}\leftarrow\texttt{array}((N,n))$ ;

$R\leftarrow\texttt{array}((N,))$ ;

$x\leftarrow x_{0}$; $\dot{x}\leftarrow\dot{x}_{0}$, $t\leftarrow t_{max}$ ;

for  _$i=0:N$_ do

while  _$t >0$_ do

$x$, $\dot{x}$, $\tau$, $j$ $\leftarrow$ Evolve($t$, $j$, $x$, $\dot{x}$, $M$,
$r$, $A$, $y$, $L$, $F$, $g$, mean) ;

$t\leftarrow t-\tau$ ;

end while

$X[i,{:}]\leftarrow x$ ;

$\dot{X}[i,{:}]\leftarrow\dot{x}$ ;

$R[i]\leftarrow j$ ;

$t\leftarrow t_{max}$ ;

$\dot{x}\leftarrow Null$ ;

end for

Algorithm 14 HMC1

Data: $N$, $t_{max}$, $j$, $x_{0}$, $\dot{x}_{0}$, $M$, $r$, $A$, $y$, $L$,
$F$, $g$, mean

Result: $x$, $\dot{x}$, $t_{remain}$

$n\leftarrow\texttt{len}(x_{0})$ ;

$X\leftarrow\texttt{array}((N,n))$ ;

$\dot{X}\leftarrow\texttt{array}((N,n))$ ;

$R\leftarrow\texttt{array}((N,))$ ;

$I\leftarrow\texttt{array}((N,))$ ;

$x\leftarrow x_{0}$; $\dot{x}\leftarrow\dot{x}_{0}$, $t\leftarrow t_{max}$ ;

for  _$i=0:N$_ do

$x$, $\dot{x}$, $\tau$, $j$ $\leftarrow$ Evolve($t$, $j$, $x$, $\dot{x}$, $M$,
$r$, $A$, $y$, $L$, $F$, $g$, mean) ;

$t\leftarrow t-\tau$ ;

$X[i,{:}]\leftarrow x$ ;

$\dot{X}[i,{:}]\leftarrow\dot{x}$ ;

$R[i]\leftarrow j$ ;

$I[i]\leftarrow t=0$;

if  _ $t=0$ _ then

$t\leftarrow t_{max}$ ;

$\dot{x}\leftarrow Null$ ;

end if

end for

Algorithm 15 HMC2