d}\mu}{{\rm d}\nu}(Z)+\frac{{\rm d}\mu^{\psi_{K}^{-1}}}{{\rm
d}\nu}(Z)\big{)}^{2}\right]-2}\,,$

which can be estimated using empirical averages.

###### Proof of Theorem 18.

We apply Lemma 21 with

$A(Z^{1},\ldots,Z^{N})=\frac{1}{N}\sum_{i=1}^{N}\bar{\mu}_{/\nu}(f\mid
Z^{i})\quad
B(Z^{1},\ldots,Z^{N})=\frac{1}{N}\sum_{j=1}^{N}\bar{\mu}_{/\nu}(\mathrm{1}\mid
Z^{j})\,.$

With $Z^{i}\overset{{\rm iid}}{\sim}\nu$, we have directly

$a=\mathbb{E}\big{(}A\big{)}=\mathbb{E}\big{(}\bar{\mu}_{/\nu}(f\mid
Z)\big{)}=\mu(f)\quad b=\mathbb{E}\big{(}B\big{)}=1\,,$

and

$\displaystyle{\rm var}\big{(}B\big{)}$ $\displaystyle=\frac{1}{N}{\rm
var}\big{(}\bar{\mu}_{/\nu}(1\mid Z^{j})\big{)}$
$\displaystyle=\frac{1}{K^{2}N}\left\\{\mathbb{E}\Bigl{[}\bigl{(}{\textstyle\sum}_{k\in\llbracket
K\rrbracket}{\rm d}\mu^{\psi_{k}^{-1}}/{\rm
d}\nu(Z)\bigr{)}^{2}\Bigr{]}-K^{2}\right\\}\,.$

Now with $\|f\|_{\infty}<\infty$

$\displaystyle{\rm var}\big{(}A\big{)}$ $\displaystyle=\frac{1}{N}{\rm
var}\big{(}\bar{\mu}_{/\nu}(f\mid Z)\big{)}$
$\displaystyle=\frac{1}{K^{2}N}\Bigl{\\{}\|f\|_{\infty}^{2}\mathbb{E}\Bigl{[}\bigl{(}{\textstyle\sum}_{k\in\llbracket
K\rrbracket}{\rm d}\mu^{\psi_{k}^{-1}}/{\rm
d}\nu(Z)\bigr{)}^{2}\Bigr{]}-K^{2}\mu(f)^{2}\Bigr{\\}}$
$\displaystyle\leq\frac{\|f\|_{\infty}^{2}}{K^{2}N}\mathbb{E}\Bigl{[}\bigl{(}{\textstyle\sum}_{k\in\llbracket
K\rrbracket}{\rm d}\mu^{\psi_{k}^{-1}}/{\rm d}\nu(Z)\bigr{)}^{2}\Bigr{]}\,.$

We conclude noting that we have $|A/B|\leq\|f\|_{\infty}$. ∎

For clarity we reproduce the very useful lemma [38, Lemma 6], correcting a
couple of minor typos along the way.

######  Lemma 21.

Let $A,B$ be two integrable random variables satisfying $|A/B|\leq M$ almost
surely for some $M>0$ and let $a=\mathbb{E}(A)$ and $b=\mathbb{E}(B)\neq 0$.
Then

$\displaystyle{\rm
var}(A/B)\leq\mathbb{E}[|A/B-a/b|^{2}]\leq\frac{2}{b^{2}}\Bigl{\\{}{\rm
var}(A)+M^{2}{\rm var}(B)\Bigr{\\}}\,,$
$\displaystyle|\mathbb{E}[A/B]-a/b|\leq\frac{\sqrt{{\rm var}(A/B){\rm
var}(B)}}{b}\,.$

####  4.2.2 Optimal weights

As mentioned in Subsection 2.3 it is possible to consider the more general
scenario where unequal probability weights
$\omega=\\{\omega_{k}\in\mathbb{R},k\in\llbracket
0,T\rrbracket\colon\sum_{k=0}^{T}\omega_{k}=1\\}$ are ascribed to the elements
of the snippets, yielding the estimator

$\displaystyle\hat{\mu}_{\omega}(f)$
$\displaystyle:=\frac{\sum_{i=1}^{N}\sum_{k=0}^{T}\omega_{k}\frac{{\rm
d}\mu^{\psi_{k}^{-1}}}{{\rm
d}\nu}(Z^{i})f\circ\psi_{k}(Z^{i})}{\sum_{j=1}^{N}\sum_{l=0}^{T}\omega_{l}\frac{{\rm
d}\mu^{\psi_{l}^{-1}}}{{\rm d}\nu}(Z^{j})}\,,$

and a natural question is that of the optimal choice of $\omega$. Note that in
the context of PISA the condition $\omega_{k}\geq 0$ is not required, as
suggested by the justification of the identity (30). However it should be
clear that this condition should be enforced in the context of integrator
snippet SMC, since the probabilistic interpretation is otherwise lost, or if
the expectation is known to be non-negative. Here we discuss optimization of
the variance upperbound provided by Lemma 21,

$\displaystyle\frac{2\|f\|_{\infty}^{2}}{N}\Bigl{\\{}{\rm
var}\biggl{(}{\textstyle\sum_{k\in\llbracket K\rrbracket}}\omega_{k}\frac{{\rm
d}\mu^{\psi_{k}^{-1}}}{{\rm
d}\nu}(Z)\frac{f\circ\psi_{k}}{\|f\|_{\infty}}(Z)\biggr{)}+$
$\displaystyle{\rm var}\biggl{(}{\textstyle\sum_{k\in\llbracket
K\rrbracket}}\omega_{k}\frac{{\rm d}\mu^{\psi_{k}^{-1}}}{{\rm
d}\nu}(Z)\biggr{)}\Bigr{\\}}$ $\displaystyle\leq\frac{2\|f\|_{\infty}^{2}}{N}$
$\displaystyle\omega^{\top}\big{(}\Sigma(f;\psi)+\Sigma_{\psi}(1;\psi)\big{)}\omega$

where for $k,l\in\llbracket K\rrbracket$,

$\Sigma_{kl}(f;\psi)=\mathbb{E}\biggl{[}\frac{{\rm d}\mu^{\psi_{k}^{-1}}}{{\rm
d}\nu}(Z)\frac{f\circ\psi_{k}}{\|f\|_{\infty}}(Z)\frac{{\rm
d}\mu^{\psi_{l}^{-1}}}{{\rm
d}\nu}(Z)\frac{f\circ\psi_{l}}{\|f\|_{\infty}}(Z)\biggr{]}-\mu(f/\|f\|_{\infty})^{2}\,.$

It is a classical result that
$\omega^{\top}\big{(}\Sigma(f;\psi)+\Sigma_{\psi}(1;\psi)\big{)}\omega\geq\lambda_{\min}\big{(}\Sigma(f;\psi)+\Sigma_{\psi}(1;\psi)\big{)}\omega^{\top}\omega$
and that minimum is reached for the eigenvector(s) $\omega_{{\rm min}}$
corresponding to the smallest eigenvalue
$\lambda_{\min}\big{(}\Sigma(f;\psi)+\Sigma_{\psi}(1;\psi)\big{)}$ of
$\Sigma(f;\psi)+\Sigma_{\psi}(1;\psi)$. If the constraint of non-negative
entries is to be enforced then a standard quadratic programming procedure
should be used. The same ideas can be applied to the function independent
upperbounds used in our definitions of efficiency.

###  4.3 More on variance reduction and optimal flow

We now focus on some properties of the estimator $\hat{\mu}(f)$ of $\mu(f)$ in
(16). To facilitate the analysis and later developments we consider the
scenario where, with $t\mapsto\psi_{t}(z)$ the flow solution of an ODE, assume
the dominating measure $\upsilon\gg\mu$ and $\upsilon$ is invariant by the
flow. We consider

$\bar{\mu}({\rm d}t,{\rm d}z)=\frac{1}{T}\mu^{\psi_{-t}}({\rm
d}z)\mathbf{1}\\{0\leq t\leq T\\}{\rm d}t\,,$

and notice that similarly to the integrator scenario, for any
$f\colon\mathsf{Z}\rightarrow\mathbb{R}$, $\bar{f}(t,z):=f\circ\psi_{t}(z)$ we
have
$\mathbb{E}_{\bar{\mu}}\left(\bar{\mu}(\bar{f}\mid\check{Z})\right)=\bar{\mu}(\bar{f})=\mu(f)$.
where for any $z\in\mathsf{Z}$, where

$\bar{\mu}(\bar{f}\mid
z)=\frac{1}{T}\int_{0}^{T}f\circ\psi_{t}(z)\frac{\mu\circ\psi_{t}(z)}{\int\mu\circ\psi_{u}(z){\rm
d}u}{\rm d}t$

the following estimator of $\mu(f)$ is considered, with
$\check{Z}_{i}\overset{{\rm iid}}{\sim}\bar{\mu}$

$\hat{\mu}(f)=\frac{1}{N}\sum_{i=1}^{N}\bar{\mu}\big{(}\bar{f}(\mathrm{T},\check{Z}_{i})\mid\check{Z}_{i}\big{)}\,.$

We now show that in the examples considered in this paper our approach can be
understood as implementing unbiased Riemann sum approximations of line
integrals along contours. Adopting a continuous time approach can be justified
as follows. For numerous integrators, the conditions of the following
proposition are satisfied; this is the case for the leapfrog integrator of
Hamilton’s equation e.g. [20, 38, Theorem 3.4] and [38, Appendix 3.1, Theorem
9] for detailed results.

######  Proposition 22.

Let $\tau>0$, for any $z\in\mathsf{Z}$ let $[0,\tau]\ni t\mapsto\psi_{t}(z)$
be a flow and for $\epsilon>0$ let $\\{\psi^{k}(z;\epsilon),0\leq
k\epsilon\leq\tau\\}$ be a discretization of $t\mapsto\psi_{t}$ such that that
for any $z\in\mathsf{Z}$ there exists $C>0$ such that for any
$(k,\epsilon)\in\mathbb{N}\times\mathbb{R}_{+}$

$|\psi^{k}(z;\epsilon)-\psi_{k\epsilon}(z)|\leq C\epsilon^{2}\,.$

Then for any continuous $g\colon\mathsf{Z}\rightarrow\mathbb{R}$ such that the
Riemann integral

$I(g):=\int_{0}^{\tau}g\circ\psi_{t}(z){\rm{\rm d}}t\,,$

exists we have

$\lim_{T\rightarrow\infty}\frac{1}{n}\sum_{k=0}^{n-1}g\circ\psi_{k\tau/n}(z;\tau/n)=I(g)\,.$

###### Proof.

We have

$\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{k=0}^{n-1}g\circ\psi_{k\tau/n}(z)=\int_{0}^{\tau}g\circ\psi_{t}(z){\rm{\rm
d}}t$
$\lim_{n\rightarrow\infty}\frac{1}{n}\biggl{|}\sum_{k=0}^{n-1}g\circ\psi_{k\tau/n}(z)-g\circ\psi^{k}(z;\tau/n)\biggr{|}=0\,,$

and we can conclude. ∎

######  Remark 23.

Naturally convergence is uniform in $z\in\mathsf{Z}$ with additional
assumptions and we note that in some scenarios dependence on $\tau$ can be
characterized, allowing in principle $\tau$ to grow with $n$ .

####  4.3.1 Hamiltonian contour decomposition

Assume $\mu$ has a probability density with respect to the Lebesgue measure
and let $\zeta\colon\mathsf{Z}\subset\mathbb{R}^{d}\rightarrow\mathbb{R}$ be
Lipschitz continuous such that $\nabla\zeta(z)\neq 0$ for all
$z\in\mathsf{Z}$, then the co-area formula states that

$\int_{\mathsf{Z}}f(z)\mu(z){\rm
d}z=\int\big{[}\int_{\zeta^{-1}(s)}f(z)|\nabla\zeta(z)|^{-1}\mu(z)\mathcal{H}_{d-1}({\rm
d}z)\big{]}{\rm d}s,$

where $\mathcal{H}_{d-1}$ is the $(d-1)$-dimensional Hausdorff measure, used
here to measure length along the contours $\zeta^{-1}(s)$. For example in the
HMC setup where $\mu(z)\propto\exp\big{(}-H(z)\big{)}$ and $z=(x,v)$ one may
choose $\zeta(z)=H(z)=-\log\big{(}\mu(z)\big{)}$, leading to a decomposition
of the expectation $\mu(f)$ according to equi-energy contours of $H(z)$

$\mu(f)=\int\exp(-s)\big{[}\int_{\zeta^{-1}(s)}f(z)|\nabla
H(z)|^{-1}\mathcal{H}_{d-1}({\rm d}z)\big{]}{\rm d}s.$

We now show how the solution of Hamilton’s equation could be used as the basis
for estimators of $\mu(f)$ mixing Riemanian sum-like and Monte Carlo
estimation techniques.

##### Favourable scenario where $d=2$.

Let $s\in\mathbb{R}$ such that $H^{-1}(s)\neq\emptyset$ and assume that for
some $(x_{0},v_{0})\in H^{-1}(s)$ Hamilton’s equations
$\dot{x}=-\nabla_{v}H(x,v)$ and $\dot{v}=\nabla_{x}H(x,v)$ have solutions
$t\mapsto(x_{t},v_{t})=\psi_{t}(z_{0})\in H^{-1}(s)$ with
$\big{(}x_{\tau(s)},v_{\tau(s)})=(x_{0},v_{0})=z_{0}$ for some $\tau(s)>0$,
that is the contours $H^{-1}(s)$ can be parametrised with the solutions of
Hamilton’s equation at corresponding level.

######  Proposition 24.

We have

$\int_{\mathsf{Z}}f(z)\mu({\rm d}z)=\int_{\mathsf{Z}}\left[[\tau\circ
H(z_{0})]^{-1}\int_{0}^{\tau\circ H(z_{0})}f(z_{t}){\rm d}t\right]\mu({\rm
d}z_{0})\,,$

implying that, assuming the integral along the path $t\mapsto z_{t}$ is
tractable,

$[\tau\circ H(z_{0})]^{-1}\int_{0}^{\tau\circ H(z_{0})}f(z_{t}){\rm
d}t\,\text{for}\,z_{0}\sim\mu\,,$ (36)

is an unbiased estimator of $\mu(f)$.

###### Proof.

Since $H^{-1}(s)$ is rectifiable we have that [9, Theorem 2.6.2] for
$g\colon\mathsf{Z}\rightarrow\mathbb{R}$,

$\int_{H^{-1}(s)}g(z)\mathcal{H}_{d-1}({\rm
d}z)=\int_{0}^{\tau(s)}g(z_{t})|\nabla H(z_{t})|{\rm d}t\,.$

Consequently

$\displaystyle\int_{\mathsf{Z}}f(z)\mu({\rm d}z)$
$\displaystyle=\int\tau(s)\exp(-s)\left[\tau(s)^{-1}\int_{0}^{\tau(s)}f(z_{s,t}){\rm
d}t\right]{\rm d}s\,,$ (37)

where the notation now reflects dependence on $s$ of the flow and notice that
since

$H^{-1}(s)\ni z_{0}\mapsto[\tau\circ H(z_{0})]^{-1}\int_{0}^{\tau\circ
H(z_{0})}f(z_{t})$

is constant, $z_{0,s}\in H^{-1}(s)$ can be arbitrary. since indeed from (37)

$\displaystyle\int_{\mathsf{Z}}\mathbf{1}\\{H(z)\in A\\}\mu({\rm d}z)=$
$\displaystyle\int\tau(s)\exp(-s)\left[\tau(s)^{-1}\int_{0}^{\tau(s)}\mathbf{1}\\{s\in
A\\}{\rm d}t\right]{\rm d}s$ $\displaystyle=$
$\displaystyle\int\tau(s)\exp(-s)\mathbf{1}\\{s\in A\\}{\rm d}s.$

Note that we can write, since for $s\in\mathbb{R}$, $H^{-1}(s)\ni
z_{0}\mapsto[\tau\circ H(z_{0})]^{-1}\int_{0}^{\tau\circ H(z_{0})}f(z_{t})$ is
constant,

$\displaystyle\int_{\mathsf{Z}}\left[[\tau\circ
H(z_{0})]^{-1}\int_{0}^{\tau\circ H(z_{0})}f(z_{t}){\rm d}t\right]\mu({\rm
d}z_{0})$ $\displaystyle=\int\tau(s)\exp(-s)$
$\displaystyle\left\\{\tau(s)^{-1}\int_{0}^{\tau(s)}\left[\big{(}\tau\circ
H(z_{s,0})\big{)}^{-1}\int_{0}^{\tau\circ H(z_{0})}f(z_{s,t}){\rm
d}t\right]\right\\}{\rm d}t{\rm d}s$ $\displaystyle=\int\tau(s)\exp(-s)$
$\displaystyle\left[\tau(s)^{-1}\int_{0}^{\tau(s)}f(z_{s,t}){\rm
d}t\right]{\rm d}t{\rm d}s\,.$

∎

A remarkable point is that the strategy developed in this manuscript provides
a general methodology to implement numerically the ideas underpinning the
decomposition of Proposition 24 by using the estimator $\check{\mu}(f)$ in
(15) and invoking Proposition 22, assuming $s\mapsto\tau(s)$ to be known. This
point is valid outside of the SMC framework and it is worth pointing out that
$\check{\mu}(f)$ in (15) is unbiased if the samples used are marginally from
$\bar{\mu}$.

Remark the promise of dimension free estimators if the one dimensional line
integrals in (36) were easily computed and sampling from the one-dimensional
energy distribution was routine – however the scenario $d\geq 3$ is more
subtle.

##### General scenario

In the scenario where $d\geq 3$ the co-area decomposition still holds, but the
solution to Hamilton’s equation can in general not be used to compute
integrals over the hypersurface $H^{-1}(s)$. This would require a form of
ergodicity [35, 40] of the form, for $z_{0}\in H^{-1}(s)$,

$\lim_{\tau\rightarrow\infty}\frac{1}{\tau}\int_{0}^{\tau}f\circ\psi_{t}(z_{0}){\rm
d}t=\bar{f}(z_{0})=\frac{\int_{\zeta^{-1}(s)}f(z)|\nabla
H(z)|^{-1}\mathcal{H}_{d-1}({\rm d}z)}{\int_{\zeta^{-1}(s)}|\nabla
H(z)|^{-1}\mathcal{H}_{d-1}({\rm d}z)}\,,$

where the limit always exists in the ${\rm L}^{2}(\mu)$ sense and constitutes
Von Neumann’s mean ergodic theorem [41, 35], and the rightmost equality forms
Boltzman’s conjecture. An interesting property in the present context is that
$\mathbb{E}_{\bar{\mu}}\bigl{(}\bar{f}(Z)\bigr{)}=\mathbb{E}_{\mu}\bigl{(}f(Z)\bigr{)}$
for $f\in\mathrm{L}^{2}(\mu)$ and one could replicate the variance reduction
properties developed earlier. Boltzman’s conjecture has long been disproved,
as numerous Hamiltonians can be shown not to lead to ergodic systems, although
some sub-classes do. However a weaker, or local, form or ergodicity can hold
on sets of a partition of $\mathsf{Z}$

######  Example 25 (Double well potential).

Consider the scenario where $\mathsf{X}=\mathsf{V}=\mathbb{R}$,
$U(x)=(x^{2}-1)^{2}$ and kinetic energy $v^{2}$ [34]. Elementary manipulations
show that satisfying Hamilton’s equations (2) imposes $t\mapsto
H(x_{t},\dot{x}_{t})=(x_{t}^{2}-1)^{2}+\dot{x}_{t}^{2}=s>0$ and therefore
requires $t\mapsto
x_{t}\in[-\sqrt{1+\sqrt{s}},-\sqrt{1-\sqrt{s}}]\cup[\sqrt{1-\sqrt{s}},\sqrt{1+\sqrt{s}}]$
– importantly the intervals are not connected for $s<1$. Rearranging terms any
solution of (2) must satisfy

$\displaystyle\dot{x}_{t}$ $\displaystyle=\pm\sqrt{s-(x_{t}^{2}-1)^{2}}\,,$

that is the velocity is a function of the position in the double well, maximal
for $x_{t}^{2}=1$, vanishing as $x_{t}^{2}\rightarrow 1\pm\sqrt{s}$ and a sign
flip at the endpoints of the intervals. Therefore the system is not ergodic,
but ergodicity trivially occurs in each well.

In general this partitioning of $\mathsf{Z}$ can be intricate but it should be
clear that in principle this could reduce variability of an estimator. In the
toy Example 25, a purely mathematical algorithm inspired by the discussions
above would choose the right or left well with probability $1/2$ and then
integrate deterministically, producing samples taking at most two values which
could be averaged to compute $\mu(f)$. We further note that in our context a
relevant result would be concerned with the limit, for any
$x_{0}\in\mathsf{X}$,

$\lim_{\tau\rightarrow\infty}\int\big{[}\frac{1}{\tau}\int_{0}^{\tau}f\circ\psi_{t}(x_{0},\rho\mathbf{e}){\rm
d}t\big{]}\varpi_{\mathbf{e}}({\rm d}\mathbf{e})\,,$

where we have used a polar reparametrization $v=\rho\mathbf{e}$
$(\rho,\mathbf{e})\in\mathbb{R}_{+}\times\mathcal{S}_{d-1}$, and whether a
marginal version of Boltzman’s conjecture holds for this limit. Example 25
indicates that this is not true in general but the cardinality of the
partition of $\mathsf{X}$ may be reduced.

####  4.3.2 Advantage of averaging and control variate interpretation

Consider the scenario where $(\mathsf{T},\check{Z})\sim\bar{\mu}(t,{\rm
d}z)=\frac{1}{\tau}\mu^{\psi_{-t}}({\rm d}z)\mathbf{1}\\{0\leq t\leq\tau\\}$,
$[0,\tau]\ni t\mapsto\psi_{t}$ for some $\tau>0$ is the flow solution of an
ODE of the form $\dot{z}_{t}=F\circ z_{t}$ for some field
$F\colon\mathsf{Z}\rightarrow\mathsf{Z}$. For
$f\colon\mathsf{Z}\rightarrow\mathbb{R}$ we have

${\rm var}_{\mu}\left(f(Z)\right)={\rm
var}_{\bar{\mu}}\left(\mathbb{E}_{\bar{\mu}}\big{(}\bar{f}(\mathsf{T},\check{Z})\mid\check{Z}\big{)}\right)+\mathbb{E}_{\bar{\mu}}\big{(}{\rm
var}_{\bar{\mu}}(\bar{f}(\mathsf{T},\check{Z})\mid\check{Z})\big{)}\,.$

and we are interested in determining $t\mapsto\psi_{t}$ (or $F$) in order to
minimize ${\rm
var}_{\bar{\mu}}\left(\mathbb{E}_{\bar{\mu}}\big{(}\bar{f}(\mathsf{T},\check{Z})\mid\check{Z}\big{)}\right)$
and maximize improvement over simple Monte Carlo. We recall that the Jabobian
determinant of the flow $t\mapsto\psi_{t}(z)$ is given by [19, p. 174108-5]

$J_{t}(z)=|\det\big{(}\nabla\otimes\psi_{t}(z)\big{)}|=\mathcal{J}\circ\psi_{t}(z)\;\text{with}\;\mathcal{J}(z):=\exp\left(\int_{0}^{t}(\nabla\cdot
F)(z)\right)\,.$

######  Lemma 26.

Let $t\mapsto\psi_{t}$ be a flow solution of $\dot{z}_{t}=F\circ z_{t}$ and
assume that with $\upsilon$ the Lebesgue measure, $\upsilon\gg\mu$. Then

$\displaystyle\mathbb{E}_{\bar{\mu}}\left(\mathbb{E}_{\bar{\mu}}\big{(}\bar{f}(\mathsf{T},\check{Z})\mid\check{Z}\big{)}^{2}\right)$
$\displaystyle=2\int_{0}^{\tau}\left\\{\int_{0}^{\tau-u}[\bar{\mu}\circ\psi_{-t}(z)]^{-1}{\rm
d}t\right\\}\Bigl{\\{}\int\mu({\rm
d}z)f(z)f\circ\psi_{u}(z)\mu\circ\psi_{u}(z)\mathcal{J}\circ\psi_{u}(z)\Bigr{\\}}{\rm
d}u$ $\displaystyle\leq
2\left\\{\int_{0}^{\tau}[\bar{\mu}\circ\psi_{-t}(z)]^{-1}{\rm
d}t\right\\}\int_{0}^{\tau}\Bigl{\\{}\int\mu({\rm
d}z)f(z)f\circ\psi_{u}(z)\mu\circ\psi_{u}(z)\mathcal{J}\circ\psi_{u}(z)\Bigr{\\}}{\rm
d}u\,,$

with

$\displaystyle\bar{\mu}\circ\psi_{-t}(z)$
$\displaystyle=\int_{-t}^{\tau-t}\mu\circ\psi_{u}(z)\,\mathcal{J}\circ\psi_{u}(z){\rm
d}u\,.$

In particular for $t\mapsto\psi_{t}$ the flow solution of Hamilton’s equations
associated with $\mu$ we have

$\displaystyle\mathbb{E}_{\bar{\mu}}\left(\mathbb{E}_{\bar{\mu}}\big{(}\bar{f}(\mathsf{T},\check{Z})\mid\check{Z}\big{)}^{2}\right)$
$\displaystyle=2\int\big{(}1-t/\tau\big{)}\int\mu({\rm
d}z)f(z)f\circ\psi_{t}(z){\rm d}t\,\cdot$

###### Proof.

Using Fubini’s theorem we have

$\displaystyle\int\left\\{\int_{0}^{\tau}f\circ\psi_{t}(z)\frac{{\rm
d}\mu^{\psi_{-t}}}{{\rm d}\bar{\mu}}(z){\rm d}t\right\\}^{2}\bar{\mu}({\rm
d}z)$ $\displaystyle=$
$\displaystyle=2\int_{0}^{\tau}\int_{0}^{\tau}\mathbf{1}\\{t^{\prime}\geq
t\\}\int\bar{\mu}({\rm d}z)$ $\displaystyle f\circ\psi_{t}(z)\frac{{\rm
d}\mu^{\psi_{-t}}}{{\rm
d}\bar{\mu}}(z)f\circ\psi_{t^{\prime}-t+t}(z)\frac{{\rm
d}\mu^{\psi_{-t^{\prime}}}}{{\rm
d}\bar{\mu}}\circ\psi_{-t}\circ\psi_{t}(z){\rm d}t^{\prime}{\rm d}t$
$\displaystyle=2\int_{0}^{\tau}\int_{0}^{\tau}\mathbf{1}\\{t^{\prime}\geq
t\\}\int\mu({\rm d}z)$ $\displaystyle
f(z)f\circ\psi_{t^{\prime}-t}(z)\frac{{\rm d}\mu^{\psi_{-t^{\prime}}}}{{\rm
d}\bar{\mu}}\circ\psi_{-t}(z){\rm d}t^{\prime}{\rm d}t\,.$

Further, we have

$\displaystyle\frac{{\rm d}\mu^{\psi_{-t^{\prime}}}}{{\rm d}\bar{\mu}}(z)$
$\displaystyle=\frac{\mu\circ\psi_{t^{\prime}}(z)\,\mathcal{J}\circ\psi_{t^{\prime}}(z)}{\bar{\mu}(z)}\,,$

with
$\bar{\mu}(z)=\int_{0}^{\tau}\mu\circ\psi_{u}(z)\,\mathcal{J}\circ\psi_{u}(z){\rm
d}u$. It is straightforward that

$\displaystyle\bar{\mu}\circ\psi_{-t}(z)$
$\displaystyle=\int_{0}^{\tau}\mu\circ\psi_{u-t}(z)\,\mathcal{J}\circ\psi_{u-t}(z){\rm
d}u$
$\displaystyle=\int_{-t}^{\tau-t}\mu\circ\psi_{u^{\prime}}(z)\,\mathcal{J}\circ\psi_{u^{\prime}}(z){\rm
d}u^{\prime}\,.$

Consequently

$\displaystyle\int\left\\{\int_{0}^{\tau}f\circ\psi_{t}(z)\frac{{\rm
d}\mu^{\psi_{-t}}}{{\rm d}\bar{\mu}}(z){\rm d}t\right\\}^{2}\bar{\mu}({\rm
d}z)=$ $\displaystyle=2\int_{0}^{\tau}\int_{0}^{\tau}$
$\displaystyle\frac{\mathbf{1}\\{\tau-t\geq u\geq
0\\}}{\bar{\mu}\circ\psi_{-t}(z)}\int\mu({\rm
d}z)f(z)f\circ\psi_{u}(z)\mu\circ\psi_{u}(z)\mathcal{J}\circ\psi_{u}(z){\rm
d}u{\rm d}t$ $\displaystyle=2\int_{0}^{\tau}$
$\displaystyle\left\\{\int_{0}^{\tau-u}[\bar{\mu}\circ\psi_{-t}(z)]^{-1}{\rm
d}t\right\\}\int\mu({\rm
d}z)f(z)f\circ\psi_{u}(z)\mu\circ\psi_{u}(z)\mathcal{J}\circ\psi_{u}(z){\rm
d}u$

For the second statement we have $\bar{\mu}\circ\psi_{-t}(z)=\tau\,\mu(z)$ and

$\displaystyle\int\left\\{\int_{0}^{\tau}f\circ\psi_{t}(z)\frac{{\rm
d}\mu^{\psi_{-t}}}{{\rm d}\bar{\mu}}(z){\rm d}t\right\\}^{2}\bar{\mu}({\rm
d}z)=$ $\displaystyle=2\int_{0}^{\tau}$
$\displaystyle\bigl{(}1-u/\tau\bigr{)}\int\mu({\rm
d}z)f(z)f\circ\psi_{u}(z){\rm d}t^{\prime\prime}$

which is akin to the integration correlation time encountered in MCMC. ∎

This is somewhat reminiscent of what is advocated in the literature in the
context of HMC or randomized HMC where the integration time $t$ (or $T$ when
using an integrator) is randomized [29, 21].

######  Example 27.

Consider the Gaussian scenario where $\mathsf{X}=\mathsf{V}=\mathbb{R}^{d}$,
$\pi(x)=\mathcal{N}(x;0,\Sigma)$ with diagonal covariance matrix such that for
$i\in\llbracket d\rrbracket$, $\Sigma_{ii}=\sigma_{i}^{2}$ and
$\varpi(v)=\mathcal{N}(v;0,{\rm Id})$. Using the reparametrization
$(x_{0}(i),v_{0}(i))=(a_{i}\sigma_{i}\sin(\phi_{i}),a_{i}\cos\phi_{i})$ the
solution of Hamilton’s equations is given for $i\in\llbracket d\rrbracket$ by

$\displaystyle x_{t}(i)$
$\displaystyle=a_{i}\sigma\sin(t/\sigma_{i}+\phi_{i})$ $\displaystyle
v_{t}(i)$ $\displaystyle=a_{i}\cos(t/\sigma_{i}+\phi_{i})\,.$

We have for each component
$\mathbb{E}_{\mu}[X_{0}(i)X_{t}(i)]=\sigma_{i}^{2}\cos(t/\sigma_{i})$, which
clearly does not vanish as $t$, or $\tau$, increase: this is a particularly
negative result for standard HMC where the final state of the computed
integrator is used. Worse, as noted early on in [29] this implies that for
heterogeneous values $\\{\sigma_{i}^{2},i\in\llbracket d\rrbracket\\}$ no
integration time may be suitable simultaneously for all coordinates. This
motivated the introduction of random integration times [29] which leads to the
average correlation

$\mathbb{E}_{\mu}\left\\{\frac{1}{\tau}\int_{0}^{\tau}X_{0}(i)X_{t}(i){\rm
d}t\right\\}=\frac{\sin(\tau/\sigma_{i})}{\tau/\sigma_{i}}\,,$

where it is assumed that $\tau$ is independent of the initial state. This
should be contrasted with the fixed integration time scenario since as $\tau$
increases this vanishes and is even negative for some values of $\tau$ (the
minimum is here reached for about $\tau_{i}=4.5\sigma_{i}$ for a given
component).

The example therefore illustrates that our approach implements this averaging
feature, and therefore shares its benefits, within the context of an iterative
algorithm. The example also highlights a control variate technique
intepretation. More specifically in the discrete time scenarios
$\\{f\circ\psi^{k}(z),k\in\llbracket P\rrbracket\\}$ can be interpreted as
control variates, but can induce both positive and negative correlations.

####  4.3.3 Towards an optimal flow?

In this section we are looking to determine a flow for some $\tau>0$
$[0,\tau]\ni t\mapsto\psi_{t}$ of an ODE of the form $\dot{z}_{t}=F\circ
z_{t}$ for some field $F\colon\mathsf{Z}\rightarrow\mathsf{Z}$ which defines
as above the probability model

$\bar{\mu}({\rm d}t,{\rm d}z)=\frac{1}{\tau}\mu^{\psi_{-t}}({\rm
d}z)\mathbf{1}\\{0\leq t\leq\tau\\}{\rm d}t\,,$

which has the property that for any $f\colon\mathsf{Z}\rightarrow\mathbb{R}$,
defining $\bar{f}(t,z):=f\circ\psi_{t}(z)$, then $\bar{\mu}(\bar{f})=\mu(f)$.
This suggests the use of a Rao-Blackwellized estimator inspired by
$\mathbb{E}_{\bar{\mu}}\left(\bar{\mu}(\bar{f}\mid\check{Z})\right)$. Assuming
$\mathsf{Z}=\mathbb{R}^{d}\times\mathbb{R}^{d}$ and that $\mu$ has a density
with respect to the Lebesgue measure, then for $z\in\mathsf{Z}$

$\bar{\mu}(\bar{f}\mid
z)=\frac{1}{\tau}\int_{0}^{\tau}f\circ\psi_{t}(z)\frac{\mu\circ\psi_{t}(z)\,\mathcal{J}\circ\psi_{t}(z)}{\int\mu\circ\psi_{u}(z)\,\mathcal{J}\circ\psi_{u}{\rm
d}u}{\rm d}t$ (38)

In the light of Lemma 26 we aim to find for any $z\in\mathsf{Z}$ the flow
solutions $t\mapsto\psi_{t}(z)$ of ODEs $\dot{z}_{t}=F_{z}(z_{t})$ such that
the function $t\mapsto\int\mu({\rm
d}z)f(z)f\circ\psi_{t}(z)\mu\circ\psi_{t}(z)\,\mathcal{J}\circ\psi_{t}(z)$
decreases as fast as possible. This is motivated by the fact that the integral
on $[0,\tau]$ of this mapping appears in the variance upper bound for (38) in
Lemma 26, which we want to minimize. Note that we also expect this mapping to
be smooth under general conditions not detailed in this preliminary work. For
smooth enough flow and $f$ we have, with $g:=f\times\mu\times\mathcal{J}$,

$\displaystyle\frac{{\rm d}}{{\rm
d}t}\mathbb{E}_{\mu}\left[f(Z)\frac{g\circ\psi_{t}(Z)}{\bar{\mu}(Z)}\right]$
$\displaystyle=\mathbb{E}_{\mu}\left[\frac{f(Z)}{\bar{\mu}(Z)}\langle\nabla
g\circ\psi_{t}(Z),\dot{\psi}_{t}(Z)\rangle\right]$

Pointwise, the steepest descent direction is given by

$\dot{\psi}_{t}(z)=-C_{t}(z)\frac{f(z)\,\nabla
g\circ\psi_{t}(z)}{|f(z)\,\nabla g\circ\psi_{t}(z)|}$

for a positive function $C_{t}(z)$ to be determined optimally. In this
scenario we therefore have

$\frac{{\rm d}}{{\rm
d}t}\mathbb{E}_{\mu}\left[f(Z)\frac{g\circ\psi_{t}(Z)}{\bar{\mu}(Z)}\right]=-\mathbb{E}_{\mu}\left[\frac{C_{t}(z)}{\bar{\mu}(Z)}|f(Z)\,\nabla
g\circ\psi_{t}(Z)|\right]$

and by Cauchy-Schwartz the (positive) expectation is maximized for

$C_{t}(z)=C_{t}\frac{|f(z)\,\nabla g\circ\psi_{t}(z)|}{\bar{\mu}(z)}\,.$

and the trajectories we are interested in must be such that, for some $C>0$,

$\dot{\psi}_{t}(z)=-\frac{C}{\bar{\mu}(z)}f(z)\,\nabla
g\circ\psi_{t}(z)=-\frac{C}{\bar{\mu}(z)}f(z)F\circ\psi_{t}(z)\,.$

Note that for any $z$ the term $|f(z)/\bar{\mu}(z)|$ is only a change of speed
and that the trajectory of $\mathbb{R}_{+}\ni t\mapsto\psi_{t}(z)$ is
independent of this factor, despite the remarkable fact that $\bar{\mu}(z)$
depends on this flow. The result seems fairly natural.

##  5 MCMC with integrator snippets

We restrict this discussion to integrator snippet based algorithms, but more
general Markov snippet algorithms could be considered.

Consider again the target distribution

$\bar{\mu}({\rm d}z)=\sum_{k=0}^{T}\bar{\mu}(k,{\rm d}z)\,,$

with

$\bar{\mu}(k,{\rm d}z)=\frac{1}{T+1}\mu_{k}({\rm d}z)\,,$

where $\mu_{k}({\rm d}z):=\mu^{\psi^{-k}}({\rm d}z)$ for $k\in\llbracket
0,T\rrbracket$. Assume that we are in the context of Example 1, dropping $n$
for simplicity, the HMC algorithm using the integrator $\psi^{s}$ is as
follows

$\bar{P}(z,{\rm d}z^{\prime})=\alpha(z)\delta_{\psi^{s}(z)}({\rm
d}z)+\bar{\alpha}(z)\delta_{\sigma(z)}({\rm d}z^{\prime})$ (39)

with here

$\alpha(z)=\min\left\\{1,\frac{{\rm d}\bar{\mu}^{\sigma\circ\psi^{s}}}{{\rm
d}\bar{\mu}}(z)\right\\}=\min\left\\{1,\frac{\bar{\mu}\circ\psi^{w}(z)}{\bar{\mu}(z)}\right\\}$

where the last equality holds when $\upsilon\gg\bar{\mu}$,
$\upsilon^{\psi^{s}}=\upsilon$ and we let $\bar{\mu}(z)={\rm d}\bar{\mu}/{\rm
d}\upsilon(z)$. In other works the snippet
$\mathsf{z}=\big{(}z,\psi(z),\psi^{2}(z),\ldots,\psi^{T}(z)\big{)}$ is shifted
along the orbit $\\{\psi^{k}(z),k\in\mathbb{Z}\\}$ by $\psi^{s}$. Naturally
this needs to be combined with updates of the velocity to lead to a viable
ergodic MCMC algorithm. This can be achieved with the following kernel, with
$\psi^{k}(z)=\big{(}\psi_{x}^{k}(z),\psi_{v}^{k}(z)\big{)}$ and
$A\in\mathscr{Z}$,

$\bar{Q}(z,A):=\sum_{k,l=0}^{T}\bar{\mu}(k\mid
z)\frac{1}{T+1}\int\mathbf{1}\\{\psi^{-l}\big{(}\psi_{x}^{k}(z),v^{\prime}\big{)}\in
A\\}\varpi({\rm d}v^{\prime})\,,$

described algorithmically in Alg. 6. Indeed for any
$(l,v^{\prime})\in\llbracket 0,T\rrbracket\times\mathsf{V}$, using identity
(11),

$\int\sum_{k=0}^{T}\bar{\mu}(k\mid
z)\int\mathbf{1}\\{\psi^{-l}\big{(}\psi_{x}^{k}(z),v^{\prime}\big{)}\in
A\\}\bar{\mu}({\rm
d}z)=\int\sum_{k=0}^{T}\int\mathbf{1}\\{\psi^{-l}\big{(}x,v^{\prime}\big{)}\in
A\\}\mu({\rm d}z)\,,$

and hence

$\displaystyle\bar{\mu}\bar{Q}(A)$
$\displaystyle=\frac{1}{T+1}\sum_{l=0}^{T}\int\mathbf{1}\\{\psi^{-l}\big{(}x,v^{\prime}\big{)}\in
A\\}\varpi({\rm d}v^{\prime})\mu({\rm d}z)$
$\displaystyle=\frac{1}{T+1}\sum_{l=0}^{T}\int\mathbf{1}\\{\psi^{-l}\big{(}x,v\big{)}\in
A\\}\mu({\rm d}z)$ $\displaystyle=\int\mathbf{1}\\{z\in
A\\}\frac{1}{T+1}\sum_{l=0}^{T}\mu^{\psi^{-l}}({\rm d}z)$
$\displaystyle=\bar{\mu}(A)\,.$

Again samples from this MCMC targetting $\bar{\mu}$ can be used to estimate
expectations with respect to $\mu$ using the identity (11). This is closely
related to the “windows of states” approach of [28, 29, 31], where a window of
states is what we call a Hamiltonian snippet in the present manuscript. Indeed
the windows of states approach corresponds to the Markov update

$P(z,A)=\sum_{k,l=0}^{T}\frac{1}{T}\int\bar{P}\big{(}\psi^{-k}(z),{\rm
d}z^{\prime}\big{)}\bar{\mu}(l\mid
z^{\prime})\mathbf{1}\\{\psi^{l}(z^{\prime})\in A\\}\,,$ (40)

which, we show below, leaves $\mu$ invariant. Indeed, note that for
$k,l\in\llbracket 0,T\rrbracket$ and $A\in\mathscr{Z}$,

$\displaystyle\int\mu({\rm d}z)\int\bar{P}\big{(}\psi^{-k}(z),{\rm
d}z^{\prime}\big{)}\bar{\mu}(l\mid
z^{\prime})\mathbf{1}\\{\psi^{l}(z^{\prime})\in A\\}$
$\displaystyle=\int_{A}\mu^{\psi^{-k}}({\rm d}z)\int\bar{P}\big{(}z,{\rm
d}z^{\prime}\big{)}\bar{\mu}(l\mid
z^{\prime})\mathbf{1}\\{\psi^{l}(z^{\prime})\in A\\}\,,$

and therefore

$\displaystyle\int\mu({\rm d}z)\int P(z,{\rm
d}z^{\prime})\mathbf{1}\\{z^{\prime}\in A\\}=$
$\displaystyle\int\bar{\mu}({\rm d}z)\bar{P}(z,{\rm
d}z^{\prime})\sum_{l=0}^{T}\bar{\mu}(l\mid
z^{\prime})\mathbf{1}\\{\psi^{l}(z^{\prime})\in A\\}$ $\displaystyle=$
$\displaystyle\int\bar{\mu}({\rm d}z^{\prime})\sum_{l=0}^{T}\bar{\mu}(l\mid
z^{\prime})\mathbf{1}\\{\psi^{l}(z^{\prime})\in A\\}$ $\displaystyle=$
$\displaystyle\int\mu({\rm d}z)\mathbf{1}\\{z\in A\\}\,,$

where we obtain the last line from (11). $P$ is not $\mu$-reversible in
general, making theoretical comparisons challenging.

1Given $z$

2Sample $k\sim\bar{\mu}(\cdot\mid z)$ and $l\sim\mathcal{U}(\llbracket
0,T\rrbracket)$.

3Compute $\psi^{k}(z)=\big{(}\psi_{x}^{k}(z),\psi_{v}^{k}(z)\big{)}$.

4Sample $v^{\prime}\sim\varpi(\cdot)$

5Return $\psi^{-l}(\psi_{x}^{k}(z),v^{\prime})$

Algorithm 6 Kernel $\bar{Q}$

##  6 Discussion

We have shown how mappings used in various Monte Carlo schemes relying on
numerical integrators of ODE can be implemented to fully exploit all
computations to design robust and efficient sampling algorithms. Numerous
questions remain open, including the tradeoff between $N$ and $T$. A precise
analysis of this question is made particularly difficult by the fact that
integration along snippets is straightforwardly parallelizable, while
resampling does not lend itself to straightforward parallelisation.

Another point is concerned with the particular choice of mutation Markov
kernel $\bar{M}_{n}$, or $\bar{M}$, in (12) or (23). Indeed such a kernel
starts with a transition from samples approximating the snippet distribution
$\bar{\mu}_{n-1}$ to $\mu_{n-1}$, which is then followed by a reweighting of
samples leading to a representation of $\bar{\mu}_{n}$. Instead, for
illustration, one could suggest using an SMC sampler with (39) as mutation
kernel.

In relation to the discussion in Remark 20, a natural question is how our
scheme would compare with a “Rao-Blackwellized” SMC where weights of the type
(35), derived from (6) are used.

We leave all these questions for future investigations.

## Acknowledgements

The authors would like to thank Carl Dettman for very useful discussions on
Boltzman’s conjecture. Research of CA and MCE supported by EPSRC grant
‘CoSInES (COmputational Statistical INference for Engineering and Security)’
(EP/R034710/1), and EPSRC grant Bayes4Health, ‘New Approaches to Bayesian Data
Science: Tackling Challenges from the Health Sciences’ (EP/R018561/1).
Research of CZ was supported by a CSC Scholarship.

##  Appendix A Notation and definitions

We will write $\mathbb{N}=\left\\{0,1,2,\dots\right\\}$ for the set of natural
numbers and $\mathbb{R}_{+}=\left(0,\infty\right)$ for positive real numbers.
Throughout this section $(\mathsf{E},\mathscr{E})$ is a generic measurable
space.

  * •

For $A\subset\mathsf{E}$ we let $A^{\complement}$ be its complement.

  * •

$\mathcal{M}(\mathsf{E},\mathscr{E})$ (resp.
$\mathcal{P}(\mathsf{E},\mathscr{E})$ is the set of measures (resp.
probability distributions) on $(\mathsf{E},\mathscr{E})$

  * •

For a set $A\in\mathscr{E}$, its complement in $\mathsf{E}$ is denoted by
$A^{\complement}$. We denote the corresponding indicator function by
$\mathbf{1}_{A}:\mathsf{E}\to\left\\{0,1\right\\}$ and may use the notation
$\mathbf{1}\\{z\in A\\}:=\mathbf{1}_{A}(z)$.

  * •

For $\mu$ a probability measure on $(\mathsf{E},\mathscr{E})$ and a measurable
function $f\colon\mathsf{E}\rightarrow\mathbb{R}$ and , we let $\mu(f):=\int
f(x)\mu({\rm d}x)$.

  * •

For two probability measures $\mu$ and $\nu$ on $(\mathsf{E},\mathscr{E})$ we
let $\mu\otimes\nu$ be a measure on
$(\mathsf{E}\times\mathsf{E},\mathscr{E}\otimes\mathscr{E})$ such that
$\mu\otimes\nu(A\times B)=\mu(A)\nu(B)$ for $A,B\in\mathscr{E}$.

  * •

For a Markov kernel $P(x,{\rm d}y)$ on $\mathsf{E}\times\mathscr{E}$, we write

    * –

$\mu\otimes P$ for the probability measure on
$(\mathsf{E}\times\mathsf{E},\mathscr{E}\otimes\mathscr{E})$ such that for
$\bar{A}\in\mathscr{E}\otimes\mathscr{E}$, the minimal product
$\sigma$-algebra, $\mu\otimes P(\bar{A})=\int_{\bar{A}}\mu({\rm d}x)P(x,{\rm
d}y)$.

    * –

$\mu\accentset{\curvearrowleft}{\otimes}P$ for the probability measure on
$(\mathsf{E}\times\mathsf{E},\mathscr{E}\otimes\mathscr{E})$ such that for
$A,B\in\mathscr{E}$ $\mu\accentset{\curvearrowleft}{\otimes}P(A\times
B)=\mu\otimes P(B\times A)$.

  * •

For $\mu,\nu$ probability distributions on $(\mathsf{E},\mathscr{E})$ and
kernels $M,L\colon\mathsf{E}\times\mathscr{E\rightarrow}[0,1]$ such that
$\mu\otimes M\gg\nu\accentset{\curvearrowleft}{\otimes}L$ then we denote

$\frac{{\rm d}\nu\accentset{\curvearrowleft}{\otimes}L}{{\rm d}\mu\otimes
M}(z,z^{\prime})$

the corresponding Radon-Nikodym derivative such that for
$f\colon\mathsf{E}\times\mathsf{E}\rightarrow\mathbb{R}$,

$\int f(z,z^{\prime})\frac{{\rm
d}\nu\accentset{\curvearrowleft}{\otimes}L}{{\rm d}\mu\otimes
M}(z,z^{\prime})\mu\otimes M\big{(}{\rm d}(z,z^{\prime})\big{)}=\int
f(z,z^{\prime})\nu\accentset{\curvearrowleft}{\otimes}L\big{(}{\rm
d}(z^{\prime},z)\big{)}\,.$

  * •

A point mass distribution at $x$ will be denoted by $\delta_{x}({\rm d}y)$; it
is such that for $f\colon\mathsf{E}\rightarrow\mathbb{R}$

$\int f(x)\delta_{x}({\rm d}y)=f(x)$

  * •

In order to alleviate notation, for $M\in\mathbb{N}$,
$(z^{(i)},w_{i})\in\mathsf{E}\times[0,\infty)$, $i\in\llbracket M\rrbracket$,
we refer to $\big{\\{}(z^{(i)},w_{i}),i\in\llbracket M\rrbracket\big{\\}}$ as
weighted samples to mean $\big{\\{}(z^{(i)},\tilde{w}_{i}),i\in\llbracket
M\rrbracket\big{\\}}$ where $\tilde{w}_{i}\propto w_{i}$ but
$\sum_{i=1}^{M}\tilde{w}_{i}=1$.

  * •

We say that a set of weighted samples, or particles,
$\\{(z_{i},w_{i})\in\mathsf{Z}\times\mathbb{R}_{+}\colon i\in\llbracket
N\rrbracket\\}$ for $N\geq 1$ represents a distribution $\mu$ whenever for $f$
$\mu$-integrable

$\sum_{i=1}^{N}\frac{w_{i}}{\sum_{j=1}^{N}w_{j}}f(z_{i})\approx\mu(f)\,,$

in either in the $L^{p}$ sense for some $p\geq 1$.

  * •

For $M\in\mathbb{N}$, $w_{i}\in[0,\infty)$, $i\in\llbracket M\rrbracket$, we
let $K\sim{\rm Cat}\left(w_{1},w_{2},\ldots,w_{M}\right)$ mean that
$\mathbb{P}(K=k)\propto w_{k}$.

  * •

For $M,N\in\mathbb{N}$, $w_{ij}\in[0,\infty)$, $i\in\llbracket
M\rrbracket\times\llbracket N\rrbracket$, we let $K\sim{\rm
Cat}\left(w_{ij},i\in\llbracket M\rrbracket\times\llbracket
N\rrbracket\right)$ mean that $\mathbb{P}(K=(k,l))\propto w_{kl}$.

  * •

for $f\colon\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$ we let

    * –

$\nabla\otimes f$ be the transpose of the Jacobian

    * –

for $n=1$ we let $\nabla f=(\nabla\otimes f)^{\top}$ be the gradient,

    * –

$\nabla\cdot f$ be the divergence.

##  Appendix B Radon-Nikodym derivative

The general formalism required and used throughout the paper relies on a
unique measure theoretic tool, the Radon-Nikodym derivative. We gather here
definitions and intermediate results used throughout, pointing out the
simplicity of the tools involved and and the benefits they bring.

######  Definition 28 (Pushforward).

Let $\mu$ be a measure on $(\mathsf{Z},\mathscr{Z})$ and
$\psi:(\mathsf{Z},\mathscr{Z})\to(\mathsf{Z}^{\prime},\mathscr{F}^{\prime})$ a
measurable function. The pushforward of $\mu$ by $\psi$ is defined by

$\mu^{\psi}(A):=\mu(\psi^{-1}(A)),\qquad A\in\mathscr{Z}^{\prime},$

where $\psi^{-1}(A)=\\{z\in\mathsf{Z}:\psi(z)\in A\\}$ is the preimage of $A$
under $\psi$.

If $\mu$ is a probability distribution then $\mu^{\psi}$ is the probability
measure associated with $\psi(Z)$ when $Z\sim\mu$.

######  Definition 29 (Dominating and equivalent measures).

For two measures $\mu$ and $\nu$ on the same measurable space
$(\mathsf{Z},\mathscr{Z})$,

  1. 1.

$\nu$ is said to dominate $\mu$ if for all measurable $A\in\mathscr{Z}$,
$\mu(A)>0\Rightarrow\nu(A)>0$ – this is denoted $\nu\gg\mu$.

  2. 2.

$\mu$ and $\nu$ are equivalent, written $\mu\equiv\nu$, if $\mu\gg\nu$ and
$\nu\gg\mu$.

We will need the notion of Radon-Nikodym derivative [7, Theorems 32.2 &
16.11]:

######  Theorem 30 (Radon–Nikodym).

Let $\mu$ and $\nu$ be $\sigma$-finite measures on $(\mathsf{Z},\mathscr{Z})$.
Then $\nu\ll\mu$ if and only if there exists an essentially unique,
measurable, non-negative function $f\colon\mathsf{Z}\rightarrow[0,\infty)$
such that

$\int_{A}f(z)\mu({\rm d}z)=\nu(A),\qquad A\in\mathscr{E}.$

Therefore we can view ${\rm d}\nu/{\rm d}\mu:=f$ as the density of $\nu$ w.r.t
$\mu$ and in particular if $g$ is integrable w.r.t. $\nu$ then

$\int g(z)\frac{{\rm d}\nu}{{\rm d}\mu}(z)\mu({\rm d}z)=\int g(z)\nu({\rm
d}z)\,.$

If $\mu$ is a measure and $f$ a non-negative, measurable function then
$\mu\cdot f$ is the measure $(\mu\cdot f)(A)=\int{\bf 1}_{A}(z)f(z)\mu({\rm
d}z)$, i.e. the measure $\nu=\mu\cdot f$ such that $f$ is the Radon–Nikodym
derivative ${\rm d}\nu/{\rm d}\mu=f$.

The following establishes the expression of an expectation with respect to the
pushforward $\mu^{\psi}$ in terms of expectations with respect to $\mu$ [7,
Theorem 16.13].

######  Theorem 31 (Change of variables).

A function $f:\mathsf{Z}^{\prime}\to\mathbb{R}$ is integrable w.r.t.
$\mu^{\psi}$ if and only if $f\circ\psi$ is integrable w.r.t. $\mu$, in which
case

$\int_{\mathsf{Z}^{\prime}}f(z)\mu^{\psi}({\rm
d}z)=\int_{\mathsf{Z}}f\circ\psi(z)\mu({\rm d}z)\,.$ (41)

We now establish results useful throughout the manuscript. The central
identity used throughout the manuscript is a direct application of Theorem 31
for $\psi\colon\mathsf{Z}\rightarrow\mathsf{Z}$ invertible

$\int_{\mathsf{Z}^{\prime}}f\circ\psi(z)\mu^{\psi^{-1}}({\rm
d}z)=\int_{\mathsf{Z}}f(z)\mu({\rm d}z)$

which seems tautological since it can be summarized as follows: for
$Z\sim\mu$, then $\psi^{-1}(Z)\sim\mu^{\psi^{-1}}$ and
$\psi\circ\psi^{-1}(Z)\sim\mu$! However the interest of the approach stems
from the following properties.

######  Lemma 32.

Let $\psi\colon\mathsf{Z}\rightarrow\mathsf{Z}$ be measurable and integrable,
$\mu$ and $\upsilon$ be $\sigma$-finite measures on $(\mathsf{Z},\mathscr{Z})$
such that $\upsilon\gg\mu$ and $\upsilon\gg\upsilon^{\psi^{-1}}$. Then

  1. 1.

$\upsilon^{\psi^{-1}}\gg\mu^{\psi^{-1}}$ and therefore
$\upsilon\gg\mu^{\psi^{-1}}$,

  2. 2.

for $\upsilon$-almost all $z\in\mathsf{Z}$,

$\frac{{\rm d}\mu^{\psi^{-1}}}{{\rm d}\upsilon}(z)=\frac{{\rm d}\mu}{{\rm
d}\upsilon}\circ\psi(z)\frac{{\rm d}\upsilon^{\psi^{-1}}}{{\rm d}\upsilon}$

  3. 3.

we have

$\mu\gg\mu^{\psi^{-1}}\iff\upsilon\big{(}\big{\\{}z\in\mathsf{Z}\colon{\rm
d}\mu^{\psi^{-1}}/{\rm d}\upsilon(z)>0,{\rm d}\mu/{\rm
d}\upsilon(z)=0\big{\\}}\big{)}=0$

in which case for $\upsilon$-almost all $z\in\mathsf{Z}$

$\frac{{\rm d}\mu^{\psi^{-1}}}{{\rm
d}\mu}(z)=\begin{cases}\frac{\mu\circ\psi}{\mu}(z)&\mu(z)>0\\\
0&\text{otherwise}\end{cases}$

and therefore

$\int_{\mathsf{Z}^{\prime}}f\circ\psi(z)\mu^{\psi^{-1}}({\rm
d}z)=\int_{\mathsf{Z}^{\prime}}f\circ\psi(z)\frac{\mu\circ\psi(z)}{\mu(z)}\frac{{\rm
d}\upsilon^{\psi}}{{\rm d}\upsilon}\mu({\rm d}z).$

###### Proof.

For the first part of the first statement, let $A\in\mathscr{Z}$ such that
$\upsilon^{\psi^{-1}}(A)=\upsilon\big{(}\psi(A)\big{)}>0$, then since
$\psi(A)\in\mathscr{Z}$ and $\upsilon\gg\mu$ we deduce
$\mu\big{(}\psi(A)\big{)}=\mu^{\psi^{-1}}(A)>0$ and we conclude; the second
parts follows from $\upsilon\gg\upsilon^{\psi^{-1}}$. For the second statement
for $f\colon\mathsf{Z}\rightarrow\mathbb{R}$ bounded and measurable,

$\displaystyle\int f(z)\frac{{\rm d}\mu^{\psi^{-1}}}{{\rm
d}\upsilon}(z)\upsilon({\rm d}z)$ $\displaystyle=\int f(z)\mu^{\psi^{-1}}({\rm
d}z)$ $\displaystyle=\int f\circ\psi^{-1}(z)\mu({\rm d}z)$ $\displaystyle=\int
f\circ\psi^{-1}(z)\frac{{\rm d}\mu}{{\rm d}\upsilon}(z)\upsilon({\rm d}z)$
$\displaystyle=\int f(z)\frac{{\rm d}\mu}{{\rm
d}\upsilon}\circ\psi(z)\upsilon^{\psi^{-1}}({\rm d}z)$ $\displaystyle=\int
f(z)\frac{{\rm d}\mu}{{\rm d}\upsilon}\circ\psi(z)\frac{{\rm
d}\upsilon^{\psi^{-1}}}{{\rm d}\upsilon}(z)\upsilon({\rm d}z)$

The third statement is given as [7, Problem 32.6.], which we solve in Lemma
34. ∎

######  Corollary 33.

In the scenario when $\psi\colon\mathsf{Z}\rightarrow\mathsf{Z}$ and
$\upsilon$ are such that $\upsilon^{\psi}=\upsilon$ then ${\rm
d}\upsilon^{\psi}/{\rm d}\upsilon\equiv 1$.

######  Lemma 34 (Billingsley, Problem 32.6.).

Assume $\mu,\nu$ and $\upsilon$ are $\sigma$\- finite and that
$\upsilon\gg\nu,\mu$. Then $\mu\gg\nu$ if and only if
$\upsilon\big{(}\big{\\{}z\in\mathsf{Z}\colon{\rm d}\nu/{\rm
d}\upsilon(z)>0,{\rm d}\mu/{\rm d}\upsilon(z)=0\big{\\}}\big{)}=0$, in which
case

$\frac{{\rm d}\nu}{{\rm d}\mu}(z)=\mathbf{1}\\{z\in\mathsf{Z}:{\rm d}\mu/{\rm
d}\upsilon(z)>0\\}\frac{{\rm d}\nu}{{\rm d}\upsilon}/\frac{{\rm d}\mu}{{\rm
d}\upsilon}(z)\,.$

###### Proof.

Let $S:=\big{\\{}z\in\mathsf{Z}\colon{\rm d}\nu/{\rm d}\upsilon(z)>0,{\rm
d}\mu/{\rm d}\upsilon(z)=0\big{\\}}$. For
$f\colon\mathsf{Z}\rightarrow\mathbb{R}$ integrable we always have

$\displaystyle\int f(z)\nu({\rm d}z)$ $\displaystyle=\int\mathbf{1}\\{z\in
S\\}f(z)\frac{{\rm d}\nu}{{\rm d}\upsilon}(z)\upsilon({\rm
d}z)+\int\mathbf{1}\big{\\{}z\in S^{\complement}\big{\\}}f(z)\frac{{\rm
d}\nu}{{\rm d}\upsilon}/\frac{{\rm d}\mu}{{\rm d}\upsilon}(z)\mu({\rm d}z)\,.$

Assume $\mu\gg\nu$ then from above for any
$f\colon\mathsf{Z}\rightarrow\mathbb{R}$ integrable

$\displaystyle\int\mathbf{1}\\{z\in S\\}f(z)\frac{{\rm d}\nu}{{\rm
d}\mu}(z)\mu({\rm d}z)=$ $\displaystyle\int\mathbf{1}\\{z\in
S\\}f(z)\frac{{\rm d}\nu}{{\rm d}\mu}(z)\frac{{\rm d}\mu}{{\rm
d}\upsilon}(z)\upsilon({\rm d}z)$ $\displaystyle=$ $\displaystyle 0\,,$

and therefore $\upsilon\big{(}S\big{)}=0$ and we conclude from the first
identity above. Now assume that $\upsilon\big{(}S\big{)}=0$, then

$\displaystyle\int f(z)\nu({\rm d}z)$ $\displaystyle=\int\mathbf{1}\\{z\in
S\\}f(z)\frac{{\rm d}\nu}{{\rm d}\upsilon}(z)\upsilon({\rm
d}z)+\int\mathbf{1}\big{\\{}z\in S^{\complement}\big{\\}}f(z)\frac{{\rm
d}\nu}{{\rm d}\upsilon}/\frac{{\rm d}\mu}{{\rm d}\upsilon}(z)\mu({\rm d}z)$
$\displaystyle=\int\mathbf{1}\big{\\{}z\in
S^{\complement}\big{\\}}f(z)\frac{{\rm d}\nu}{{\rm d}\upsilon}/\frac{{\rm
d}\mu}{{\rm d}\upsilon}(z)\mu({\rm d}z)\,.$

The equivalence is therefore established and when either conditions is
satisfied we have

$\frac{{\rm d}\nu}{{\rm d}\mu}(z)=\mathbf{1}\big{\\{}z\in
S^{\complement}\big{\\}}\frac{{\rm d}\nu}{{\rm d}\upsilon}/\frac{{\rm
d}\mu}{{\rm d}\upsilon}(z)$

and we conclude. ∎

## References

  * [1] Christophe Andrieu, Anthony Lee, and Sam Livingstone.  A general perspective on the Metropolis-Hastings kernel.  arXiv preprint arXiv:2012.14881, 2020.
  * [2] Christophe Andrieu, James Ridgway, and Nick Whiteley.  Sampling normalizing constants in high dimensions using inhomogeneous diffusions, 2018.
  * [3] Mark A Beaumont.  Approximate bayesian computation in evolution and ecology.  Annual review of ecology, evolution, and systematics, 41:379–406, 2010.
  * [4] Alexandros Beskos, Natesh Pillai, Gareth Roberts, Jesus-Maria Sanz-Serna, and Andrew Stuart.  Optimal tuning of the hybrid Monte Carlo algorithm.  Bernoulli, 19(5A):1501–1534, 2013.
  * [5] Michael Betancourt.  Nested sampling with constrained Hamiltonian Monte Carlo.  In AIP Conference Proceedings, volume 1305, pages 165–172. American Institute of Physics, 2011.
  * [6] Joris Bierkens and Gareth Roberts.  A piecewise deterministic scaling limit of lifted Metropolis–Hastings in the Curie–Weiss model.  The Annals of Applied Probability, 27(2):846–882, 2017.
  * [7] P Billingsley.  Probability and measure. 3rd wiley.  New York, 1995.
  * [8] Alexandre Bouchard-Côté, Sebastian J Vollmer, and Arnaud Doucet.  The bouncy particle sampler: A nonreversible rejection-free Markov chain monte carlo method.  Journal of the American Statistical Association, 113(522):855–867, 2018.
  * [9] Dmitri Burago, Yuri Burago, and Sergei Ivanov.  A course in metric geometry, volume 33.  American Mathematical Society, 2022.
  * [10] Mari Paz Calvo, Daniel Sanz-Alonso, and Jesús María Sanz-Serna.  HMC: reducing the number of rejections by not using leapfrog and some results on the acceptance rate.  Journal of Computational Physics, 437:110333, 2021.
  * [11] Cédric M Campos and Jesús María Sanz-Serna.  Extra chance generalized hybrid Monte Carlo.  Journal of Computational Physics, 281:365–374, 2015.
  * [12] Joseph T Chang and David Pollard.  Conditioning as disintegration.  Statistica Neerlandica, 51(3):287–317, 1997.
  * [13] Nicolas Chopin and James Ridgway.  Leave pima indians alone: binary regression as a benchmark for Bayesian computation.  Statistical Science, 32(1):64–87, 2017.
  * [14] Hai-Dang Dau and Nicolas Chopin.  Waste-free sequential Monte Carlo.  arXiv preprint arXiv:2011.02328, 2020.
  * [15] Pierre Del Moral and Pierre Del Moral.  Feynman-kac formulae.  Springer, 2004.
  * [16] Pierre Del Moral, Arnaud Doucet, and Ajay Jasra.  Sequential Monte Carlo samplers.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(3):411–436, 2006.
  * [17] Simon Duane, Anthony D Kennedy, Brian J Pendleton, and Duncan Roweth.  Hybrid Monte Carlo.  Physics letters B, 195(2):216–222, 1987.
  * [18] Mauro Camara Escudero.  Approximate Manifold Sampling: Robust Bayesian Inference for Machine Learning.  PhD thesis, School of Mathematics, January 2024.
  * [19] Youhan Fang, J. M. Sanz-Serna, and Robert D. Skeel.  Compressible generalized hybrid Monte Carlo.  The Journal of Chemical Physics, 140(17):174108, 05 2014.
  * [20] Ernst Hairer, SP Nörsett, and G. Wanner.  Solving ordinary differential equations I, Nonstiff problems.  Springer-Verlag, 1993.
  * [21] Matthew Hoffman, Alexey Radul, and Pavel Sountsov.  An adaptive-MCMC scheme for setting trajectory lengths in Hamiltonian Monte Carlo.  In International Conference on Artificial Intelligence and Statistics, pages 3907–3915. PMLR, 2021.
  * [22] Matthew D Hoffman and Pavel Sountsov.  Tuning-Free Generalized Hamiltonian Monte Carlo.  In International Conference on Artificial Intelligence and Statistics, pages 7799–7813. PMLR, 2022.
  * [23] Paul B Mackenze.  An improved hybrid Monte Carlo method.  Physics Letters B, 226(3-4):369–371, 1989.
  * [24] Florian Maire and Pierre Vandekerkhove.  On markov chain monte carlo for sparse and filamentary distributions.  ArXiv e-prints, 2018.
  * [25] Florian Maire and Pierre Vandekerkhove.  Markov kernels local aggregation for noise vanishing distribution sampling.  SIAM Journal on Mathematics of Data Science, 4(4):1293–1319, 2022\.
  * [26] Youssef Marzouk, Tarek Moselhy, Matthew Parno, and Alessio Spantini.  Sampling via Measure Transport: An Introduction, pages 1–41.  Springer International Publishing, Cham, 2016.
  * [27] Xiao-Li Meng and Stephen Schilling.  Warp bridge sampling.  Journal of Computational and Graphical Statistics, 11(3):552–586, 2002.
  * [28] Radford M Neal.  An improved acceptance procedure for the hybrid Monte Carlo algorithm.  Journal of Computational Physics, 111(1):194–203, 1994.
  * [29] Radford M. Neal.  MCMC Using Hamiltonian Dynamics, chapter 5.  CRC Press, 2011.
  * [30] Ari Pakman and Liam Paninski.  Exact Hamiltonian Monte Carlo for truncated multivariate Gaussians.  Journal of Computational and Graphical Statistics, 23(2):518–542, 2014.
  * [31] Zhaohui S. Qin and Jun S. Liu.  Multipoint Metropolis method with application to Hybrid Monte Carlo.  Journal of Computational Physics, 172(2):827–840, 2001.
  * [32] James Ridgway.  Computation of gaussian orthant probabilities in high dimension.  Statistics and computing, 26(4):899–916, 2016.
  * [33] Grant M Rotskoff and Eric Vanden-Eijnden.  Dynamical computation of the density of states and Bayes factors using nonequilibrium importance sampling.  Physical review letters, 122(15):150602, 2019.
  * [34] Gabriel Stoltz, Mathias Rousset, et al.  Free energy computations: A mathematical perspective.  World Scientific, 2010.
  * [35] D Szasz.  Boltzmann’s ergodic hypothesis, a conjecture for centuries?, chapter Hard ball systems and the Lorentz gas, pages 421–446.  Springer, 2000.
  * [36] Esteban G Tabak and Eric Vanden-Eijnden.  Density estimation by dual ascent of the log-likelihood.  Communications in Mathematical Sciences, 8(1):217–233, 2010.
  * [37] Erik H Thiede, Brian Van Koten, Jonathan Weare, and Aaron R Dinner.  Eigenvector method for umbrella sampling enables error analysis.  The Journal of chemical physics, 145(8), 2016.
  * [38] Achille Thin, Yazid Janati El Idrissi, Sylvain Le Corff, Charles Ollion, Eric Moulines, Arnaud Doucet, Alain Durmus, and Christian X Robert.  Neo: Non equilibrium sampling on the orbits of a deterministic transform.  Advances in Neural Information Processing Systems, 34:17060–17071, 2021.
  * [39] G.M. Torrie and J.P. Valleau.  Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling.  Journal of Computational Physics, 23(2):187–199, 1977.
  * [40] Paul F Tupper.  Ergodicity and the numerical simulation of Hamiltonian systems.  SIAM Journal on Applied Dynamical Systems, 4(3):563–587, 2005.
  * [41] J. von Neumann.  Proof of the quasi-ergodic hypothesis.  Proc. Natl. Acad. Sci. USA, 18:70–82, 1932.
  * [42] Chang Zhang.  On the Improvements and Innovations of Monte Carlo Methods.  PhD thesis, School of Mathematics, https://research-information.bris.ac.uk/en/studentTheses/on-the-improvements-and-innovations-of-monte-carlo-methods, June 2022.


# Proton s-resonance states of 12C and 14,15O within the Skyrme Hartree-Fock
mean-field framework

Nguyen Le Anh<EMAIL_ADDRESS>Department of Theoretical Physics, Faculty
of Physics and Engineering Physics, University of Science, Ho Chi Minh City,
Vietnam  Vietnam National University, Ho Chi Minh City, Vietnam  Department of
Physics, Ho Chi Minh City University of Education, 280 An Duong Vuong,
District 5, Ho Chi Minh City, Vietnam  Young-ho Song<EMAIL_ADDRESS>Rare
Isotope Science Project, Institute for Basic Science, Daejeon 34047, Korea
Bui Minh Loc<EMAIL_ADDRESS>Center for Exotic Nuclear Studies,
Institute for Basic Science (IBS), Daejeon 34126, Korea

###### Abstract

The excitation functions of proton elastic scattering on 12C and 14,15O nuclei
at the energies near the proton-emission threshold are calculated using the
Skyrme Hartree-Fock (SHF) in continuum approach. For each excitation function,
the first resonance is identified as the $s$-state resonance of the mean-field
theory. For 15O, whose ground-state spin is nonzero, the $s$-state resonance
splits into two resonances via the spin-spin component of the optical
potential. With a slight adjustment of the strength of central potential,
which is obtained from the SHF in continuum approach, the excitation functions
of proton elastic scattering for the three nuclei can be explained with high
accuracy. The proposed framework can provide a practical method to explain
nuclear scattering at the energies near the proton-emission threshold with
minimal experimental input.

##  I INTRODUCTION

Unbound atomic nuclei, which can be experimentally observed as resonance peaks
in the nuclear excitation functions, have been widely studied through resonant
elastic scattering at low energy. Experiments of this kind are performed at
radioactive beam facilities, such as CERN-ISOLDE (Europe), GANIL (France), GSI
(Germany), RIKEN (Japan), and TRIUMF (Canada). The properties of these
resonances reveal the coupling of discrete states to a scattering continuum
and information about spectroscopic properties of unbound nuclei [1]. Besides,
the data collected from these experiments, particularly those of low-lying
resonances, are crucial for accurately computing cross sections and nuclear
reaction rates in astrophysics, such as radiative-capture and transfer
reactions [2].

Light nuclei, especially oxygen isotopes, merit further discussion. Both the
neutron-rich and proton-rich oxygen isotopes offer an important testing ground
for understanding the structure of exotic nuclei. On the proton-rich side, the
properties of the ground states and low-lying excited states of unbound nuclei
15,16F were measured and discussed from proton elastic scattering on proton-
rich isotopes 14,15O at sub-MeV energies [3, 4, 5, 6, 7], and recently [8, 9].
As a light nucleus at the proton drip line, 16F is an ideal case for studying
the structure and reaction of exotic nuclei. In addition, the structure of 16F
has been discussed in the context of isospin-symmetry breaking [10, 11, 12].

As resonances play an essential role not only in nuclear physics but also in
physics in general, theoretical analysis was developed back in the early day
of the field with the $R$-matrix theory [13, 14]. Nowadays, the properties of
resonances are widely analyzed with the $R$-matrix method, see Refs. [15, 16].
Besides, new theoretical models have been developed. The coupled cluster model
[17] is specialized for describing the low-lying resonances. The shell model
with coupling to the continuum known as the Gamow shell model [18, 19, 20] has
been applied with success to reproduce the experimental results including
nuclear states with a complicated configuration [12]. In the shell-model
approach, low-lying resonances near the emission threshold are suggested to be
closely related to the configuration of the core and the proton in a single-
particle state. In the framework of the mean-field theory for scattering
problems (i.e. the optical potential), they can be simply identified as
single-particle resonant states.

The optical-potential analysis for nucleon scattering from the target with
nonzero spin provides information on spin dependence of the nucleon-nucleus
optical potential such as the spin-spin component proposed long ago in Ref.
[21]. It is crucial in the low-energy elastic scattering on light nuclei. As
the absorption is absent at the energy of a few MeV, the elastic scattering
reflects the nuclear structure as a function of energy, and different spin
states are observed in the energy spectrum. A number of experiments were
designed to measure the properties of the spin-spin component [22, 23, 24, 25,
26, 27, 28]. The work in Ref. [29] discussed that uncertainties still remand
in experimental and theoretical analyses because of the experimental precision
and complicated nuclear mechanism. This subject was left aside for a while.
Nowadays, as radioactive beam facilities produce more exotic nuclei, a highly
accurate optical model potential including spin-dependent term is essential.
The microscopic nucleon optical model for target nuclei with nonzero spin
within the folding model was intensively re-examined [30]. Within the folding
calculation, the spin dependence of the nucleon-nucleus optical potential is
directly linked to that of the effective nucleon-nucleon interaction. The
recent study of spin-polarized neutron stars emphasized the necessity of
precise measurements for the strength of the spin dependence of effective
nucleon-nucleon interactions [31].

In our present work, proton elastic scattering was analyzed using the Skyrme
Hartree-Fock (SHF) optical potential, where the SHF calculation is extended
for the continuum. The same mean-field potential that originally describes the
single-particle bound-state also accounts for the single-particle scattering
states [32, 33, 34]. The method of calculation was applied successfully on
several low-energy radiative-capture reactions [35, 36, 37]. Recently, the new
low-lying narrow resonant state found in the 10Be($p,p$) scattering [38, 39]
was also explained as an $s$-state resonance [40].

The description of the SHF optical potentials for target nuclei with various
ground-state spins-parities is presented in the next section. Our present work
is restricted to the $s$-wave scattering. The contribution of the spin-orbit
potential vanishes which is convenient for the study of the spin-spin
potential. The method was applied on a stable nucleus 12C and proton-rich
oxygen isotopes 14,15O. The $s$-wave scattering data and resonances are well
reproduced with minimal adjustment of the SHF central potential strength. Our
result gives information about the properties of the nuclear central
potential, especially the spin-spin potential.

##  II Skyrme Hartree-Fock optical potential

As only the $s$-wave is of interest, the scattering equation is written as

$\displaystyle\left(-\dfrac{\hbar^{2}}{2m}\dfrac{d^{2}}{dr^{2}}+\lambda V_{\rm
cent}+V_{\rm Coul}\right)\chi_{s}=E\chi_{s}.$ (1)

The nuclear central part $V_{\rm cent}$ and the Coulomb part $V_{\rm Coul}$ of
the optical potential are obtained from the SHF calculation [33, 34]:

$\displaystyle V_{\rm cent}$ $\displaystyle=\dfrac{m^{*}}{m}\left[V^{\rm
HF}_{\rm
cent}+\dfrac{1}{2}\dfrac{d^{2}}{dr^{2}}\dfrac{\hbar^{2}}{2m^{*}}-\dfrac{m^{*}}{2\hbar^{2}}\left(\dfrac{d}{dr}\dfrac{\hbar^{2}}{2m^{*}}\right)^{2}\right]+\left(1-\dfrac{m^{*}}{m}\right)E,$
(2) $\displaystyle V_{\rm Coul}$ $\displaystyle=\dfrac{m^{*}}{m}V^{\rm
HF}_{\rm Coul}.$ (3)

The one-body potentials including $V^{\rm HF}_{\rm cent}$ and $V^{\rm HF}_{\rm
Coul}$, and the effective mass $m^{*}$ are given in the SHF formalism as
functions of $r$ only. One adjustable parameter $\lambda$ for the nuclear
central potential in Eq. (1) is introduced to reproduce the exact location of
the low-energy $s$-state resonance. Note that the optical potential is purely
real at this low-energy range. There is a slight energy dependence in $V_{\rm
cent}$ caused by the second term in Eq. (2).

In the case the target has a nonzero spin $I$ as for 15O($p,p$), the spin-spin
interaction plays a role. The optical potential for the target with nonzero
spin is

$\displaystyle V=\lambda V_{\rm cent}+V_{\rm
Coul}+V_{sI}({\bm{s}}\cdot{\bm{I}}).$ (4)

The last term in Eq. (4) is the spin-spin potential [21, 41]. To minimize the
number of parameters, the spin-spin potential has the same form factor as the
central term [42, 43], i.e. $V_{sI}=\alpha V_{\rm cent}$. The validity of this
assumption will be justified by the comparison with the experimental data.
Under this simple assumption, the potential is rewritten as

$\displaystyle V_{J}=\lambda_{J}V_{\rm cent}+V_{\rm Coul},$ (5)

where $\lambda_{J}\equiv\left[\lambda+\alpha({\bm{s}}\cdot{\bm{I}})\right]$.
In particular, the 15O($p,p$) reaction has parameters
$\lambda_{J=0}=\lambda-\frac{3}{4}\alpha$ for the $0^{-}$ resonance and
$\lambda_{J=1}=\lambda+\frac{1}{4}\alpha$ for the $1^{-}$ resonance in 16F.
The parameters $\lambda_{J=0,1}$ are adjusted according to the experimental
energies of the resonances. The final result is obtained by averaging the
parameters for the two resonance states with ($2J+1$) weight. All potentials
are deduced from the SHF calculation using the computer program skyrme_rpa
provided in Ref. [44]. The scattering equation is solved by the ECIS06 code
[45]. The single-particle width of the $s$-wave resonance is extracted from
the inverse of the first derivative of phase shifts $\delta(E)$ at the
resonance energy $E_{R}$ in the center of mass frame

$\Gamma_{\text{cal}}=2\left[\dfrac{d\delta(E)}{dE}\bigg{|}_{E=E_{R}}\right]^{-1}.$
(6)

##  III Results and discussions

###  III.1 Spinless targets 12C and 14O

First, the 12C($p,p$) scattering is considered, which exhibits the first
resonance at $420$ keV above the proton-emission threshold of 13N ($1.94$
MeV). It is in good agreement with the observed resonance with $\Gamma=31.7$
keV [46] at the same energy in the proton radiative-capture reaction. In the
mean-field theory, this resonance is identified as the $s$-state resonance
corresponding to the state $1/2^{+}$ of 13N at $2.36$ MeV. Fig. 1 shows our
results with and without the adjustment of the central potential from the SLy4
interaction [47]. The slight increase of $\lambda$ from $1.00$ to $1.03$ led
to a great agreement with the experimental data, especially the shape of the
resonance. The SLy4 interaction without the adjustment of $\lambda$ also
predicts the existence of the $s$-state resonance; however, the location is
not exact as shown in Fig. 1. For the purely single-particle state such as
this state, the SHF calculation reproduces the energy distribution of the
cross section that agrees with the experimental data in Ref. [48] and
phenomenological $R$-matrix [49]. The angular distribution of the cross
section is also obtained with the same potential as shown in Fig. 2. The
necessity of adjustment of the parameter $\lambda$ is significant near the
resonance. The calculated angular distribution of the cross section at
energies outside the resonance region does not change (Figs. 2(a) and 2(b))
while there is a remarkable difference in the vicinity of the resonance (Figs.
2(c) and 2(d)). Note that while the reported SHF calculation is based on Ref.
[44] which is assumed spherical symmetry and excluded pairing correlations,
the calculation without adjustment of the strength of the potential still
provides a good description of the nonresonant cross section and hints of
possible resonances in the energy range of interest for the investigated
nuclei.

Figure 1: The excitation function of the 12C($p,p$) scattering in the
$300$–$1000$ keV range. The resonance at $420$ keV is the $s$-state resonance.
The experimental data at the scattering angle of $170^{\circ}$ are taken from
Ref. [48]. Figure 2: Angular distributions of differential cross section of
12C($p,p$) around the resonance of $420$ keV. The slight adjustment of the
scaling parameter $\lambda$ from $1.00$ to $1.03$ (solid lines) gives the best
agreement with literature, especially at $453$ keV. The experimental data are
taken from Ref. [50].

The same method can be applied not only for stable nuclei but also for exotic
nuclei without additional modifications. The keV-proton elastic scattering
from 14O was first reported in Refs. [3, 4]. It was then re-measured with
higher precision [5, 8, 9]. The first two resonances in 15F were described as
the single-particle resonances using the Woods-Saxon potential in Ref. [5].
Our calculation confirms that the $1/2^{+}$ ground state of 15F is the result
of the single-particle resonances of the $s$-scattering wave. Fig. 3 presents
the calculated differential cross section of the 14O($p,p$) scattering using
the SHF optical potential with $\lambda=1.02$, which is consistent with the
two experimental data sets [8, 9].

Figure 3: Excitation function of 14O($p,p$) below $2$ MeV. The resonance at
$1.27$ MeV is the $s$-state resonance. The calculated different cross section
was obtained with $\lambda=1.02$. The experimental data at $\theta_{\rm
lab}=180^{\circ}$ are from Ref. [8] (circles) and Ref. [9] (squares).

###  III.2 Nonzero spin target 15O and the role of spin-spin potential

The unbound nucleus 16F was experimentally examined and discussed in Refs.
[51, 7, 8]. The first two states $0^{-}$ and $1^{-}$ are the $s$-resonance
states. To explain this splitting, the spin-spin interaction as described
earlier is included. Indeed, when taking into account the effect of the 15O
target spin of $1/2^{-}$, the $s$-wave scattering generates two resonances
with the total spins of $J^{\pi}=0^{-}$ and $J^{\pi}=1^{-}$. The
experimentally observed two resonances are well reproduced by SHF optical
potential using $\lambda_{J=0}=0.99$ for the first at $0.57$ MeV (dashed line)
and $\lambda_{J=1}=0.97$ for the second at $0.78$ MeV (dotted line) as shown
in Fig. 4. The values of the parameters show the adjustment of the central
potential is 0.98, while the parameter $\alpha$ which describes the relative
strength of the spin-spin potential compared to the central potential is
$\alpha=0.02$. Without any adjustment ($\lambda=1,\alpha=0$), the existence of
the resonances can still be predicted by SHF optical potential. Although the
strength of the spin-spin potential is small, it is responsible for the
splitting of the $0^{-}$ and $1^{-}$ resonances.

Figure 4: Excitation function of 15O($p,p$) below $1000$ keV with two
$s$-state resonances. $\lambda_{J=0}=0.99$ and $\lambda_{J=1}=0.97$ were
adopted for the calculation at the two resonances. The experimental data at
$\theta_{\rm lab}=180^{\circ}$ are from Ref. [8].

The agreement between the experimental data and SHF calculation points out
that the resonant states are in purely single-particle configurations. In many
other cases such as $p$-, $d$-states, and even the $s$-resonance state in 11B
in which the $\alpha$-cluster configuration is important [40], single-particle
resonances are mixed with other configurations. Therefore, the $s$-state
resonances in the study are the sensitive probes of the single-particle mean-
field potential. Furthermore, its splitting is invaluable for studying the
spin-spin component of the optical potential.

In the past, the test of the spin-spin potential required sensitive
measurements in which the beam was polarized in the depolarized measurements,
or the beam and the target were both polarized in the transmission
measurements. The theoretical formalism for the past experiments was
established using the distorted wave Born approximation [41], and the coupled-
channel formalism [52, 53]. Because of the complicated nuclear reaction
mechanism, theoretical calculations were significantly improved [22, 54, 55].
The microscopic approach within the folding model using the realistic nucleon-
nucleon interaction was applied in Refs. [56, 30]. It has been considered that
these studies provided a unique test of the spin dependence in the optical
potential.

In our work, the analysis of experimental data measured in Ref. [8] within the
SHF approach provides another reliable method to determine the sign,
magnitude, and shape of the spin-spin component of the optical potential. As
the elastic scattering was analyzed with the optical potential, the
complicated reaction mechanism is completely avoided. The result shows that
the strength of the spin-spin potential is 2% of the central term. It is
consistent with the previous conclusions that the strength of spin-spin
interaction is of the order of 1 MeV. Although the spin-spin interaction is
weak, it is important for understanding its shape. The assumption that the
spin-spin potential takes the same form as the central term following Ref.
[42] provides an excellent result in our work.

The calculated widths of resonances are presented in Table 1 which are not far
from the observed width. For the narrow resonance in the stable nucleus 12C,
the location is measured precisely at 420 keV and the width is $31.7\pm 0.8$
keV [46]. The width of 26 keV is computed from Eq. (6) at where the phase
shift is $\pi/2$. The unbound states of 16F are observed as narrow resonances.
The experimental proton-emission widths are given in Ref. [8] with high
accuracy. The calculated widths for $0^{-}$ and $1^{-}$ of 16F are smaller by
10 keV compared to the measurement. The results of narrow widths pointed out
that proton-decay channels are mainly dominated in the unbound states of 13N
and 15,16F. In the broad resonance of the unbound state of 15F, the analysis
of the data using the potential model gave the width of approximately 700 keV
with positions depending on the different definitions for the location of the
resonance [6]. The estimated ground state decay width for 15F was $800\pm 300$
keV presented in Refs. [3, 57]. The relation of the width and the diffuseness
parameter within the potential model was discussed carefully in Ref. [6]. Our
calculation shows the width is 534 or 773 keV corresponding to the definition
of resonance position, 1270 or 1400 keV.

Scattering | $J^{\pi}$ |  $E_{R}$ (c.m.) | $\lambda_{J}$ | $\Gamma_{\text{cal}}$ | $\Gamma_{\text{exp}}$
---|---|---|---|---|---
|  | [keV] |  | [keV] | [keV]
p + 12C | $1/2^{+}$ | $420$ | $1.03$ | $26$ |  $31.7\pm 0.8$ [46]
p + 14O | $1/2^{+}$ | $1270$ | $1.02$ | $534$ |  $\sim 700$ [6]
|  | $1400$ | $1.00$ | $773$ |  $800\pm 300$ [3, 57],
|  |  |  |  |  $>900$ [4]
p + 15O | $0^{-}$ | $536$ | $0.99$ | $15$ |  $25.6\pm 4.6$ [8]
| $1^{-}$ | $729$ | $0.97$ | $63$ |  $76\pm 5$ [8]
Table 1: Level parameters for the resonances produced in the Skyrme Hartree-
Fock approach. The differences between $\Gamma_{\text{cal}}$ and
$\Gamma_{\text{exp}}$ are discussed in the text.

Finally, considering the precision of the experimental data used in the
analysis, it is worth mentioning the distortion relates to the energy
resolution with the radioactive beam. It is about 20 keV in the laboratory
frame according to Refs. [9, 8]. The experimental energy resolution on our
results is negligible.

##  IV Conclusions and outlooks

The SHF formalism is a powerful tool for the investigation of not only nuclear
structure but also low-energy elastic scattering. With just one parameter, the
optical potentials for the elastic scattering of protons on stable and
unstable targets 12C and 14,15O are obtained from the standard SHF calculation
and applied successfully for the resonant elastic scattering in the study of
unbound nuclei. The $s$-state resonance and its splitting are sensitive probes
of the nuclear central and the spin-spin potential terms, respectively.

In order to improve the method of calculation, one possibility is to use a new
Skyrme interaction which is suitable for both bound and unbound nuclei and can
be used to study both nuclear structure and scattering. This could be used to
further understand how the central potentials of unbound nuclei differ from
those of bound nuclei. Additionally, it would be of interest to investigate
whether the relatively small magnitude of the spin-spin interaction compared
to the central potential, as seen in the 15O($p,p$) scattering, is a general
property for other nuclei. This would provide additional insights into the
role of spin-spin interaction in nuclei. In general, future efforts to refine
the calculation based on a more microscopic approach are ongoing.

## Acknowledgments

The authors would like to thank Prof. Naftali Auerbach for the fruitful
discussions, Dr. Laszlo Stuhl, and Dr. Jason Park for their discussions and
for carefully reading the manuscript. The work of B.M.L. was supported by the
Institute for Basic Science (IBS-R031-D1). The work of Y.H.S. was supported by
the Rare Isotope Science Project of Institute for Basic Science, funded by the
Ministry of Science and ICT (MSICT) and by the National Research Foundation of
Korea (2013M7A1A1075764).

## References

  * Johnson  _et al._ [2020] C. W. Johnson, K. D. Launey, N. Auerbach, S. Bacca, B. R. Barrett, C. R. Brune, M. A. Caprio, P. Descouvemont, W. H. Dickhoff, C. Elster, P. J. Fasano, K. Fossez, H. Hergert, M. Hjorth-Jensen, L. Hlophe, B. Hu, R. M. I. Betan, A. Idini, S. König, K. Kravvaris, D. Lee, J. Lei, A. Mercenne, R. N. Perez, W. Nazarewicz, F. M. Nunes, M. Płoszajczak, J. Rotureau, G. Rupak, A. M. Shirokov, I. Thompson, J. P. Vary, A. Volya, F. Xu, R. G. T. Zegers, V. Zelevinsky, and X. Zhang, J. Phys. G: Nucl. Part. Phys. 47, 123001 (2020).
  * Descouvemont [2020] P. Descouvemont, Front. Astron. Space Sci. 7, 9 (2020).
  * KeKelis  _et al._ [1978] G. J. KeKelis, M. S. Zisman, D. K. Scott, R. Jahn, D. J. Vieira, J. Cerny, and F. Ajzenberg-Selove, Phys. Rev. C 17, 1929 (1978).
  * Benenson  _et al._ [1978] W. Benenson, E. Kashy, A. G. Ledebuhr, R. C. Pardo, R. G. H. Robertson, and L. W. Robinson, Phys. Rev. C 17, 1939 (1978).
  * Peters  _et al._ [2003] W. A. Peters, T. Baumann, D. Bazin, B. A. Brown, R. R. C. Clement, N. Frank, P. Heckman, B. A. Luther, F. Nunes, J. Seitz, A. Stolz, M. Thoennessen, and E. Tryggestad, Phys. Rev. C 68, 034607 (2003).
  * Goldberg  _et al._ [2004] V. Z. Goldberg, G. G. Chubarian, G. Tabacaru, L. Trache, R. E. Tribble, A. Aprahamian, G. V. Rogachev, B. B. Skorodumov, and X. D. Tang, Phys. Rev. C 69, 031302 (2004).
  * Lee  _et al._ [2007] D. W. Lee, K. Peräjärvi, J. Powell, J. P. O’Neil, D. M. Moltz, V. Z. Goldberg, and J. Cerny, Phys. Rev. C 76, 024314 (2007).
  * Stefan  _et al._ [2014] I. Stefan, F. de Oliveira Santos, O. Sorlin, T. Davinson, M. Lewitowicz, G. Dumitru, J. C. Angélique, M. Angélique, E. Berthoumieux, C. Borcea, R. Borcea, A. Buta, J. M. Daugas, F. de Grancey, M. Fadil, S. Grévy, J. Kiener, A. Lefebvre-Schuhl, M. Lenhardt, J. Mrazek, F. Negoita, D. Pantelica, M. G. Pellegriti, L. Perrot, M. Ploszajczak, O. Roig, M. G. Saint Laurent, I. Ray, M. Stanoiu, C. Stodel, V. Tatischeff, and J. C. Thomas, Phys. Rev. C 90, 014307 (2014).
  * de Grancey  _et al._ [2016] F. de Grancey, A. Mercenne, F. de Oliveira Santos, T. Davinson, O. Sorlin, J. C. Angélique, M. Assié, E. Berthoumieux, R. Borcea, A. Buta, I. Celikovic, V. Chudoba, J. M. Daugas, G. Dumitru, M. Fadil, S. Grévy, J. Kiener, A. Lefebvre-Schuhl, N. Michel, J. Mrazek, F. Negoita, J. Okołowicz, D. Pantelica, M. Pellegriti, L. Perrot, M. Płoszajczak, G. Randisi, I. Ray, O. Roig, F. Rotaru, M. G. Saint Laurent, N. Smirnova, M. Stanoiu, I. Stefan, C. Stodel, K. Subotic, V. Tatischeff, J. C. Thomas, P. Ujić, and R. Wolski, Phys. Lett. B 758, 26 (2016).
  * Fortune [2006] H. T. Fortune, Phys. Rev. C 74, 054310 (2006).
  * Fortune [2018] H. T. Fortune, Phys. Rev. C 97, 044314 (2018).
  * Michel  _et al._ [2022] N. Michel, J. G. Li, L. H. Ru, and W. Zuo, Phys. Rev. C 106, L011301 (2022).
  * Wigner and Eisenbud [1947] E. P. Wigner and L. Eisenbud, Phys. Rev. 72, 29 (1947).
  * Lane and Thomas [1958] A. M. Lane and R. G. Thomas, Rev. Mod. Phys. 30, 257 (1958).
  * Descouvemont and Baye [2010] P. Descouvemont and D. Baye, Rep. Prog. Phys. 73, 036301 (2010).
  * Azuma  _et al._ [2010] R. E. Azuma, E. Uberseder, E. C. Simpson, C. R. Brune, H. Costantini, R. J. de Boer, J. Görres, M. Heil, P. J. LeBlanc, C. Ugalde, and M. Wiescher, Phys. Rev. C 81, 045805 (2010).
  * Baye  _et al._ [2005] D. Baye, P. Descouvemont, and F. Leo, Phys. Rev. C 72, 024309 (2005).
  * Hagen  _et al._ [2006] G. Hagen, M. Hjorth-Jensen, and N. Michel, Phys. Rev. C 73, 064307 (2006).
  * Tsukiyama  _et al._ [2009] K. Tsukiyama, M. Hjorth-Jensen, and G. Hagen, Phys. Rev. C 80, 051301 (2009).
  * Sun  _et al._ [2017] Z. H. Sun, Q. Wu, Z. H. Zhao, B. S. Hu, S. J. Dai, and F. R. Xu, Phys. Lett. B 769, 227 (2017).
  * Feshbach [1958] H. Feshbach, Annu. Rev. Nucl. Sci. 8, 49 (1958).
  * Nagamine  _et al._ [1970] K. Nagamine, A. Uchida, and S. Kobayashi, Nucl. Phys. A 145, 203 (1970).
  * Batty [1971] C. J. Batty, Nucl. Phys. A 178, 17 (1971).
  * Fisher  _et al._ [1972] T. R. Fisher, H. A. Grench, D. C. Healey, J. S. McCarthy, D. Parks, and R. Whitney, Nucl. Phys. A 179, 241 (1972).
  * Birchall  _et al._ [1974] J. Birchall, H. E. Conzett, J. Arvieux, W. Dahme, and R. M. Larimer, Phys. Lett. B 53, 165 (1974).
  * von Przewoski  _et al._ [1991] B. von Przewoski, P. D. Eversheim, F. Hinterberger, U. Lahr, J. Campbell, J. Götz, M. Hammans, R. Henneck, G. Masson, I. Sick, and W. Bauhoff, Nucl. Phys. A 528, 159 (1991).
  * Henneck  _et al._ [1994] R. Henneck, G. Masson, P. D. Eversheim, R. Gebel, F. Hinterberger, U. Lahr, H. W. Schmitt, J. Schleef, and B. von Przewoski, Nucl. Phys. A 571, 541 (1994).
  * Hannen  _et al._ [2003] V. M. Hannen, K. Amos, A. M. van den Berg, R. K. Bieber, P. K. Deb, F. Ellinghaus, D. Frekers, M. Hagemann, M. N. Harakeh, J. Heyse, M. A. de Huu, B. A. M. Krüsemann, S. Rakers, R. Schmidt, S. Y. van der Werf, and H. J. Wörtche, Phys. Rev. C 67, 054321 (2003).
  * Nagadi  _et al._ [2004] M. M. Nagadi, G. J. Weisel, R. L. Walter, J. P. Delaroche, and P. Romain, Phys. Rev. C 70, 064604 (2004).
  * Cunningham  _et al._ [2013] E. S. Cunningham, J. S. Al-Khalili, and R. C. Johnson, Phys. Rev. C 87, 054601 (2013).
  * Khoa  _et al._ [2022] N. H. D. Khoa, N. H. Tan, and D. T. Khoa, Phys. Rev. C 105, 065802 (2022).
  * Vautherin and Vénéroni [1968] D. Vautherin and M. Vénéroni, Phys. Lett. B 26, 552 (1968).
  * Dover and Van Giai [1971] C. B. Dover and N. Van Giai, Nucl. Phys. A 177, 559 (1971).
  * Dover and Van Giai [1972] C. B. Dover and N. Van Giai, Nucl. Phys. A 190, 373 (1972).
  * Le Anh and Minh Loc [2021] N. Le Anh and B. Minh Loc, Phys. Rev. C 103, 035812 (2021).
  * Le Anh  _et al._ [2021] N. Le Anh, P. Nhut Huan, and B. Minh Loc, Phys. Rev. C 104, 034622 (2021).
  * Le Anh and Minh Loc [2022] N. Le Anh and B. Minh Loc, Phys. Rev. C 106, 014605 (2022).
  * Ayyad  _et al._ [2022] Y. Ayyad, W. Mittig, T. Tang, B. Olaizola, G. Potel, N. Rijal, N. Watwood, H. Alvarez-Pol, D. Bazin, M. Caamaño, J. Chen, M. Cortesi, B. Fernández-Domínguez, S. Giraud, P. Gueye, S. Heinitz, R. Jain, B. P. Kay, E. A. Maugeri, B. Monteagudo, F. Ndayisabye, S. N. Paneru, J. Pereira, E. Rubino, C. Santamaria, D. Schumann, J. Surbrook, L. Wagner, J. C. Zamora, and V. Zelevinsky, Phys. Rev. Lett. 129, 012501 (2022).
  * Lopez-Saavedra  _et al._ [2022] E. Lopez-Saavedra, S. Almaraz-Calderon, B. W. Asher, L. T. Baby, N. Gerken, K. Hanselman, K. W. Kemper, A. N. Kuchera, A. B. Morelock, J. F. Perello, E. S. Temanson, A. Volya, and I. Wiedenhöver, Phys. Rev. Lett. 129, 012502 (2022).
  * Le Anh  _et al._ [2022] N. Le Anh, B. Minh Loc, N. Auerbach, and V. Zelevinsky, Phys. Rev. C 106, L051302 (2022).
  * Davies and Satchler [1964] K. T. R. Davies and G. R. Satchler, Nucl. Phys. 53, 1 (1964).
  * Amos  _et al._ [2003] K. Amos, L. Canton, G. Pisent, J. Svenne, and D. van der Knijff, Nucl. Phys. A 728, 65 (2003).
  * Amos  _et al._ [2021] K. Amos, P. R. Fraser, S. Karataglidis, and L. Canton, Eur. Phys. J. A 57, 165 (2021).
  * Colò  _et al._ [2013] G. Colò, L. Cao, N. Van Giai, and L. Capelli, Comput. Phys. Commun. 184, 142 (2013).
  * Raynal [1971] J. Raynal, ECIS06 code is distributed by the NEA DATA Bank, Paris, France; ECIS72 described in “Optical model and coupled-channels calculations in nuclear physics,”, in  _Computing as a language of physics. ICTP International Seminar Course, Trieste, Italy_ (1971) p. 281.
  * Ajzenberg-Selove [1991] F. Ajzenberg-Selove, Nucl. Phys. A 523, 1 (1991).
  * Chabanat  _et al._ [1998] E. Chabanat, P. Bonche, P. Haensel, J. Meyer, and R. Schaeffer, Nucl. Phys. A 635, 231 (1998).
  * Zhengmin  _et al._ [1993] L. Zhengmin, L. Beijing, D. Zhenzhong, and H. Huimin, Nucl. Instrum. Methods Phys. Res. B 74, 439 (1993).
  * Meyer [1976] H. O. Meyer, Z. Physik A 279, 41 (1976).
  * Mazzoni  _et al._ [1998] S. Mazzoni, M. Chiari, L. Giuntini, P. Mandò, and N. Taccetti, Nucl. Instrum. Methods Phys. Res. B 136-138, 86 (1998).
  * Fortune [1995] H. T. Fortune, Phys. Rev. C 52, 2261 (1995).
  * Tamura [1965] T. Tamura, Rev. Mod. Phys. 37, 679 (1965).
  * Stamp [1967] A. P. Stamp, Phys. Rev. 153, 1052 (1967).
  * Hussein and Sherif [1973] A. H. Hussein and H. S. Sherif, Phys. Rev. C 8, 518 (1973).
  * McAbee  _et al._ [1990] T. L. McAbee, W. J. Thompson, and H. Ohnishi, Nucl. Phys. A 509, 39 (1990).
  * Cunningham  _et al._ [2011] E. S. Cunningham, J. S. Al-Khalili, and R. C. Johnson, Phys. Rev. C 84, 041601 (2011).
  * Lépine-Szily  _et al._ [2003] A. Lépine-Szily, J. M. Oliveira, D. Galante, G. Amadio, V. Vanin, R. Lichtenthäler, V. Guimarães, G. F. Lima, H. G. Bohlen, A. N. Ostrowski, A. Di Pietro, A. M. Laird, L. Maunoury, F. de Oliveira Santos, P. Roussel-Chomaz, H. Savajols, W. Trinder, A. C. C. Villari, and A. de Vismes, Nucl. Phys. A 722, C512 (2003).