# On the optimal relaxation rate for the Metropolis algorithm in one dimension

A. Patrón<EMAIL_ADDRESS>Física Teórica, Universidad de Sevilla, Apartado de
Correos 1065, E-41080 Sevilla, Spain  A. D. Chepelianskii
<EMAIL_ADDRESS>Laboratoire de Physique des
Solides, Université Paris-Saclay, CNRS, 91405, Orsay, France  A. Prados
<EMAIL_ADDRESS>Física Teórica, Universidad de Sevilla, Apartado de Correos 1065,
E-41080 Sevilla, Spain  E. Trizac  LPTMS, Université Paris-Saclay, CNRS,
91405, Orsay, France  Ecole Normale Supérieure de Lyon, F-69364 Lyon, France

###### Abstract

We study the relaxation of the Metropolis Monte Carlo algorithm corresponding
to a single particle trapped in a one-dimensional confining potential, with
even jump distributions that ensure that the dynamics verifies detailed
balance. Previous work suggested that, for smooth jump distributions, the
fastest relaxation rate is obtained as a result of the competition between
diffusive and rejection-dominated dynamics. In this work, we show that a new
regime comes into play for two-peaked jump distributions, where the relaxation
dynamics is neither dominated by diffusion nor rejection: the eigenmodes adopt
an oscillatory form, reminiscent of charge density waves (CDW)—thus we term
this new regime the CDW regime. Using a combination of numerical and
analytical techniques, the parameter regions corresponding to diffusion,
rejection, and CDW are characterised, as well as the transition lines between
them—i.e. a phase diagram is built. The optimal relaxation rate is located at
the triple point of phase coexistence, where the transition lines (diffusive-
rejection, diffusive-CDW, and CDW-rejection) intersect. Our theoretical
framework is checked versus the numerical diagonalisation of the master
equation. We also briefly discuss more sophisticated attempts at optimising
the relaxation rate to equilibrium.

##  I Introduction

Computational algorithms based on Monte Carlo techniques are widely known for
their versatility in solving many scientific problems. For instance, in
molecular simulations, they are essential for determining the equilibrium
properties of physical systems involving many degrees of freedom with high
accuracy [1]. Other branches of knowledge that range from the natural sciences
of physics, chemistry, and biology [2, 3] to the fields of economy, machine
learning, and archaeology [4, 5, 6] benefit from these techniques when
approaching problems at the deterministic and stochastic levels of
description, hence the interest in the study of the convergence of such
algorithms, as it provides useful insights on the optimal performance they may
reach.

The focus of our work concerns the Metropolis-Hastings algorithm, or simply
Metropolis [7], which belongs to the class of Markov Chain Monte Carlo
techniques. These are based on creating a sequence of correlated steps from a
random walk that, given sufficient time, makes it possible to sample a desired
target probability distribution—which usually corresponds to the equilibrium
distribution. In the Metropolis algorithm, the random displacements $\eta$ in
parameter space giving the new steps of the sequence are drawn from a so-
called jump distribution $\omega(\eta)$, which is an even function of the
displacements. Once a random displacement is chosen, the associated energy
change $\Delta U$ is calculated and the step is accepted or rejected according
to the Metropolis rule: if $\Delta U<0$, it is accepted with probability
$p=1$; if $\Delta U>0$, it is accepted with probability $p=\exp(-\beta\Delta
U)$, where $\beta$ is the inverse temperature. This rule verifies detailed
balance, and therefore it ensures that the system equilibrates in the long
time limit [8, 9]—provided that the jump distribution allows for the sampling
of all configurations, i.e. ergodicity is not broken.

A key issue concerns the overall convergence time; it quantifies the time
needed by the algorithm to attain the target distribution, which greatly
impacts its performance. The longer the convergence time, the larger the error
in the computed observables. On the one hand, if the typical size $a$ of the
attempted displacements is too small, most of the jumps are accepted but the
phase space is not explored sufficiently, thus involving a long convergence
time. On the other hand, if the size $a$ of these displacements is too large,
most jumps are rejected because they are likely to involve a high energy cost.
This discussion suggests that there must exist an optimal amplitude
$a_{\text{opt}}$ that minimises the convergence time for the Metropolis
algorithm. Although there have been many attempts to derive $a_{\text{opt}}$
for specific models [10, 11, 12], the general rule of thumb corresponds to
choosing $a$ so that the acceptance and rejection rates are comparable, i.e.
the probability of accepting or rejecting the new step is close to $0.5$ [1,
13, 14].

In previous work [15], the relaxation dynamics of the Metropolis algorithm was
studied for a particle trapped in a one-dimensional confining potential
$U(x)$, and the jump distribution that optimises the convergence rate to
equilibrium was investigated. For different potentials and choices of the jump
distribution, it was shown that a critical jump length $a^{*}$ appears where a
localisation transition occurs, which radically changes the relaxation
behaviour from $a<a^{*}$ to $a>a^{*}$. This localisation transition mainly
affects the eigenmodes of the master equation characterising the algorithm:
for $a<a^{*}$, the eigenmodes resemble the diffusive modes obtained at the
limit $a\rightarrow 0^{+}$, while for $a>a^{*}$, the eigenmodes become
spatially localised at discrete points $x_{0}$, where the rejection
probability is maximal. Interestingly, it was shown that
$a_{\text{opt}}=a^{*}$ for most of the cases considered, implying that the
optimal convergence time is obtained as a balance between the diffusive and
rejection modes.

Here, we show that the diffusion-rejection picture described above breaks down
when two-peaked jump distributions are considered for one-dimensional
potentials. In such a scenario, keeping the diffusion-rejection picture leads
to a collapse instability, in which the optimal jump distribution tends to a
double delta peak at $\eta=\pm a$. This is problematic, since this jump
distribution breaks ergodicity and thus does not allow the system to
equilibrate. We focus our study on the aforementioned phenomenon and reveal
that a new regime emerges, in which the spatial dependence of the eigenmodes
of the Metropolis algorithm varies drastically to an oscillatory behaviour,
reminiscent of the charge density wave (CDW) transition in correlated electron
systems [16, 17, 18]. Therefore, we will refer to this regime, which competes
with both the diffusion and rejection ones, as the CDW regime.

One of the main goals of this paper is to characterise the novel CDW regime
and analyse its implications for the optimisation of the convergence rate of
the Metropolis algorithm. To do so, we introduce two order parameters, the
fidelity $\mathcal{F}$ and the inverse participation ratio IPR, which allow us
to characterise the parameter regions dominated by diffusion, rejection, and
CDW. A key finding of our work is that the optimal convergence rate is placed
at the triple point of coexistence of the diffusive, rejection, and CDW
phases.

The structure of this paper is as follows. Section II presents the
mathematical framework behind the Metropolis algorithm. It also provides a
brief account of the diffusion-rejection picture from Ref. [15] and discusses
the emergence of the collapse instability. Section III is devoted to the
derivation and characterisation of the CDW regime. In addition, we numerically
check the emergence of these oscillatory modes. The phase diagram of the
Metropolis algorithm is analysed in Sec. IV, specifically for a bi-parametric
family of jump distributions. Then, in Sec. V we numerically show that
optimality is attained as the competition between the CDW, diffusion, and
rejection regimes. Section VI presents some preliminary ideas for further
optimising the convergence rate of the Metropolis algorithm. Finally, we
provide some concluding remarks in Sec. VII. The appendix discusses some
technicalities that are omitted in the main text.

##  II Relaxation dynamics for the Metropolis algorithm

In this work, we study the relaxation dynamics of a particle confined in a
one-dimensional confining potential $U(x)$, with $x$ accounting for its
position on the real line. At each discrete time step $n$, the particle’s
position evolves according to the so-called Metropolis rule

$x_{n}=\begin{dcases}x_{n-1}+\eta_{n}&\text{with prob.}\
p\equiv\text{min}\left(1,e^{-\beta\Delta U}\right),\\\ x_{n-1}&\text{with
prob.}\ 1-p,\end{dcases}$ (1)

where $\Delta U\equiv U(x_{n})-U(x_{n-1})$ and $\beta\equiv(k_{B}T)^{-1}$ is
the inverse temperature. Independent random variables $\eta_{n}$ are drawn
from a normalised distribution $\omega(\eta)$, which is usually termed the
jump distribution, note that $\Delta x=x_{n}-x_{n-1}=\eta_{n}$. We further
assume the jump distribution to be symmetric, i.e.
$\omega(-\eta)=\omega(\eta)$, which ensures that detailed balance holds.

At the ensemble level, the system is described via the position probability
distribution $P_{n}(x)$. Its dynamical evolution is governed by the master
equation

$P_{n}[\omega](x)=\int_{-\infty}^{+\infty}dx^{\prime}\
F[\omega](x,x^{\prime})P_{n-1}[\omega](x^{\prime}),$ (2)

in which

$F[\omega](x,x^{\prime})\equiv\delta(x-x^{\prime})R[\omega](x)+\omega(x-x^{\prime})\left[\theta(U(x^{\prime})-U(x))+e^{-\beta(U(x)-U(x^{\prime}))}\theta(U(x)-U(x^{\prime}))\right],$
(3)

is the (temperature-dependent 111To simplify our notation, we do not show the
explicit dependence on the temperature. In simple cases like the harmonic
potential, the temperature can be absorbed in the definition of the unit of
length.) integral kernel of the master equation, and

$R[\omega](x)\equiv\int_{-\infty}^{+\infty}dx^{\prime}\
\omega(x^{\prime}-x)\left[1-e^{-\beta(U(x^{\prime})-U(x))}\right]\theta(U(x^{\prime})-U(x)).$
(4)

We remark that $R[\omega](x)$ gives the rejection probability at position $x$,
which is a central quantity in this work. In the above expressions,
$\theta(x)$ corresponds to the Heaviside step-function; $\theta(x)=1$ for
$x\geq 0$ and $\theta(x)=0$ for $x<0$. Our notation emphasises the
(functional) dependence of $P_{n}$, $F$, and $R$ on the jump distribution
$\omega(\eta)$, in addition to their spatial dependence.

The master equation (2) has the canonical equilibrium distribution

$P_{\text{eq}}(x)=\frac{1}{Z}e^{-\beta U(x)},\quad
Z\equiv\int_{-\infty}^{+\infty}dx\ e^{-\beta U(x)},$ (5)

as a stationary solution. The Metropolis rule (1) ensures that the equilibrium
distribution is always attained in the long time limit $n\to\infty$,
regardless of the shape of the initial probability distribution $P_{0}(x)$,
i.e.

$P_{\infty}(x)\equiv\lim_{n\to\infty}P_{n}[\omega](x)=P_{\text{eq}}(x),$ (6)

provided that the jump distribution allows for the sampling of all positions
$x$, i.e. the system is ergodic.

We stress that our analysis on the convergence rate of the Metropolis
algorithm is general, encompassing a wide class of jump distributions
$\omega(\eta)$. Still, for illustration purposes, we consider some specific
choices—the following bi-parametric families of jump distributions—throughout
this work:

  * •

The $(a,\sigma)$-family of Gaussian jump distributions

$\omega^{(a,\sigma)}(\eta)=C^{(a,\sigma)}\text{exp}\left(-\frac{(|\eta|-a)^{2}}{2\sigma^{2}}\right),$
(7)

where $C^{(a,\sigma)}$ is a normalisation constant, and $a,\sigma>0$.

  * •

The $(a,\alpha)$-family of algebraic jump distributions

$\omega^{(a,\alpha)}(\eta)=\frac{\alpha+1}{2a^{\alpha+1}}|\eta|^{\alpha}\theta(a-|\eta|),$
(8)

with $a,\alpha>0$. In the above, $\alpha$ accounts for the width of the
distribution, while $a$ stands for the position of its maximum.

  * •

The $(a,b)$-family of box-like jump distributions:

$\omega^{(a,b)}(\eta)=\frac{1}{4b}\theta(|\eta|-a+b)\theta(a+b-|\eta|),\quad
a>b>0.$ (9)

This distribution has a simple physical interpretation: only jumps with
$a-b<|\eta|<a+b$ are allowed, with a flat distribution thereof.

All these families tend to a double Dirac-delta peak

$\omega_{0}(\eta)=\frac{1}{2}\Big{[}\delta(\eta-a)+\delta(\eta+a)\Big{]},$
(10)

when their respective widths tend to zero, i.e. when (i) $\sigma\to 0^{+}$ in
Eq. (7), (ii) $\alpha\to+\infty$ in Eq. (8), and (iii) $b\to 0^{+}$ in Eq.
(9). The double delta peak does not allow for the sampling of all positions,
since the length of the jumps is fixed to $\pm a$. This will be important in
connection with the collapse instability discussed at the end of Sec. II.2.

###  II.1 Convergence rate

To study the relaxation rate to equilibrium, we put our focus on the spectrum
of eigenvalues of the integral kernel $F[\omega](x,x^{\prime})$ of the master
equation (2). For that purpose, it is useful to introduce a new integral
kernel $K[\omega](x,x^{\prime})$ that share the spectrum of eigenvalues with
$F[\omega](x,x^{\prime})$,

$K[\omega](x,x^{\prime})\equiv
e^{-\beta(U(x)-U(x^{\prime}))/2}F[\omega](x,x^{\prime})=\omega(x-x^{\prime})e^{-\beta|U(x)-U(x^{\prime})|/2}+\delta(x-x^{\prime})R[\omega](x).$
(11)

Moreover, $K[\omega](x,x^{\prime})$ is symmetric under the exchange of
$x\leftrightarrow x^{\prime}$—note the absolute value in the exponential
factor multiplying $\omega(x-x^{\prime})$, which makes it self-adjoint and
simplifies our mathematical description below.

Now we define the eigenvalues $\lambda_{\nu}$ and eigenfunctions
$\phi_{\nu}(x)$ of the integral operator $\hat{K}$,

$\hat{K}[\omega]\phi_{\nu}(x)\equiv\int_{-\infty}^{+\infty}dx^{\prime}\
K[\omega](x,x^{\prime})\,\phi_{\nu}(x^{\prime})=\lambda_{\nu}[\omega]\,\phi_{\nu}[\omega](x),\quad\nu=0,1,2,\ldots$
(12)

Since $K[\omega](x,x^{\prime})$ is self-adjoint, all the eigenvalues of Eq.
(12) are real. Moreover, for jump distributions $\omega$ such that the Perron-
Frobenius theorem [20] applies, this theorem ensures that (i) the maximum
eigenvalue is $\omega$ independent, $\lambda_{0}=1$, which corresponds to the
stationary solution being nondegenerate, and (ii) the remaining eigenvalues
$\lambda_{\nu}[\omega]$, with $\nu>0$, verify the strict inequality
$|\lambda_{\nu}|<1$. The eigenfunction $\phi_{0}(x)$ for the eigenvalue
$\lambda_{0}=1$, corresponding to the equilibrium distribution, is also
$\omega$ independent,

$\phi_{0}(x)\propto e^{-\beta U(x)/2}.$ (13)

Explicitly, we can write

$\int_{-\infty}^{+\infty}dx^{\prime}\
K[\omega](x,x^{\prime})\phi_{0}(x^{\prime})=\phi_{0}(x),$ (14a)
$\int_{-\infty}^{+\infty}dx^{\prime}\
K[\omega](x,x^{\prime})\phi_{\nu}[\omega](x^{\prime})=\lambda_{\nu}[\omega]\,\phi_{\nu}[\omega](x),\quad\nu>0.$
(14b)

Our notation emphasises the fact that, under the conditions of the Perron-
Frobenius theorem, both $\lambda_{0}$ and $\phi_{0}$ are independent of the
jump distribution $\omega$, whereas the remainder of the eigenvalues and
eigenfunctions are functionals thereof. Making use of the explicit form of the
kernel $K[\omega](x,x^{\prime})$ in Eq. (11), we have that

$\int_{-\infty}^{+\infty}dx^{\prime}\
\omega(x-x^{\prime})e^{-\beta|U(x)-U(x^{\prime})|/2}\phi_{\nu}(x^{\prime})+R[\omega](x)\,\phi_{\nu}[\omega](x)=\lambda_{\nu}[\omega]\,\phi_{\nu}[\omega](x),\quad\nu>0.$
(15)

The solution of the master equation can always be written as a linear
combination of the eigenfunctions $\phi_{\nu}(x)$:

$P_{n}[\omega](x)=P_{\text{eq}}(x)+e^{-\beta
U(x)/2}\sum_{\nu>0}\mathcal{A}_{\nu}[\omega]\,\phi_{\nu}[\omega](x)\,(\lambda_{\nu}[\omega])^{n},$
(16)

where $\mathcal{A}_{\nu}[\omega]$ are constants determined by the initial
condition $P_{0}(x)$. For long times, the system always reaches equilibrium
because $|\lambda_{\nu}[\omega]|<1$ for $\nu>0$. The leading convergence rate
to equilibrium of the Metropolis algorithm is thus given by the next-to-
largest eigenvalue $\lambda_{1}[\omega]$, since

$P_{n}[\omega](x)-P_{\text{eq}}(x)\sim e^{-\beta
U(x)/2}\mathcal{A}_{1}[\omega]\,\phi_{1}[\omega](x)\,(\lambda_{1}[\omega])^{n},\quad
n\gg 1.$ (17)

Analytically determining both $\lambda_{1}[\omega]$ and its corresponding
eigenfunction $\phi_{1}[\omega](x)$ from the integral equation (12) is
generally not possible. Still, the developed framework makes it possible to
estimate the convergence rate with a variational procedure, as explained
below.

The next-to-largest eigenvalue can always be written as

$\Lambda[\omega]\equiv\lambda_{1}[\omega]=\underset{\phi\perp\phi_{0}}{\max}\
\int dx\int dx^{\prime}\ \phi(x)K[\omega](x,x^{\prime})\phi(x^{\prime}).$ (18)

Our final goal is to select the optimal $\omega^{*}(\eta)$ that minimises
$\Lambda[\omega]$, in order to find the optimal convergence rate to
equilibrium

$\Lambda^{*}=\min_{\omega}\Lambda[\omega].$ (19)

Therefore, we are actually dealing with a min-max problem. The minimisation in
Eq. (19) should ideally be carried out over the whole space of possible jump
distributions but, for practical purposes, it may also be performed within a
given family of jump distributions, such as those in Eqs. (7)-(9).

Now, we may resort to a variational approach: Let us consider a certain known
basis of orthogonal functions $\psi_{n}(x)$, and expand the function $\phi(x)$
in Eq. (18) in this basis. Truncating the expansion after $N_{s}$ states,
which must be perpendicular to the ground state, we get a lower bound for
$\Lambda[\omega]$, i.e.

$\Lambda[\omega]\geq\tilde{\Lambda}[\omega]\equiv\underset{c_{1},c_{2},...,c_{N_{s}}}{\max}\
\int dx\int dx^{\prime}\
\phi(x)K[\omega](x,x^{\prime})\phi(x^{\prime}),\quad\phi(x)=\sum_{n=1}^{N_{s}}c_{n}\psi_{n}(x).$
(20)

In order to calculate our estimate $\tilde{\Lambda}[\omega]$ for the
convergence rate, we must find the optimal values of the coefficients
$\\{c_{n}\\}_{n=1}^{N_{s}}$ that maximise the integral in Eq. (20). Such
problem is equivalent to finding the largest eigenvalue of the reduced
$N_{s}\times N_{s}$ matrix $\tilde{K}[\omega]$ with elements

$\tilde{K}_{nm}[\omega]=\int dx\int dx^{\prime}\
\psi_{n}(x)K[\omega](x,x^{\prime})\psi_{m}(x^{\prime}),$ (21)

corresponding to the truncation $1\leq n,m\leq N_{s}$. Given a fixed value of
$N_{s}$, the choice of the basis affects how the lower bound of Eq. (20)
approaches the true eigenvalue from Eq. (18). It is worth emphasising that for
symmetric confining potentials, the eigenbasis of our problem may be split
into symmetric and antisymmetric eigenfunctions. The integral kernel
$K[\omega](x,x^{\prime})$ inherits these symmetries, implying that the matrix
elements $\tilde{K}_{nm}[\omega]$ are zero for those values of $n$ and $m$ for
which the corresponding eigenfunctions have different parity. Thus, we may
divide the basis into two bases: one containing the even eigenfunctions, and
the other one containing the odd eigenfunctions. In this way, the matrix
$\tilde{K}_{nm}$ can be written as a direct sum of even and odd sectors.

Using different bases, we can ensure that our lower bound is sufficiently
close to the true convergence rate for every possible jump distribution
$\omega(\eta)$. In fact, using $M$ different bases, we may define the
functional

$\Lambda_{M}[\omega]\equiv\max\left[\Lambda^{(1)}[\omega],\Lambda^{(2)}[\omega],...,\Lambda^{(M)}[\omega]\right]$
(22)

to estimate the convergence rate, for a sufficiently high number of retained
eigenfunctions $N_{s}$ for each basis. Throughout this work, we mainly use the
estimator above with different values of $M$, typically with $N_{s}=100$
retained eigenfunctions per basis. In order to approximately solve the min-max
problem presented above, we thus resort to a two-step procedure. First, for
each considered jump distribution $\omega(\eta)$, we use the functional (22)
to find an accurate estimate for the leading eigenvalue $\lambda_{1}[\omega]$,
i.e. the maximum $\Lambda_{M}[\omega]$ of its lower bounds
$\Lambda^{(i)}[\omega]$, $i=1,\ldots,M$, for the $M$ bases considered. Second,
we vary the jump distribution to find the minimum value of such estimates
$\Lambda_{M}^{*}$, which is attained for a specific jump distribution
$\omega_{M}^{*}$, i.e. $\Lambda_{M}^{*}\equiv\Lambda_{M}[\omega_{M}^{*}]$.

One could use bases without a clear physical meaning and the lower bound
obtained in this way might still be quite useful for numerical purposes, but
this “blind” procedure would not provide meaningful insights into the
behaviour of the eigenfunctions of the master equation and the physics
underlying the optimal value for the convergence rate.

###  II.2 Localisation transition and collapse instability

In order to get more information about the physics of the problem, a
sufficiently narrow jump distribution $\omega(\eta)=\delta_{\sigma}(\eta)$
around $\eta=0$ with standard deviation $\sigma$ was considered in Ref. [15].
In this way, Eq. (12) could be expanded around $x^{\prime}=x$ due to the
smallness of $\sigma$, giving as a result a Schrödinger-like equation with an
effective potential,

$\sigma^{2}\left\\{-\frac{\phi_{n}^{\prime\prime}(x)}{2}+\left[\frac{\beta^{2}}{8}(U^{\prime}(x))^{2}-\frac{\beta}{4}U^{\prime\prime}(x)\right]\phi_{n}(x)\right\\}=(1-\lambda_{n})\phi_{n}(x).$
(23)

The eigenmodes $\phi_{n}(x)$ are the eigenfunctions of this Schrödinger-like
equation, which have a well-defined parity for even potentials, such that
$U(x)=U(-x)$. This makes it possible to split this Schrödinger basis into two
parts that contain, respectively, the even and odd eigenfunctions, as already
stated in Sec. II.1.

In Ref. [15], it was also argued that the convergence rate is limited by the
rejection probability $R[\omega](x)$. The main idea is the following: starting
from an initial distribution $\delta(x-x_{0})$ at step $n=0$, the probability
distribution function at any later time is

$P_{n}(x)=R[\omega](x_{0})^{n}\delta(x-x_{0})+p_{n}(x|x_{0}),$ (24)

where $p_{n}(x|x_{0})$ is a smooth function. The system cannot relax faster
than this exponential decay, i.e., the convergence rate should always be
larger than $R[\omega](x_{0})$, for any $x_{0}$, i.e.
$\Lambda\geq\underset{x}{\max}\ R[\omega](x)$. It is only upon the
localisation transition that the above inequality saturates [15].

Therefore, the convergence rate $\Lambda$ was estimated in Ref. [15] by using
the 3-bases functional

$\Lambda_{3}[\omega]\equiv\max\left[\Lambda^{(\text{Sch.})}_{\text{o}}[\omega],\Lambda^{(\text{Sch.})}_{\text{e}}[\omega],\underset{x}{\max}\
R[\omega](x)\right],$ (25)

where $\Lambda^{(\text{Sch.})}_{\text{o}}[\omega]$ and
$\Lambda^{(\text{Sch.})}_{\text{e}}[\omega]$ are the estimates for $\Lambda$
obtained from the odd and even sectors of the Schrödinger basis, respectively.
Estimate (25) is based on the assumption that the convergence rate is obtained
as a balance between the diffusion modes of the Schrödinger equation and the
maximum value of the rejection rate. The validity of such an argument was
checked for several choices of the confining potential $U(x)$ and different
jump distributions $\omega(\eta)$, even for a few number $N_{s}$ of retained
eigenmodes.

The main result of Ref. [15] is the optimal convergence rate being obtained at
the transition point between the diffusion and rejection regimes, which is
where the localisation transition takes place. The term localisation stems
from the fact that, in addition to the discrete spectrum of well-defined
eigenfunctions, a continuum of eigenvalues appears at the rejection regime,
with singular localised eigenfunctions—as suggested by Eq. (24), these
singular eigenfunctions are peaked at any point $x_{0}$, with eigenvalue
$R[\omega](x_{0})$.

Here we show that, although the above diffusion-rejection picture is quite
robust for a wide class of jump distributions, this relatively simple scenario
fails when considering two-peaked jump distributions—regardless of the chosen
confining potential. Let us consider the above optimisation problem for a
particle confined in a harmonic potential,

$U(x)=\chi x^{2}/2,$ (26)

with stiffness $\chi$. Specifically, we employ the bi-parametric $(a,\sigma)$
family of Gaussian jump distributions defined in Eq. (7) for the Metropolis
algorithm. That is, we look for the optimal values of the parameters
$(a^{*},\sigma^{*})$ that minimise the leading eigenvalue $\Lambda$.

In the numerical scheme for the optimisation 222The following numerical
procedure has been employed for all the examples presented throughout the
paper. To obtain the eigenvalues and eigenfuntions of the master equation, we
have employed the FORTRAN 90 subroutine dsyev.f90 from LAPACK [53]. To avoid
spurious boundary effects, the spatial variable $x$ is restricted to a finite
interval $[-L,L]$, with a sufficiently large $L$. Moreover, space is
discretised with a mesh of $N_{d}$ points. For the harmonic potential, we have
fixed units by choosing $k_{B}T=1$ and $\chi=1$, and taken $L=10$ and
$N_{d}=2000$. To numerically estimate $\Lambda$, as defined in Eq. (25), we
have employed the open source nonlinear optimisation libraries from NLopt [54]
with $N_{s}=100$ modes for each basis., a lower cut-off in the width $\sigma$
must be imposed, i.e. $\sigma\geq\sigma_{\text{cutoff}}$, since very small
values of $\sigma$ lead to numerical difficulties 333Since the numerics
involve a spatial discretisation, $\sigma$ must be larger than the lattice
spacing in this discretisation.. If, when decreasing $\sigma_{\text{cutoff}}$,
$\sigma^{*}$ consistently reached the cut-off value, this would hint at an
incomplete convergence of $\sigma$ due to the cut-off, i.e. it would hint at
actually having $\sigma^{*}\to 0$. We refer to this problematic scenario as
the collapse instability, since the optimal jump distribution seems to
converge to a double Dirac-delta peak, which breaks ergodicity as discussed at
the end of Sec. II.

The failure of the rejection-diffusion picture is illustrated in Fig. 1. We
employ three different values for the cutoff, namely
$\sigma_{\text{cutoff}}=0.5$, $0.2$, and $0.05$. Using the diffusion-rejection
estimator (25), we find the corresponding optimal jump distribution
$\omega_{3}^{*}(\eta)=\omega^{(a_{3}^{*},\sigma_{3}^{*})}(\eta)$, and the
associated convergence rate
$\Lambda_{3}^{*}\equiv\Lambda_{3}[\omega_{3}^{*}]$. The values obtained for
each cutoff are shown in the legend. As $\sigma_{\text{cutoff}}$ is decreased,
$a^{*}$ goes to a definite value, around $1.38$, but $\sigma_{3}^{*}$ always
reaches the cutoff, $\sigma^{*}=\sigma_{\text{cutoff}}$. The optimal jump
distribution thus shows the collapse instability, tending to a Dirac-delta in
the limit $\sigma\to 0$, and the diffusion-rejection picture breaks down. This
hints at the emergence of new physics: clarifying this issue is the main goal
of the next section. We defer the discussion of the relevance of the values
found for the convergence rate $\Lambda_{3}^{*}$ to Sec. V.

Figure 1: Optimal jump distributions for the bi-parametric Gaussian family in
Eq. (7). The plotted curves correspond to three different values for the
cutoff, $\sigma_{\text{cutoff}}=0.5$ (dotted), $0.2$ (dashed), and $0.05$
(solid). In all cases, the numerically obtained optimal jump distribution
obtained from the functional $\Lambda_{3}[\omega]$ verifies
$\sigma_{3}^{*}=\sigma_{\text{cutoff}}$—thus showing the collapse instability
phenomenon.

##  III Collapse instability and charge density waves

It is the tendency of the optimal jump distribution towards a Dirac-delta-like
distribution in our optimisation scheme that we term the collapse instability.
As just discussed at the end of Sec. II, this is problematic: Dirac-delta jump
distributions break ergodicity in the one-dimensional case and the system does
not equilibrate. In this section, we provide a rigorous characterisation of
the intermediate CDW regime that emerges for two-peaked jump distributions in
1D, which solves this issue.

###  III.1 Analytical derivation of the CDW eigenfunctions from the master
equation

First, we derive the corresponding CDW eigenfunctions directly from the master
equation and study their properties. We consider here a two-peak family of
jump distributions $\omega(\eta;a,\sigma)$, the maxima of which are at
$\eta=\pm a$ with a certain small width, measured by a parameter $\sigma$. The
parameters $(a,\sigma)$ define the two characteristic lengths of the jump
distribution. Let us introduce the auxiliary distribution

$\delta_{\sigma}(\eta-a)\equiv 2\omega(\eta)\theta(\eta).$ (27)

By definition, (i) the maximum of $\delta_{\sigma}(\eta)$ is located at
$\eta=0$ and (ii) $\delta(\eta)$ identically vanishes for $\eta<-a$. The width
of the peak of $\delta_{\sigma}$ is measured by its standard deviation. Using
this definition, we write the jump distribution as follows:

$\omega(\eta;a,\sigma)=\frac{1}{2}\left[\delta_{\sigma}(\eta-a)+\delta_{\sigma}(-\eta-a)\right].$
(28)

For $\sigma=0$, $\delta_{0}(\eta)$ corresponds to a Dirac-delta distribution
and Eq. (28) reduces to Eq. (10), which—as discussed above—breaks ergodicity
and does not allow the system to equilibrate. Here, we will be mainly
interested in the regime $\sigma\ll a$—but non-zero, where the jump
distribution is narrowly peaked around $\pm a$, and thus
$\omega(0;a,\sigma)=\delta_{\sigma}(-a)\approx 0$.

Now we proceed to analyse the eigenfunctions $\phi_{\nu}(x;a,\sigma)$ of the
master equation (2), i.e. the non-trivial solutions of Eq. (15). It is
convenient to introduce the change of variables

$\phi_{\nu}(x;a,\sigma)=\phi_{0}(x)\rho_{\nu}(x;a,\sigma)\propto e^{-\beta
U(x)/2}\rho_{\nu}(x;a,\sigma),$ (29)

such that the equilibrium distribution corresponds to a constant,
$\rho_{0}(x;\cancel{a,\sigma})=\text{const.}$ Note that functionals of
$\omega$ translate here into a parametric dependence on $(a,\sigma)$, which we
omit henceforth to simplify the notation. For the sake of concreteness and
simplicity, we present the analysis for confining potentials $U(x)$ that are
convex, i.e. $U^{\prime\prime}(x)>0$, and thus have only one minimum—without
loss of generality, we assume that the minimum is at $x=0$. We stress that the
convexity condition can be relaxed and our conclusions still hold, but it
allows us to simplify the presentation while keeping the main ideas.

We start by analysing the second term on the rhs of Eq. (15), which involves
the rejection probability defined in Eq. (4). By defining $\ell(x)$ and $h(x)$
as the $x^{\prime}$-values at which $U(x^{\prime})=U(x)$, the Heaviside
function restricts the integral over $x^{\prime}$ to a certain interval
$(-\infty,\ell(x))\cup(h(x),+\infty)$. Depending on the sign of
$U^{\prime}(x)$, either $\ell(x)$ or $h(x)$ equals $x$—in the following we
focus on the latter case, $U^{\prime}(x)>0$ ($x>0$), where $h(x)=x$. With
these definitions,

$R(x;a,\sigma)=\left(\int_{-\infty}^{\ell(x)}+\int_{x}^{+\infty}\right)dx^{\prime}\omega(x^{\prime}-x;a,\sigma)\left[1-e^{-\beta\left(U(x^{\prime})-U(x)\right)}\right].$
(30)

The first integral is subdominant, since $\omega$ is basically zero over the
considered interval, as long as $\ell(x)-x<-a$, i.e. $x-a>\ell(x)$. On the
other hand, in the second integral we have that $x^{\prime}-x>0$ and therefore
only one of the $\delta_{\sigma}$’s contribute:

$R(x;a,\sigma)\approx\frac{1}{2}\int_{x}^{+\infty}dx^{\prime}\delta_{\sigma}(x^{\prime}-x-a)\left[1-e^{-\beta\left(U(x^{\prime})-U(x)\right)}\right]\approx\frac{1}{2}\left[1-e^{-\beta\left(U(x+a)-U(x)\right)}\right]\approx\frac{1}{2},$
(31)

assuming that

$e^{-\beta\left(U(x+a)-U(x)\right)}\ll 1.$ (32)

This condition basically implies that $a$ is large enough, see Sec. III.2 for
more details. Now we move on to analysing the first term on the rhs of Eq.
(15), splitting the integral into two slices: one in which
$U(x)>U(x^{\prime})$, i.e., the interval $(\ell(x),x)$ and another in which
$U(x)<U(x^{\prime})$, i.e. $(-\infty,\ell(x))\cup(x,+\infty)$. For the first
slice, $U(x)>U(x^{\prime})$, we have

$\displaystyle\int_{-\infty}^{+\infty}dx^{\prime}$
$\displaystyle\theta(U(x)-U(x^{\prime}))\,\omega(x-x^{\prime};a,\sigma)\,e^{-\beta|U(x)-{U(x^{\prime})}|/2}\rho_{\nu}(x^{\prime})e^{-\beta
U(x^{\prime})/2}$
$\displaystyle=\int_{\ell(x)}^{x}dx^{\prime}\omega(x-x^{\prime};a,\sigma)e^{-\beta(U(x)-\cancel{U(x^{\prime})})/2}\rho_{\nu}(x^{\prime})e^{-\beta\cancel{U(x^{\prime})}/2}\approx\frac{1}{2}e^{-\beta
U(x)/2}\int_{\ell(x)}^{x}dx^{\prime}\delta_{\sigma}(x-x^{\prime}-a)\rho_{\nu}(x^{\prime})$
$\displaystyle\approx\frac{1}{2}e^{-\beta
U(x)/2}\int_{-\infty}^{+\infty}dx^{\prime}\delta_{\sigma}(x-x^{\prime}-a)\rho_{\nu}(x^{\prime})=\frac{e^{-\beta
U(x)/2}}{2}\int_{-\infty}^{+\infty}dx^{\prime}\delta_{\sigma}(x-x^{\prime})\rho_{\nu}(x^{\prime}-a)$
$\displaystyle\approx\frac{e^{-\beta
U(x)/2}}{2}\int_{-\infty}^{+\infty}dx^{\prime}\delta_{\sigma}(x-x^{\prime})\left[\rho_{\nu}(x-a)+\frac{1}{2}\rho_{\nu}^{\prime\prime}(x-a)(x-x^{\prime})^{2}\right]$
$\displaystyle=\frac{1}{2}e^{-\beta
U(x)/2}\left[\rho_{\nu}(x-a)+\frac{\sigma^{2}}{2}\rho_{\nu}^{\prime\prime}(x-a)\right],$
(33)

where we have taken into account again that (i) only one of the
$\delta_{\sigma}$’s contribute, and (ii) the interval of integration may be
extended from $(\ell(x),x)$ to $(-\infty,+\infty)$, once the peak of the
$\delta_{\sigma}$ is included in the considered interval—again, as long as
$x-a>\ell(x)$. For the second slice, $U(x^{\prime})>U(x)$, we have

$\displaystyle\int_{-\infty}^{+\infty}dx^{\prime}$
$\displaystyle\theta(U(x^{\prime})-U(x))\,\omega(x-x^{\prime};a,\sigma)\,e^{-\beta|U(x)-{U(x^{\prime})}|/2}\rho_{\nu}(x^{\prime})e^{-\beta
U(x^{\prime})/2}$
$\displaystyle=\left(\int_{-\infty}^{\ell(x)}+\int_{x}^{+\infty}\right)dx^{\prime}\omega(x-x^{\prime};a,\sigma)e^{-\beta(U(x^{\prime})-U(x))/2}\rho_{\nu}(x^{\prime})e^{-\beta
U(x^{\prime})/2}$ $\displaystyle\approx\frac{e^{\beta
U(x)/2}}{2}\int_{x}^{+\infty}dx^{\prime}\delta_{\sigma}(x+a-x^{\prime})\rho_{\nu}(x^{\prime})e^{-\beta
U(x^{\prime})}\approx\frac{e^{\beta
U(x)/2}}{2}\int_{-\infty}^{+\infty}dx^{\prime}\delta_{\sigma}(x+a-x^{\prime})\rho_{\nu}(x^{\prime})e^{-\beta
U(x^{\prime})}$ $\displaystyle\approx
e^{-\beta(U(x+a)-U(x))}\frac{1}{2}e^{-\beta
U(x)/2}\rho_{\nu}(x+a)+O(\sigma^{2}).$ (34)

Note that Eq. (34) is subdominant against Eq. (33), due to the presence of the
exponential factor (32), so it is neglected in the following. The above
analysis implies that the eigenvalue equation (15) can be approximated by

$\frac{1}{2}\cancel{e^{-\beta
U(x)/2}}\left[\rho_{\nu}(x-a)+\frac{\sigma^{2}}{2}\rho_{\nu}^{\prime\prime}(x-a)\right]+\frac{1}{2}\cancel{e^{-\beta
U(x)/2}}\rho_{\nu}(x)=\lambda_{\nu}\cancel{e^{-\beta U(x)/2}}\rho_{\nu}(x),$
(35)

where $U^{\prime}(x)>0$, i.e. $x>0$, provided that $x-a>\ell(x)$. A similar
analysis for $U^{\prime}(x)<0$ leads to a completely analogous expression for
$x<0$, with the change $a\to-a$, provided that $x+a<h(x)$—note that
$\ell(x)=x$ for $x<0$.

Summarising our findings, in the limit $\sigma\ll a$ we have the approximate
equations

$\frac{\sigma^{2}}{4}\rho_{\nu}^{\prime\prime}(x-a)+\frac{1}{2}\Big{[}\rho_{\nu}(x-a)+\rho_{\nu}(x)\Big{]}=\lambda_{\nu}\rho_{\nu}(x),\quad
x>0,\;x-a>\ell(x);$ (36a)
$\frac{\sigma^{2}}{4}\rho_{\nu}^{\prime\prime}(x+a)+\frac{1}{2}\Big{[}\rho_{\nu}(x+a)+\rho_{\nu}(x)\Big{]}=\lambda_{\nu}\rho_{\nu}(x),\quad
x<0,\;x+a<h(x).$ (36b)

For even potentials, $\ell(x)=-x<0$ for $x>0$ ($h(x)=-x>0$ for $x<0$). Then,
for even potentials Eq. (36a) (Eq. (36b)) is limited to $x>a/2$ ($x<-a/2$). A
more careful analysis would be needed in the region $-a/2<x<a/2$ for even
potentials—and a certain region of small $x$ for potentials without symmetry.
This could have been anticipated after the approximations made in Eq. (33):
whereas the starting expression clearly vanishes in the limit as $x\to 0$, the
final approximate expression does not.

As a check of the consistency of Eqs. (36), note that a uniform solution
corresponds to the equilibrium eigenvalue $\lambda_{0}=1$. Moreover, for
infinitely sharp jump distributions, i.e. in the case $\sigma=0$, the
Metropolis algorithm loses ergodicity and any function $\phi_{\nu}(x)$ of the
form (29) with $\rho_{\nu}(x)$ being periodic with period $a$ is expected to
be a stationary solution of the master equation, i.e. an eigenfunction
corresponding to the eigenvalue $\lambda_{0}=1$. Note that this is indeed
predicted by Eqs. (36), as setting $\lambda_{\nu}=1$ for $\sigma=0$ yields
$\rho_{\nu}(x\pm a)=\rho_{\nu}(x)$. Interestingly, Eqs. (36) still admit
spatially periodic solutions, with period $a$, for $\sigma\neq 0$. Therefore,
our theoretical approach below can be viewed as a singular perturbation theory
in this subspace. By inserting the charge density wave (CDW) ansatz

$\rho_{\nu}(x\pm a)=\rho_{\nu}(x),\quad\forall x,$ (37)

we get the equation

$\rho_{\nu}^{\prime\prime}(x)+\frac{4(1-\lambda_{\nu})}{\sigma^{2}}\rho_{\nu}(x)=0,\quad\forall
x,$ (38)

which provides the solutions

$\rho_{\nu}(x)=c_{\nu}\cos\left(\frac{2\pi\nu}{a}x\right)+d_{\nu}\sin\left(\frac{2\pi\nu}{a}x\right),\quad\lambda_{\nu}=1-\nu^{2}\pi^{2}\left(\frac{\sigma}{a}\right)^{2}.$
(39)

Note that these CDW modes only make sense for $\sigma\neq 0$: as stated above,
for $\sigma=0$ any $a$-periodic $\rho_{\nu}(x)$ belongs in the subspace
corresponding to $\lambda_{0}=1$, which becomes degenerate due to the breaking
of ergodicity. In addition, the CDW modes cease to make sense for very large
values of $\nu$: we know that the exact eigenvalues of the master equation
cannot become smaller than $-1$, whereas
$\lim_{\nu\to\infty}\lambda_{\nu}\to-\infty$ 444A similar issue was found in
Ref. [15] for the Schrödinger eigenvalues..

For $\sigma\neq 0$, we have that $\lambda_{0}=1$ is clearly non-degenerate,
since the sine mode identically vanishes for $\nu=0$ and the cosine mode
becomes a constant. On the other hand, all the other eigenvalues
$\lambda_{\nu}$, with $\nu>0$, seem to be degenerate, since both the sine and
the cosine modes are non-zero. Henceforth, we write

$\phi_{\nu}^{\text{CDW},e}(x)=c_{\nu}^{e}\,\phi_{0}(x)\cos\left(k_{\nu}x\right)\propto
e^{-\beta U(x)/2}\cos\left(k_{\nu}x\right),$ (40a)
$\phi_{\nu}^{\text{CDW},o}(x)=c_{\nu}^{o}\,\phi_{0}(x)\sin\left(k_{\nu}x\right)\propto
e^{-\beta U(x)/2}\sin\left(k_{\nu}x\right),$ (40b)

with

$k_{\nu}=\frac{2\pi\nu}{a},\quad\lambda_{\nu}=1-\frac{\sigma^{2}}{4}k_{\nu}^{2},$
(41)

for the CDW modes that approximate the exact eigenfunctions $\phi_{\nu}(x)$ of
the master equation for $\sigma\ll a$. In the above, the indexes $e$ and $o$
refer to the even and odd sectors of the CDW basis, respectively. The
normalisation constant $c_{\nu}$ is determined by imposing
$\norm{\phi_{\nu}^{\text{CDW}}}=1$, i.e.

$\int_{-\infty}^{+\infty}dx\,\left|\phi_{\nu}^{\text{CDW}}(x)\right|^{2}=1.$
(42)

###  III.2 Robustness and numerical check of the CDW eigenfunctions

Now we proceed to numerically check the validity of our approximate derivation
of the CDW eigenfunctions. Before doing that, it is convenient to remark on
the limitations of our derivation, in order to have better insights when
comparing with the numerics. Here, we will restrict ourselves to the case of
even potentials, since the numerical check is done for the harmonic potential
(26).

First, it must be stressed that Eq. (36a) is only valid for $x>a/2$ and Eq.
(36b) for $x<-a/2$, and we have extended the equation to all $x$. Then, we
have that $a$ must be “small” enough for the CDW to emerge. On the other hand,
Eq. (32) entails that $a$ must be “large” enough. Since
$U^{\prime\prime}(x)>0$, for $x>0$ one has

$1\gg e^{-\beta(U(x+a)-U(x))}=e^{-\beta aU^{\prime}(\tilde{x})}>e^{-\beta
aU^{\prime}(x+a)},$ (43)

where $\tilde{x}\in[x,x+a]$. The convexity condition implies that it is enough
to ensure the above relation for $x=0$, i.e. it is enough to have that
$e^{-\beta aU^{\prime}(a)}\ll 1$. For example, in the case of the harmonic
potential, this is equivalent to having $e^{-\beta\chi a^{2}}\ll 1$.

Second, we may consider the higher-order terms in the width of the
distribution that have been neglected in Eq. (33). This would result in the
inclusion of the $O(\sigma^{4})$ contribution into the term in brackets in Eq.
(35), i.e.

$\frac{1}{2}\cancel{e^{-\beta
U(x)/2}}\left[\rho_{\nu}(x-a)+\frac{\sigma^{2}}{2}\rho_{\nu}^{\prime\prime}(x-a)+\xi\frac{\sigma^{4}}{4!}\rho_{\nu}^{\prime\prime\prime\prime}(x-a)\right]+\frac{1}{2}\cancel{e^{-\beta
U(x)/2}}\rho_{\nu}(x)=\lambda_{\nu}\cancel{e^{-\beta U(x)/2}}\rho_{\nu}(x),$
(44)

where $\xi$ is a order of unity parameter, which gives the ratio between the
fourth central moment of $\delta_{\sigma}$ and its variance squared, i.e.

$\int_{-\infty}^{+\infty}dz\,z^{4}\delta_{\sigma}(z)\equiv\xi\sigma^{4}.$ (45)

Therefore, Eqs. (36) would be generalised to

$\xi\frac{\sigma^{4}}{48}\rho_{\nu}^{\prime\prime\prime\prime}(x-a)+\frac{\sigma^{2}}{4}\rho_{\nu}^{\prime\prime}(x-a)+\frac{1}{2}\Big{[}\rho_{\nu}(x-a)+\rho_{\nu}(x)\Big{]}=\lambda_{\nu}\rho_{\nu}(x),\quad
x>a/2;$ (46a)
$\xi\frac{\sigma^{4}}{48}\rho_{\nu}^{\prime\prime\prime\prime}(x+a)+\frac{\sigma^{2}}{4}\rho_{\nu}^{\prime\prime}(x+a)+\frac{1}{2}\Big{[}\rho_{\nu}(x+a)+\rho_{\nu}(x)\Big{]}=\lambda_{\nu}\rho_{\nu}(x),\quad
x<-a/2.$ (46b)

Along the same line of reasoning followed before, it is shown that we have
periodic solutions with period $a$. It is only the dispersion relation that
changes, i.e. Eq. (40) holds but now
$\lambda_{\nu}=1-\sigma^{2}k_{\nu}^{2}/4+\xi\sigma^{4}k_{\nu}^{4}/48$. See
also the appendix for a further check of the robustness of the emergence of
periodic solutions.

In order to check the correctness of our theoretical description, we have
numerically diagonalised the master equation, i.e. numerically solved the
eigenproblem (12), for the harmonic potential (26) with different jump
distributions. Again, we have fixed units by taking $k_{B}T=1$ and $\chi=1$.
For each of the families of jump distributions, Eqs. (7)–(9), we have chosen
the parameters so that $\omega(\eta)$ adopts a two-peak structure.

In Fig. 2, the numerical leading eigenfunction $\phi_{1}(x)$ is compared with
the analytical approximation $\phi_{1}^{\text{CDW}}(x)$ in Eq. (40). When it
shows oscillatory behaviour, $\phi_{1}(x)$ is always even: therefore, we have
fitted it to a cosine-shaped eigenfunction,

$\phi_{1}^{\text{fit}}(x)=c(k_{1})\phi_{0}(x)\cos(k_{1}^{*}x),$ (47)

where $c(k_{1})$ is a normalisation constant such that
$\norm{\phi_{1}^{\text{fit}}}=1$, by minimising the distance between the two
functions,

$d(\phi_{1},\phi_{1}^{\text{fit}})\equiv\int_{-\infty}^{+\infty}dx\,\left|\phi_{1}^{\text{fit}}(x)-\phi_{1}(x)\right|^{2}$
(48)

over $k_{1}$. The obtained fits are excellent in all cases, with the largest
separation taking place near the origin, where the amplitude of the numerical
eigenfunction is slightly larger than that of the theoretical curve. This is
reasonable, since the equation for the CDW actually holds for $|x|>a/2$. On
the one hand, the fitted wave vector $k_{1}^{*}=2\pi/a^{*}$ perfectly matches
the theoretical value $k_{1}=2\pi/a$ for the Gaussian and box-like jump
distributions in Eqs. (7) and (9), respectively. On the other hand,
$k_{1}^{*}\neq 2\pi/a$ for the algebraic jump distribution in Eq. (8). We
recall that we described $\delta_{\sigma}(\eta-a)$—which is basically the
quantity shown in the lower right panel of Fig. 2—as a peaked distribution
around $\eta=a$ with a width proportional to its variance. Although this is a
fair description of Gaussian and box-like distributions, it is not so for
algebraic jump distributions, which are not symmetric around the point at
which the probability accumulates as $\sigma\to 0^{+}$. However, a good
theoretical prediction for $k_{1}^{*}$ is always given by
$k_{\text{th}}^{*}=2\pi/a_{\text{th}}^{*}$, with $a_{\text{th}}^{*}$ being the
point at which the cumulative distribution equals $1/2$, i.e.
$\int_{0}^{a_{\text{th}}^{*}}d\eta\,\delta_{\sigma}(\eta-a)=1/2$. For the
parameters employed in the figure, $a_{\text{th}}^{*}=1.58$ while the fitted
value is $a^{*}=1.56$—the difference is approximately $1\%$. Note that
$a^{*}_{\text{th}}=a$ for both the Gaussian and box-like distributions in Eqs.
(7) and (9).

Figure 2: Leading eigenfunction for the Metropolis algorithm. It has been
obtained via numerical diagonalisation of the master equation, for: (upper
left) $\omega^{(a,\sigma)}(\eta)$, with $(a,\sigma)=(1.4,0.1)$, (upper right)
$\omega^{(a,\alpha)}(\eta)$, with $(a,\alpha)=(1.75,6)$, and (lower left)
$\omega^{(a.b)}(\eta)$ with $(a,b)=(1.85,0.15)$. The insets present the ratio
between the leading eigenfunction and the equilibrium distribution
$\phi_{0}(x)$. The dashed curves correspond to the best fit of such ratios to
the function $\cos(k_{1}^{*}x)$, with $k_{1}^{*}=2\pi/a^{*}$, with (upper
left) $a^{*}=1.40$, (upper right), $a^{*}=1.56$, and (lower left)
$a^{*}=1.85$, respectively. The corresponding jump distributions—only for
$\eta>0$, recall that they are even functions of $\eta$—are shown in the lower
right panel: $\omega^{(a,\sigma)}(\eta)$ in red, $\omega^{(a,\alpha)}(\eta)$
in blue, and $\omega^{(a,b)}(\eta)$ in purple.

It is interesting that the leading eigenfunction always seems to be even when
the CDW oscillatory behaviour emerges. This fact hints at the terms neglected
in our derivation of the CDW modes breaking the degeneracy of the sine and
cosine eigenfunctions. A more refined approach would be necessary to
analytically explain the even parity of the leading eigenfunction.

##  IV Phase-diagrams

In this section, we apply our findings in the previous section to investigate
the optimisation of the convergence rate for one family of bi-parametric jump
distributions, and explore the different phases or regimes that emerge when
varying the system parameters. Specifically, we apply our general results to
the specific case of harmonic confinement, as defined in Eq. (26), with the
algebraic family of jump distributions of Eq. (8).

Our goal is to draw the phase diagram of the system in the $(a,\alpha)$ plane
of system parameters. In that way, we will elucidate the interplay between the
diffusion, rejection, and CDW phases. The diffusive phase corresponds to the
region in the $(a,\alpha)$ plane where the leading eigenvalue $\Lambda$ is
given by the diffusive modes—analogously, we speak of the rejection and CDW
phases. As one moves in the $(a,\alpha)$ plane, the jump distribution changes
shape. For distributions that allow only small jumps, the diffusive phase is
expected; for distributions that allow for large jumps, the rejection phase is
expected; for two-peak distributions, the CDW phase is expected. One key
finding of our analysis is that the point $(a^{*},\alpha^{*})$ at which the
leading eigenvalue $\Lambda$ takes its minimum value in the plane is the
triple point of phase coexistence.

In order to characterise the different phases, i.e. diffusion-dominated,
rejection-dominated, and CDW-dominated, we introduce the following two order
parameters, the fidelity $\mathcal{F}$ and the inverse participation ratio
IPR:

  * •

Fidelity: Similarly to that employed in the context of coherent states in
quantum mechanics, we define the fidelity for the problem of our concern as

$\mathcal{F}\equiv\left|\innerproduct{\phi_{1}}{\phi_{1}^{\text{CDW}}}\right|^{2}=\underset{k}{\max}\left|\int_{-\infty}^{+\infty}dx\
c_{1}\phi_{0}(x)\cos(kx)\phi_{1}(x)\right|^{2},$ (49)

which measures the affinity between the actual leading eigenfunction
$\phi_{1}(x)$ of the Metropolis algorithm with the CDW ansatz
$\phi_{1}^{\text{CDW}}$ from Eq. (40) with the corresponding optimal wave
number $k^{*}$. Let us note that the fidelity is closely related to the
distance about the involved functions, $d=2(1-\mathcal{F})$—compare Eq. (49)
with Eq. (48).

The case $\mathcal{F}=0$ corresponds to null affinity, while $\mathcal{F}=1$
results in a perfect match. For $k^{*}\rightarrow 0$, the fidelity tends to
zero, since the actual leading eigenfunction $\phi_{1}(x)$ is orthogonal to
the equilibrium solution.

  * •

Inverse Participation Ratio: This quantity, introduced in Ref. [15], measures
the emergence of the localisation transition in the eigenfunctions of the
master equation. In our numerical diagonalisation of the master equation,
space is discretised into a lattice with $N_{d}$ sites. Within this
discretisation scheme, the eigenfunction $\phi_{\nu}(x)$ simplifies into a set
of $N_{d}$ values $\phi_{\nu}(i)$, $i=1,\ldots,N_{d}$. The IPR for the leading
eigenfunction is

$\text{IPR}\equiv\sum_{i=1}^{N_{d}}|\phi_{1}(i)|^{4}\bigg{/}\left(\sum_{i=1}^{N_{d}}|\phi_{1}(i)|^{2}\right)^{2}.$
(50)

If $\phi_{1}(x)$ is completely delocalised, then $\text{IPR}\rightarrow 0$ for
large $N_{d}$, implying that it has a well-defined continuum limit, with no
delta peaks. If part of the eigenfunction gets localised, then
$\text{IPR}=O(1)$.

In Figs. 3 and 4, the heat maps for the fidelity $\mathcal{F}$, the IPR and
the convergence rate $\Lambda$ are shown on the $(a,\alpha)$ plane for the bi-
parametric family of algebraic jump distributions in Eq. (8). These diagrams
have been obtained by numerically diagonalising the master equation. We recall
that the parameter $a$ gives the typical size of the attempted jumps while the
parameter $\alpha$ controls the width of the distribution, which becomes a
delta peak at $\eta=a$ in the limit as $\alpha\to\infty$.

In Fig. 3, we show the heat map of the order parameters, the fidelity
$\mathcal{F}$ (left) and the IPR (right). First, we analyse the behaviour of
the fidelity. For low values of $\alpha$, the fidelity is always zero. The
jump distribution does not have the double-peak structure, and then the CDW
modes do not emerge. Instead, for $\alpha>\alpha_{c}=1.22$, $\mathcal{F}$
abruptly varies from close-to-zero values to close-to-unity values, when
increasing $a$. For $a<a_{1}(\alpha)$, the Schrödinger modes dominate and we
are in the diffusive regime, for which $\mathcal{F}\approx 0$, while for
$a>a_{1}(\alpha)$, the CDW modes emerge and dominate, and thus
$\mathcal{F}\approx 1$. The fidelity order parameter is discontinuous at the
line $a=a_{1}(\alpha)$, which is thus a first-order transition line. As we
continue increasing $a$, the fidelity begins to smoothly decrease back towards
zero around a certain value $a_{2}(\alpha)$ that is less clearly defined. This
is related to the change from the CDW regime to the rejection regime—as we
will see more clearly when considering the IPR. Second, we analyse the
behaviour of the IPR. For $\alpha<\alpha_{c}$, for which the CDW regime does
not emerge, the localisation transition takes place along the line
$a=a_{3}(\alpha)$. Thereat, the IPR abruptly changes from almost null values
to $O(1)$ ones, and $a=a_{3}(\alpha)$ is another first-order transition
line—as already reported in Ref. [15]. For $\alpha>\alpha_{c}$, the IPR
vanishes both in the regions where the diffusive modes and the CDW modes
dominate. As $a$ increases and approaches the region where the fidelity
smoothly decreases back from unity to zero, the IPR in turn smoothly increases
from zero to unity—and retains non-zero, $O(1)$, values when $a$ is further
increased. The transition line is not as clear-cut as $a_{1}(\alpha)$ and
$a_{3}(\alpha)$, since both $\mathcal{F}$ and the IPR vary smoothly. As a
first approach, we could define the second-order transition line
$a_{2}(\alpha)$ as the locus at which $\mathcal{F}=$IPR. The three so-defined
transition lines cross at the tricritical point
$(a_{c},\alpha_{c})=(2.08,1.22)$.

Figure 3: Heat maps of the fidelity $\mathcal{F}$ (left) and the IPR (right).
Both are plotted in the $(a,\alpha)$ plane of parameters for the algebraic
jump distribution in Eq. (8). Both order parameters have been evaluated by
numerically obtaining the leading eigenfunction of the master equation. The
transition lines between the different phases are plotted with dashed lines
and intersect at the tricritical point $(a_{c},\alpha_{c})=(2.08,1.22)$, which
is represented by a black dot. These lines correspond to (i) $a_{1}(\alpha)$
for $\alpha>\alpha_{c}$ and $a<a_{c}$, (ii) $a_{2}(\alpha)$ for
$\alpha>\alpha_{c}$ and $a>a_{c}$, and (iii) $a_{3}(\alpha)$ for
$\alpha<\alpha_{c}$ and $a>a_{c}$, respectively, as defined in the text.

Figure 4 shows the corresponding heat map of the convergence rate $\Lambda$.
Interestingly, $\Lambda$ attains its global minimum, $\Lambda_{\min}\simeq
0.6170$, at the tricritical point. As is neatly seen, there exists a valley-
like region for $\alpha\leq 2$ in which $\Lambda$ has values that are very
close to the minimum. The valley comprises points close to the transition line
(i) $a_{1}(\alpha)$ between the diffusive and CDW phases, for
$\alpha_{c}<\alpha<2$, and (ii) $a_{3}(\alpha)$ between the diffusive and
rejection phases, for $\alpha<\alpha_{c}$.

Figure 4: Heat map of the convergence rate $\Lambda$ in the $(a,\alpha)$
plane. As in Fig. 3, the dashed curves correspond to the transition lines
between phases. The black point indicates the coordinates of the global
minimum, $(a_{c},\alpha_{c})=(2.08,1.22)$, which coincides with the
tricritical point at which the transition lines intersect.

Figure 5 gives us further information about the phase transitions described
above. The abrupt change of the fidelity across the diffusion-CDW interphase
curve suggests a radical change in the behaviour of the dominant eigenmode of
the Metropolis algorithm. This can be clearly discerned on the left panel of
Fig. 5, where we plot the eigenvalues $\lambda_{\nu}$, $1\leq\nu\leq 9$, as
functions of $a$ for a fixed value of $\alpha=3>\alpha_{c}$. At the transition
point $a_{1}(\alpha)$ between the diffusive and CDW regimes, an eigenvalue
crossing takes place, such that the dominant eigenfunction $\phi_{1}(x)$
changes parity: from odd, for $a<a_{1}$, to even, for $a>a_{1}$. This is
consistent with the behaviour of the dominant eigenfunctions in the respective
phases: in the diffusion regime, the leading eigenfunction $\phi_{1}(x)$ is
provided by the odd sector, whereas, in the CDW regime, the leading
eigenfunctions is even—as discussed in Sec. III.2. Upon further increase in
$a$, the eigenvalues smoothly tend to the rejection curve, marking the onset
of the rejection phase. The right panel shows the corresponding behaviour of
the order parameters as a function of $a$, for the same value of $\alpha$.
First, the transition between the diffusive and the CDW regimes constitutes a
first-order phase transition, characterised by an discontinuous change of
fidelity, which is due to the fact that the takeover of the leading relaxation
mode involves a change of symmetry in the leading eigenfunction. Second, the
transition from the CDW regime to the rejection regime constitutes a higher-
order—at least second-order—phase transition, characterised by the smooth
changes of both $\mathcal{F}$ and IPR, which is due to the also smooth
variation of the leading eigenvalue between the two regimes.

Figure 5: Leading eigenvalues and order parameters as a function of $a$, for a
fixed value of $\alpha=3>\alpha_{c}$. (Left) Eigenvalues $\lambda_{\nu}$,
$\nu=1,\ldots,9$. (Right) Fidelity $\mathcal{F}$ and IPR. In both panels, we
have employed the jump distribution $\omega^{(a,\alpha)}(\eta)$ in Eq. (8).
The purple curve from the left panel corresponds to the maximal rejection
probability, with the light purple region corresponding to the continuum of
eigenvalues obtained upon the localisation transition. The dashed, vertical
lines delimit the three different phases explored in this work: for
$a<a_{1}(\alpha=3)=1.37$, the Schrödinger diffusion-like regime, for
$a_{1}(\alpha=3)=1.37<a<a_{2}(\alpha=3)=2.72$, upon eigenvalue crossing, the
CDW regime, and for $a>a_{2}(\alpha=3)=2.72$, the rejection regime.

##  V Estimating the optimal convergence rate

Given the results previously presented concerning the characterisation of the
CDW regime, in this section, we apply them to estimate the optimal convergence
rate for different families of jump distributions. This is done by considering
different estimators for the convergence rate of the form
$\Lambda_{M}[\omega]$ with $M$ bases, introduced in Eq. (22).

First, we construct an estimator $\Lambda_{4}[\omega]$ for the convergence
rate with $M=4$ bases, by adding the CDW basis to the diffusion-rejection
estimator $\Lambda_{3}[\omega]$, given by Eq. (25), which only included $M=3$
bases. The incorporation of the CDW basis fixes the collapse instability
described at the end of Sec. II.2, i.e. the tendency of the optimal jump
distribution $\omega_{3}^{*}$ for the diffusion-rejection estimator
$\Lambda_{3}[\omega]$ towards a Dirac delta-like distribution, as illustrated
in Fig. 1. That is, the optimal jump distribution $\omega_{4}^{*}$ for
$\Lambda_{4}[\omega]$ no longer presents this issue. Moreover, the estimates
for the optimal convergence rate $\Lambda_{3}^{*}=\Lambda_{3}[\omega_{3}^{*}]$
and $\Lambda_{4}^{*}=\Lambda_{4}[\omega_{4}^{*}]$ are compared. Finally, we
also consider a more complex estimator with $M=6$ bases, with its associated
optimal jump distribution $\omega_{6}^{*}$ and estimate for the optimal
convergence rate $\Lambda_{6}^{*}=\Lambda_{6}[\omega_{6}^{*}]$.

###  V.1 Estimators for the convergence rate

Our starting point is the diffusion-rejection estimator $\Lambda_{3}[\omega]$
introduced in Eq. (25). We build another estimator for the convergence rate by
adding the CDW basis, as given by Eq. (40) thereto, i.e.

$\Lambda_{4}[\omega]\equiv\max\left[\Lambda^{(\text{Sch.})}_{\text{o}}[\omega],\Lambda^{(\text{Sch.})}_{\text{e}}[\omega],\Lambda^{(\text{CDW})}[\omega],\underset{x}{\max}\
R(x;\omega)\right],$ (51)

where $\Lambda^{(\text{CDW})}[\omega]$ is the estimate given by the CDW basis
for the convergence rate.

Let us also consider the localised Gaussian modes

$\phi_{n}(x)=\mathcal{N}_{n}\left[e^{-\frac{(x-a_{n})^{2}}{2\sigma^{2}}}+\epsilon
e^{-\frac{(x+a_{n})^{2}}{2\sigma^{2}}}\right],\quad
a_{n}=\frac{L}{N_{s}}\left(n+\frac{1}{2}\right),\quad\sigma=\frac{L}{N_{s}},\quad
n=0,1,...,N_{s}-1,$ (52)

where $\mathcal{N}_{n}$ is a normalisation constant. The parameter $\epsilon$
is equal to either $+1$ or $-1$ and sets the symmetry of the eigenfunction.
Therefore, this basis can also be split into even and odd sectors. By adding
the even and odd bases of localised Gaussian modes to $\Lambda_{4}$, we build
the more refined estimator with $M=6$ bases:

$\Lambda_{6}[\omega]\equiv\max\left[\Lambda^{(\text{Sch.})}_{\text{o}}[\omega],\Lambda^{(\text{Sch.})}_{\text{e}}[\omega],\Lambda^{(\text{Gau.})}_{\text{e}}[\omega],\Lambda^{(\text{Gau.})}_{\text{o}}[\omega],\Lambda^{(\text{CDW})}[\omega],\underset{x}{\max}\
R(x;\omega)\right],$ (53)

where $\Lambda^{(\text{Gau.})}_{\text{e}}[\omega]$ and
$\Lambda^{(\text{Gau.})}_{\text{o}}[\omega]$ are the estimates for the
convergence rate provided by the even and odd sectors, respectively, of the
localised Gaussian modes basis.

We recall the min-max problem that we solve for estimating the convergence
rate. For each jump distribution $\omega$, the above estimators have a
definite value—the maximum over the $M$ bases considered in each estimator,
which we denote by $\Lambda_{M}[\omega]$. Afterwards, by varying the jump
distribution, we find out the optimal jump distribution that minimises
$\Lambda_{M}[\omega]$: this is what we denote by $\omega_{M}^{*}$, which in
turn provides the estimate $\Lambda_{M}^{*}\equiv\Lambda_{M}[\omega_{M}^{*}]$
for the considered set of bases. Note that we always have
$\Lambda_{6}^{*}\geq\Lambda_{4}^{*}\geq\Lambda_{3}^{*}$; the functional
$\Lambda_{6}[\omega]$ provides the most accurate lower bound for the
convergence rate. The details of the numerical procedure were already given in
Sec. II.2, when discussing the collapse instability.

###  V.2 Fixing the collapse instability

To start with, we show that the collapse instability, which was illustrated in
Fig. 1 for the Gaussian family of jump distributions (7), is fixed by
including the CDW basis. In Fig. 6, we display the optimal jump distributions
obtained via the functionals $\Lambda_{3}$, $\Lambda_{4}$ and $\Lambda_{6}$
for the convergence rate—note that the latter two functionals include the CDW
basis. The cutoff $\sigma_{\text{cutoff}}=0.1$ has been employed when finding
the optimal parameters $(a,\sigma)$ for each estimator. It is neatly seen that
the variances of $\omega_{4}^{*}(\eta)$ and $\omega_{6}^{*}(\eta)$ no longer
reach the cutoff, so the collapse instability is fixed. Moreover,
$\omega_{4}^{*}$ and $\omega_{6}^{*}$ are close and quite different from the
Dirac-like distribution $\omega_{3}^{*}(\eta)$. On the one hand, all the
functionals give estimates for the convergence rate that are not so different,
$\Lambda_{3}^{*}=0.6114$, $\Lambda_{4}^{*}=0.6138$, and
$\Lambda_{6}^{*}=0.6155$. On the other hand, the incomplete convergence of
$\omega_{3}^{*}$ can further be verified by inserting it into
$\Lambda_{4}[\omega]$ and $\Lambda_{6}[\omega]$. This gives
$\Lambda_{4}[\omega_{3}^{*}]=0.9239$ and $\Lambda_{6}[\omega_{3}^{*}]=0.9521$,
which have a relative error of more than $50\%$ with respect to our best
estimate of the convergence rate $\Lambda_{6}^{*}$. If we do a similar check
with $\omega_{4}^{*}$, by evaluating $\Lambda_{6}[\omega_{4}^{*}]=0.6239$, the
relative error is reduced to $1.36\%$, thus clearly showing the improvement of
the convergence obtained with the incorporation of the CDW basis.

Figure 6: Optimal jump distributions for the Gaussian family in Eq. (7). These
optimal jump distributions have been obtained by employing the three
functionals $\Lambda_{3}[\omega]$, $\Lambda_{4}[\omega]$, and
$\Lambda_{6}[\omega]$ for the convergence rate. Specifically, we show
$\omega_{6}^{*}$ (purple), with associated parameters
$(a_{6}^{*},\sigma_{6}^{*})=(1.45,0.50)$ for $\Lambda_{6}[\omega]$,
$\omega_{4}^{*}$ (red), with $(a_{4}^{*},\sigma_{4}^{*})=(1.43,0.42)$ for
$\Lambda_{4}[\omega]$, and $\omega_{3}^{*}$ (blue), with
$(a_{3}^{*},\sigma_{3}^{*})=(1.38,0.10)$ for $\Lambda_{3}[\omega]$. In the
legend, we present the values of the convergence rate mentioned in the main
text.

We further support the above picture by applying the same methodology to the
other two bi-parametric families of jump distributions: (i) the algebraic
$(a,\alpha)$ one, Eq. (8), and (ii) the box-like $(a,b)$ one, Eq. (9). We
recall that the algebraic and the box-like families tend to a double Dirac-
delta peak in the limit $\alpha\to\infty$ and $b\to 0^{+}$, respectively. This
implies that we have to impose cutoffs in $\alpha$ and $b$, respectively, to
avoid convergence problems: specifically, we have employed the values
$\alpha_{\text{cutoff}}=4$ and $b_{\text{cutoff}}=0.2$. In Fig 7, we show
again that the collapse instability arises when considering the diffusion-
rejection estimator $\Lambda_{3}[\omega]$, and that this instability is fixed
by incorporating the CDW basis, i.e. by considering, at least,
$\Lambda_{4}[\omega]$ 555Of course, more complex functionals like
$\Lambda_{6}[\omega]$, which also include the CDW basis, also fix the
instability and improve the estimate for the convergence rate.. When using the
diffusion-rejection estimator $\Lambda_{3}[\omega]$, the optimal parameters in
the distributions $\omega_{3}^{*}$ reach the corresponding cutoffs, thus
showing the collapse instability. Going to the estimators
$\Lambda_{4}[\omega]$ or $\Lambda_{6}[\omega]$ solves this issue again. In
fact, the global picture is completely analogous to the one found for the
family of Gaussian distributions: (i) in all cases, $\Lambda_{3}^{*}$,
$\Lambda_{4}^{*}$ and $\Lambda_{6}^{*}$ are quite close—of course, verifying
always $\Lambda_{3}^{*}<\Lambda_{4}^{*}<\Lambda_{6}^{*}$, (ii) as compared
with $\Lambda_{6}^{*}$, $\Lambda_{6}[\omega_{3}^{*}]$ has a large error that
is radically diminished by $\Lambda_{6}[\omega_{4}^{*}]$.

Figure 7: Optimal jump distributions for the algebraic and box-like family of
distributions. Specifically, the algebraic (8) and box-like (9) cases are
considered in the left and right panels, respectively. Similarly to Fig. 6,
the optimal jump distributions depicted correspond to the functionals
$\Lambda_{3}[\omega]$, $\Lambda_{4}[\omega]$, and $\Lambda_{6}[\omega]$ for
the convergence rate. Specifically, we show: $\omega_{6}^{*}$(purple) , for
which $(a_{6}^{*},\alpha_{6}^{*})=(2.08,1.22)$ and
$(a_{6}^{*},b_{6}^{*})=(1.44,0.99)$, $\omega_{4}^{*}$ (red), for which
$(a_{4}^{*},\alpha_{4}^{*})=(1.98,1.57)$ and
$(a_{4}^{*},b_{4}^{*})=(1.43,0.63)$ and (blue) $\omega_{3}^{*}$, for which
$(a_{3}^{*},\alpha_{3}^{*})=(1.67,4.00)$ and
$(a_{3}^{*},b_{3}^{*})=(1.30,0.20)$. Again, the values of the convergence rate
mentioned in the main text are given in the legend.

##  VI Further optimising the convergence rate

In the previous sections, we have carried out the minimisation of the
convergence rate over several bi-parametric families of jump distributions,
all of which tend to a double Dirac-delta peak in a certain limit. Here, we
carry out this minimisation over more complex jump distributions. In this way,
we try to shed some light on the global optimisation problem, i.e. given a
confining potential $U(x)$, look for the jump distribution $\omega^{*}(\eta)$
that minimises the convergence rate of the dynamics.

###  VI.1 More complex jump distributions

In order to see if we can further improve our estimation of the global minimum
value of the convergence rate, we consider a more complex, discrete-point,
family of the form:

$\omega(\eta)=\sum_{j=1}^{N_{\text{P}}}a_{j}\,\delta(|\eta|-\eta_{j}),$ (54)

where both the weights $a_{j}$ and the positions $\eta_{j}$ of each of the
$N_{\text{P}}$ delta peaks are adjustable parameters, over which we carry out
the optimisation procedure. Figure 8 shows the optimal jump distribution for
the case $N_{\text{P}}=20$, employing the $\Lambda_{6}[\omega]$ functional for
the convergence rate. With this approach, we get $\Lambda_{6}^{*}=0.6146$,
which is slightly smaller than the previously found values for the two-
parametric families in Sec. V.2. But, more interestingly, although most of the
optimal jump distribution is concentrated in a definite interval,
$0.5\lesssim\eta\lesssim 2.0$, there are a couple of delta peaks for larger
jumps, around $\eta=3$. This behaviour hints at the actual optimal jump
distribution having a more complicated structure, with at least four
peaks—recall that $\omega(\eta)$ is an even function of $\eta$, so the two
peaks found for $\eta>0$ are replicated for $\eta<0$.

Figure 8:  Optimal jump distribution for the discrete-point family of Eq.
(54). Specifically, we have taken $N_{\text{P}}=20$. The vertical dashed lines
give the value of the coefficients $a_{j}$ of each of the delta peaks at
discrete positions $\eta_{j}$. The optimal jump distribution has been obtained
using the $\Lambda_{6}$ functional of the convergence rate, and the
corresponding estimate $\Lambda_{6}^{*}$ for the convergence rate is shown in
the legend.

To check the above insights, we have employed the following two-Gaussian
family of jump distributions with a cutoff,

$\omega(\eta)=C\sum_{n=1}^{2}\alpha_{n}\
\text{exp}\left(-\frac{(|\eta|-a_{n})^{2}}{2\sigma_{n}^{2}}\right)\theta(\eta_{\max}-|\eta|),\quad
a_{1}<a_{2},$ (55)

where $C$ is a normalisation constant. Note that this family has seven
parameters: $(a_{1},\mu_{1},\sigma_{1},a_{2},\mu_{2},\sigma_{2},\eta_{\max})$.
The constants $\alpha_{i}>0$ give the relative weights of the four peaks at
$\eta=\pm a_{1}$ and $\eta=\pm a_{2}$. In addition, this family includes as a
particular case the bi-parametric Gaussian family in Eq. (7), in the limit
$\eta_{\max}\to\infty$ and (e.g.) $\alpha_{1}\to 0^{+}$.

Figure 9: Optimal jump distribution for the 7-parameter family in Eq. (55).
Similarly to the previous figures, the optimal jump distributions have been
obtained using three different functionals of the convergence rate:
$\Lambda_{6}^{*}$ (purple), $\Lambda_{4}^{*}$ (red), and $\Lambda_{3}^{*}$
(blue). In the legend, we present the values of the convergence rate mentioned
in the text. The parameters characterising each jump distribution are shown in
Table 1.

Figure 9 shows the optimal jump distributions for the 7-parameter family (55),
once more employing the functionals $\Lambda_{M}[\omega]$ with $M=3$, $4$, and
$6$. Again, the CDW fixes the collapse instability shown by the estimate
$\omega_{3}^{*}$ of the diffusion-rejection estimator—note that $\sigma_{1}$
again has reached the cutoff value $\sigma_{\text{cutoff}}=0.1$.
Interestingly, this 7-parameter family slightly improves the optimal
convergence rate: here we have $\Lambda_{6}^{*}=0.6133$, which is smaller than
the values for the previous bi-parametric families, displayed in Figs. 6-7.
Apart from the small decrease in the optimal convergence rate, the optimal
jump distributions $\omega_{4}^{*}$ and $\omega_{6}^{*}$ have a four-peak
structure, analogous to that found in Fig. 8 for the discrete-point
family—recall that we only plot $\omega(\eta)$ for $\eta>0$, due to its even
parity.

The above analysis points out that the global optimal jump distribution may
present at least four peaks, corresponding to jumps centred around $\pm a_{1}$
and $\pm a_{2}$, instead of the two peaks at $\pm a$ we have considered in the
previous sections.

|  $\eta_{\max}$ |  $\alpha_{1}$ |  $a_{1}$ |  $\sigma_{1}$ |  $\alpha_{2}$ |  $a_{2}$ |  $\sigma_{2}$
---|---|---|---|---|---|---|---
$\omega_{6}^{*}$ |  $3.358$ |  $0.985$ |  $1.247$ |  $0.428$ |  $0.986$ |  $3.845$ |  $0.546$
$\omega_{4}^{*}$ |  $3.265$ |  $0.893$ |  $1.198$ |  $0.354$ |  $1.050$ |  $3.343$ |  $0.235$
$\omega_{3}^{*}$ |  $5.764$ |  $1.000$ |  $1.144$ |  $0.100$ |  $0.100$ |  $0.074$ |  $2.950$
Table 1: Parameters characterising the optimal jump distributions for the
7-parameter family in Eq. (55). As in previous cases, optimisation is carried
out for the functionals $\Lambda_{M}[\omega]$, with $M=3$, $4$ and $6$.

###  VI.2 Sequence of jumps

The fact that four-peak jump distributions improve the convergence rate as
compared with two-peak jump distributions, uncovered in Sec. VI.1, also
provides insights into the design of other strategies for further optimising
the convergence rate. As an example of such a possible strategy, we consider
the deterministic alternation between two jump distributions
$\omega_{1}(\eta)$ and $\omega_{2}(\eta)$ with peaks at $\eta=\pm a_{1}$ and
$\eta=\pm a_{2}$, respectively.

The corresponding kernel $\hat{\mathcal{K}}(\omega_{1},\omega_{2})$ for the
master equation of this modified Metropolis algorithm is

$\hat{\mathcal{K}}(\omega_{1},\omega_{2})=\hat{K}(\omega_{1})\hat{K}(\omega_{2})\hat{K}(\omega_{2})\hat{K}(\omega_{1}),$
(56)

with $\hat{K}(\omega_{i})$ being the kernel from our original Metropolis
algorithm from Eq. (11), with a fixed choice $\omega_{i}(\eta)$ for the jump
distribution. We are employing parenthesis instead of the brackets that
appeared in Eq. (11), in order to highlight that $\omega$ is no longer a
variational parameter. This new kernel describes a modified Metropolis
algorithm in which two jumps from $\omega_{1}$ are attempted, then two from
$\omega_{2}$ are attempted, and the process repeats. The sequence have been
chosen so that the kernel form in Eq. (56) remains symmetric.

In the following, we numerically compare the efficiency of this strategy
against our original Metropolis algorithm with average the jump distribution
$\omega=(\omega_{1}+\omega_{2})/2$. Now, as a consequence of the structure of
the sequential kernel in Eq. (56), the corresponding convergence rate for the
modified Metropolis algorithm is
$\Lambda_{\text{seq}}=(\Lambda_{1}\Lambda_{2})^{1/2}$, with $\Lambda_{i}$
being the convergence rate for the jump distribution $\omega_{i}$. As can be
clearly observed, the sequence of jumps improves the convergence time of the
original algorithm from $\Lambda_{6}^{*}=0.6133$ to
$\Lambda_{\text{seq}}=0.5840$, i.e. an improvement of $4.78\%$, thus
corroborating our initial intuition. We remark that this improval is just for
a reduced number of degrees of freedom: in this case, one standard deviation
and two mean values for the Gaussian ansatz. Increasing the dimension of the
parameter space, or even combining different families of jump
distributions—such as the ones employed before, given by Eqs. (7), (8) and
(9)—should further improve the minimisation of the convergence rate.

Figure 10: Comparison of optimal jump distributions. Specifically, we show the
optimal jump distributions for (i) the 7-parametric jump distribution in Eq.
(55), already shown in Fig. 9 (purple), and (red) the sequence of jumps
described in the text, alternating two Gaussian distributions with the same
standard deviation but different mean value. The corresponding optimal
convergence rates mentioned in the text are also plotted in the legend.

We may understand the better performance of the sequence-of-jumps strategy
over the original Metropolis scheme with the following intuitive argument. Let
us take the characteristic lengths $a_{1}<a_{2}$, with a sufficient gap
between the two. Then, for $\omega_{1}(\eta)$, the first two steps of the
algorithm will resemble that of free diffusion, which allows exploration of
small regions of the phase space. In order to move further, we take the second
two steps, which allow one to jump to further points in space, and thus the
rejection rate increases accordingly. Thus, each distribution balances the
problems presented for the other.

The previous argument has remarkable implications in terms of choosing between
stochastic and deterministic jump strategies. On the one hand, the former
would correspond to our original Metropolis algorithm: if we consider the
four-peak jump distribution $\omega(\eta)$ as the average of two independent
distributions $\omega_{1}(\eta)$ and $\omega_{2}(\eta)$ with different
characteristic lengths, $\omega=(\omega_{1}+\omega_{2})/2$, we may understand
the algorithm as choosing randomly between these two distributions before
attempting any jumps. On the other hand, the deterministic approach would
correspond to the sequence of jumps, since we clearly set the order in which
we draw the attempted jumps from $\omega_{1}$ and $\omega_{2}$.

##  VII Conclusions

In this paper, we have analysed the relaxation dynamics of the Metropolis
algorithm for a particle confined in a one-dimensional potential, with even
jump distributions that ensure that the dynamics verify detailed balance. The
central result of this paper is the functional optimisation of the convergence
rate of the Metropolis algorithm in one dimension, which is obtained as a
balance of three mechanisms: diffusion, rejection, and the novel CDW regime
that we unveil and characterise in this paper. In the CDW regime, the leading
eigenfunction $\phi_{1}(x)$ corresponding to the largest eigenvalue
$\lambda_{1}$—except for the equilibrium one $\lambda_{0}=1$—shows an
oscillatory behaviour, characterised by a wave number $k$. For two-parametric
jump distributions, the balance among diffusion, rejection, and CDW occurs at
the tricritical point of coexistence—at which the minimum value of the
convergence rate is attained.

The diffusion-rejection functional of Ref. [15]—$\Lambda_{3}[\omega]$ in the
notation employed in this paper—leads to the collapse instability: for two-
peaked jump distributions, the diffusion-rejection functional gives an optimal
jump distribution with a Dirac delta shape that breaks ergodicity in one-
dimensional systems. This collapse instability is fixed by adding the CDW
basis to the above functional, thus building the new functional
$\Lambda_{4}[\omega]$ that includes the diffusion, rejection, and CDW bases.
The overall take-home message is that optimality, in the sense of minimising
the convergence rate, is attained as a balance between the three physical
factors above: diffusion, rejection, and CDW.

This interplay between physical factors has been illustrated by the phase
diagram for the convergence rate. Specifically, this phase diagram has been
built for the algebraic $(a,\alpha)$ family of jump distributions of Eq. (8).
The different phases are characterised by two order parameters: fidelity
$\mathcal{F}$, which is different from zero when the new CDW regime dominates,
and IPR (inverse participation ratio), which is different from zero when the
rejection regime dominates. First, for smooth jump distributions without a
double-peak structure, we have the localisation transition between the
diffusion and rejection regions as the length of the typical jumps is
increased—which was characterised in Ref. [15]. Then, we have unveiled a new
transition between the diffusion and CDW regimes in the parameter range for
which the jump distribution approaches a double peak structure, also as the
length of the typical jumps is increased. The diffusion-CDW transition is
first-order, so that (i) $\mathcal{F}$ and the optimal wave number abruptly
change, and (ii) the parity of the leading eigenfunction of the master
equation changes—from odd to even. Finally, it is only when the length of the
typical jumps is further increased that we find a transition between the CDW
and rejection regimes. At this transition, both order parameters vary
smoothly, so the CDW-rejection transition is at least second order.

Furthermore, we have shown numerical evidence on the improvement of the
convergence rate by combining different jump distributions. The sequence of
four jumps exposed in this work is just an example of how the interplay
between more than one characteristic length may improve the relaxation to
equilibrium. Interestingly, a similar improvement of other figures of merit
when considering more than one characteristic length has been observed in
other physical systems—e.g. for optimal search strategies and resetting
processes [25, 26, 27].

The optimisation of the convergence rate can be viewed as engineering the jump
distribution to accelerate the dynamics of the system. Thus our approach
shares the objective of the general field of shortcuts, started in the quantum
field of shortcuts to adiabaticity [28] and later transposed to classical and
stochastic systems as swift state-to-state transformations [29]. Therein, to
accelerate the relaxation to equilibrium, one may engineer the potential
acting on the system; see e.g. [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
41]. Note that engineering the jump distribution is analogous to engineering
the potential, since changing the shape of the jump distribution varies the
jump moments—which provide the coefficients in the Kramers-Moyal expansion of
the master equation [8, 9].

In this paper, we have completely characterised the CDW regime in one
dimension, providing a physical interpretation thereof in terms of breaking of
ergodicity. Thus, the natural question is whether or not the CDW regime would
play a role in the global optimisation of the Metropolis algorithm in higher
dimensions. Also, it is relevant to investigate the optimisation of the
convergence rate for other dynamics, different from the Metropolis
algorithm—like Glauber dynamics [42]. We would like to stress that our
analysis has focused on the leading eigenvalue $\lambda_{1}$, which as a rule
controls the convergence rate to equilibrium. However, there are certain
specific situations in which the convergence rate to equilibrium is controlled
by the following eigenvalue $\lambda_{2}<\lambda_{1}$, because the coefficient
of the first mode vanishes, as in the strong version of the Mpemba effect [43,
44, 45, 46, 47, 48, 49, 50] or when specific strategies are employed to tune
the dynamical evolution of the system [51, 52]. Therefore, it is also worth
investigating whether the general picture found here for $\lambda_{1}$, i.e.,
the competition between diffusion, rejection, and CDW regimes, extends to
$\lambda_{2}$ or even to the full eigenvalue spectrum of the master equation.
Certainly, all these questions open the door to future research.

##  VIII Acknowledgements

We thank Satya Majumdar for fruitful discussions. A. Patrón and A. Prados
acknowledge financial support from Grant PID2021-122588NB-I00 funded by
MCIN/AEI/ 10.13039/501100011033/ and by “ERDF A way of making Europe”. A.
Patrón, A. Prados, and E. Trizac acknowledge financial support from Grant
ProyExcel_00796 funded by Junta de Andalucía’s PAIDI 2020 programme. A. Patrón
also acknowledges support from the FPU programme through Grant FPU2019-4110,
and additional support from the FPU programme through Grant EST22/00346, which
funded his research stay at Univ. Paris-Saclay during autumn 2022. A. Prados
also acknowledges the hospitality of LPTMS, which funded his stay at Univ.
Paris-Saclay in June 2022.

##  Appendix A Ergodicity breaking and periodic solutions

In this section, we provide more details on the emergence of periodic
solutions in the limit $\sigma\to 0^{+}$, analysed in Sec. III of the main
text. In this limit, the jump distribution reduces to the double Dirac-delta
peak in Eq. (10). To be consistent with the discussion on the robustness of
the CDW modes from Sec. III.2, we restrict ourselves to even potentials
$U(x)$, although a similar analysis can be carried out for a non-symmetric
potential.

For $\sigma=0$, the eigenvalue/eigenfunction equation for the modes of the
master equation (15) reduces exactly to

$\displaystyle e^{-\beta|U(x)-U(x-a)|/2}\phi_{\nu}(x-a)$
$\displaystyle+e^{-\beta|U(x)-U(x+a)|/2}\phi_{\nu}(x+a)+\left[1-e^{-\beta(U(x-a)-U(x))}\right]\theta(U(x-a)-U(x))$
$\displaystyle+\left[1-e^{-\beta(U(x+a)-U(x))}\right]\theta(U(x+a)-U(x))=2\lambda_{\nu}\phi_{\nu}(x),\quad\forall
x.$ (57)

It is clearer to split the above equation into three intervals, (i) $x>a/2$,
(ii) $-a/2<x<a/2$, and (iii) $x<-a/2$. In this way, $U(x)-U(x\pm a)$ has a
definite sign on each interval, and we get

$\displaystyle e^{\beta U(x-a)/2}\phi_{\nu}(x-a)+e^{-\beta(U(x+a)-U(x))}$
$\displaystyle\left[e^{\beta U(x+a)/2}\phi_{\nu}(x+a)-e^{\beta
U(x)/2}\phi_{\nu}(x)\right]$ $\displaystyle=(2\lambda_{\nu}-1)e^{\beta
U(x)/2}\phi_{\nu}(x),$ $\displaystyle x>a/2,$ (58a) $\displaystyle
e^{-\beta(U(x+a)-U(x))}\left[e^{\beta U(x+a)/2}\phi_{\nu}(x+a)\right.$
$\displaystyle\left.-e^{\beta U(x)/2}\phi_{\nu}(x)\right]$ $\displaystyle+\
e^{-\beta(U(x-a)-U(x))}\left[e^{\beta U(x-a)/2}\phi_{\nu}(x-a)\right.$
$\displaystyle\left.-e^{\beta U(x)/2}\phi_{\nu}(x)\right]$
$\displaystyle=2(\lambda_{\nu}-1)e^{\beta U(x)/2}\phi_{\nu}(x),$
$\displaystyle-a/2<x<a/2,$ (58b) $\displaystyle e^{\beta
U(x+a)/2}\phi_{\nu}(x+a)+e^{-\beta(U(x-a)-U(x))}$ $\displaystyle\left[e^{\beta
U(x-a)/2}\phi_{\nu}(x-a)-e^{\beta U(x)/2}\phi_{\nu}(x)\right]$
$\displaystyle=(2\lambda_{\nu}-1)e^{\beta U(x)/2}\phi_{\nu}(x),$
$\displaystyle x<-a/2.$ (58c)

Indeed, the structure of this equation suggests the change of variables in Eq.
(29). Introducing it into Eq. (58),

$\displaystyle\rho_{\nu}(x-a)+e^{-\beta(U(x+a)-U(x))}\left[\rho_{\nu}(x+a)\right.\left.-\rho_{\nu}(x)\right]$
$\displaystyle=(2\lambda_{\nu}-1)\rho_{\nu}(x),$ $\displaystyle x>a/2,$ (59a)
$\displaystyle
e^{-\beta(U(x+a)-U(x))}\left[\rho_{\nu}(x+a)-\rho_{\nu}(x)\right]+e^{-\beta(U(x-a)-U(x))}[\rho_{\nu}(x-a)-$
$\displaystyle\rho_{\nu}(x)]$ $\displaystyle=2(\lambda_{\nu}-1)\rho_{\nu}(x),$
$\displaystyle-a/2<x<a/2,$ (59b)
$\displaystyle\rho_{\nu}(x+a)+e^{-\beta(U(x-a)-U(x))}\left[\rho_{\nu}(x-a)\right.\left.-\rho_{\nu}(x)\right]$
$\displaystyle=(2\lambda_{\nu}-1)\rho_{\nu}(x),$ $\displaystyle x<-a/2.$ (59c)

It is worth stressing that Eq. (59) for $\rho_{\nu}(x)$ is exact for
$\sigma=0$. Interestingly, any periodic function $\rho_{\nu}(x)$ with period
$a$, i.e. $\rho_{\nu}(x\pm a)=\rho_{\nu}(x)$ $\forall x$, satisfies this
equation with eigenvalue $\lambda_{\nu}=1$. As already discussed in the main
text, this is intimately related to the breaking of ergodicity of the one-
dimensional Monte Carlo algorithm for $\sigma=0$: any point $x$ is only
connected to $x_{n}=x\pm na$, $n\in\mathbb{N}$, and therefore the probability
density of this lattice becomes distributed following the Boltzmann weight,
i.e. it is proportional to $\exp(-\beta U(x_{n}))$, but the lattices
corresponding to different values of $x$ are disconnected and the system does
not reach equilibrium.

The consideration of a small width $\sigma$ makes it possible to introduce an
expansion in powers of $\sigma^{2}$ around Eq. (59), as carried out in the
main text—see Eqs. (36) and (46). Therein, we showed that these additional
terms still allow for the emergence of periodic solutions with period $a$, but
restrict their shape to the cosine and sine CDW modes in Eq. (39). Note that
the shape of the CDW modes so obtained does not depend on the potential
$U(x)$, due to our assumption (32) that allowed us to neglect the terms
proportional to $e^{-\beta(U(x\pm a)-U(x))}$ in Eqs. (59a) and (59c), and
analogous terms in the higher-order terms in our expansion in powers of
$\sigma^{2}$.

## References

  * Frenkel and Smit [2023] D. Frenkel and B. Smit,  _Understanding Molecular Simulation: From Algorithms to Applications_ (Elsevier Science, 2023).
  * Newman and Barkema [1999] M. Newman and G. Barkema,  _Monte Carlo Methods in Statistical Physics_ (Clarendon Press, 1999).
  * Mode [2011] C. J. Mode,  _Applications of Monte Carlo Methods in Biology, Medicine and Other Fields of Science_ (IntechOpen, Rijeka, 2011).
  * Glasserman [2004] P. Glasserman,  _Monte Carlo Methods in Financial Engineering_, Applications of mathematics: stochastic modelling and applied probability (Springer, 2004).
  * Bishop [2006] C. Bishop,  _Pattern Recognition and Machine Learning_ (Springer, 2006).
  * Gilks  _et al._ [1995] W. Gilks, S. Richardson, and D. Spiegelhalter,  _Markov Chain Monte Carlo in Practice_, Chapman & Hall/CRC Interdisciplinary Statistics (Taylor & Francis, 1995).
  * Metropolis  _et al._ [1953] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, Equation of state calculations by fast computing machines, The Journal of Chemical Physics 21, 1087 (1953).
  * Gardiner [2009] C. Gardiner,  _Stochastic Methods: A Handbook for the Natural and Social Sciences_ , Springer Series in Synergetics (Springer Berlin Heidelberg, 2009).
  * Van Kampen [1992] N. G. Van Kampen,  _Stochastic processes in Physics and Chemistry_ (North-Holland, 1992).
  * Gelman  _et al._ [1996] A. Gelman, G. O. Roberts, and W. R. Gilks, Efficient Metropolis Jumping Rules, in  _Bayesian Statistics 5: Proceedings of the Fifth Valencia International Meeting_ (Oxford University Press, 1996) _eprint: https://academic.oup.com/book/0/chapter/422210114/chapter-pdf/52447340/isbn-9780198523567-book-part-38.pdf.
  * Gelman  _et al._ [1997] A. Gelman, W. R. Gilks, and G. O. Roberts, Weak convergence and optimal scaling of random walk Metropolis algorithms, The Annals of Applied Probability 7, 110 (1997).
  * Schmon and Gagnon [2022] S. M. Schmon and P. Gagnon, Optimal scaling of random walk Metropolis algorithms using Bayesian large-sample asymptotics, Statistics and Computing 32, 28 (2022).
  * Krauth [2006] W. Krauth,  _Statistical Mechanics: Algorithms and Computations_, Oxford Master Series in Physics (Oxford University Press, UK, 2006).
  * Talbot  _et al._ [2003] J. Talbot, G. Tarjus, and P. Viot, Optimum Monte Carlo simulations: some exact results, Journal of Physics A: Mathematical and General 36, 9009 (2003).
  * Chepelianskii  _et al._ [2023] A. D. Chepelianskii, S. N. Majumdar, H. Schawe, and E. Trizac, Metropolis Monte Carlo sampling: convergence, localization transition and optimality, Journal of Statistical Mechanics: Theory and Experiment 2023, 123205 (2023).
  * Grüner and Zettl [1985] G. Grüner and A. Zettl, Charge density wave conduction: A novel collective transport phenomenon in solids, Physics Reports 119, 117 (1985).
  * Grüner [1988] G. Grüner, The dynamics of charge-density waves, Reviews of Modern Physics 60, 1129 (1988).
  * Chen  _et al._ [2016] C.-W. Chen, J. Choe, and E. Morosan, Charge density waves in strongly correlated electron systems, Reports on Progress in Physics 79, 084505 (2016).
  * Note [1] To simplify our notation, we do not show the explicit dependence on the temperature. In simple cases like the harmonic potential, the temperature can be absorbed in the definition of the unit of length.
  * Bapat and Raghavan [1997] R. B. Bapat and T. E. S. Raghavan,  _Nonnegative Matrices and Applications_, Encyclopedia of Mathematics and its Applications (Cambridge University Press, 1997).
  * Note [2] The following numerical procedure has been employed for all the examples presented throughout the paper. To obtain the eigenvalues and eigenfuntions of the master equation, we have employed the FORTRAN 90 subroutine dsyev.f90 from LAPACK [53]. To avoid spurious boundary effects, the spatial variable $x$ is restricted to a finite interval $[-L,L]$, with a sufficiently large $L$. Moreover, space is discretised with a mesh of $N_{d}$ points. For the harmonic potential, we have fixed units by choosing $k_{B}T=1$ and $\chi=1$, and taken $L=10$ and $N_{d}=2000$. To numerically estimate $\Lambda$, as defined in Eq. (25), we have employed the open source nonlinear optimisation libraries from NLopt [54] with $N_{s}=100$ modes for each basis.
  * Note [3] Since the numerics involve a spatial discretisation, $\sigma$ must be larger than the lattice spacing in this discretisation.
  * Note [4] A similar issue was found in Ref. [15] for the Schrödinger eigenvalues.
  * Note [5] Of course, more complex functionals like $\Lambda_{6}[\omega]$, which also include the CDW basis, also fix the instability and improve the estimate for the convergence rate.
  * Bénichou  _et al._ [2011] O. Bénichou, C. Loverdo, M. Moreau, and R. Voituriez, Intermittent search strategies, Reviews of Modern Physics 83, 81 (2011).
  * Evans  _et al._ [2020] M. R. Evans, S. N. Majumdar, and G. Schehr, Stochastic resetting and applications, Journal of Physics A: Mathematical and Theoretical 53, 193001 (2020).
  * García-Valladares  _et al._ [2023] G. García-Valladares, C. A. Plata, A. Prados, and A. Manacorda, Optimal resetting strategies for search processes in heterogeneous environments, New Journal of Physics 25, 113031 (2023).
  * Guéry-Odelin  _et al._ [2019] D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Shortcuts to adiabaticity: Concepts, methods, and applications, Reviews of Modern Physics 91, 045001 (2019).
  * Guéry-Odelin  _et al._ [2023] D. Guéry-Odelin, C. Jarzynski, C. A. Plata, A. Prados, and E. Trizac, Driving rapidly while remaining in control: classical shortcuts from Hamiltonian to stochastic dynamics, Reports on Progress in Physics 86, 035902 (2023).
  * Aurell  _et al._ [2011] E. Aurell, C. Mejía-Monasterio, and P. Muratore-Ginanneschi, Optimal Protocols and Optimal Transport in Stochastic Thermodynamics, Physical Review Letters 106, 250601 (2011).
  * Deffner and Lutz [2013] S. Deffner and E. Lutz, Quantum Speed Limit for Non-Markovian Dynamics, Physical Review Letters 111, 010402 (2013).
  * Jarzynski [2013] C. Jarzynski, Generating shortcuts to adiabaticity in quantum and classical dynamics, Physical Review A 88, 040101 (2013).
  * Tu [2014] Z. C. Tu, Stochastic heat engine with the consideration of inertial effects and shortcuts to adiabaticity, Physical Review E 89, 052148 (2014).
  * Martínez  _et al._ [2016] I. A. Martínez, E. Roldán, L. Dinis, D. Petrov, J. M. R. Parrondo, and R. A. Rica, Brownian Carnot engine, Nature Physics 12, 67 (2016).
  * Li  _et al._ [2017] G. Li, H. T. Quan, and Z. C. Tu, Shortcuts to isothermality and nonequilibrium work relations, Physical Review E 96, 012144 (2017).
  * Patra and Jarzynski [2017] A. Patra and C. Jarzynski, Shortcuts to adiabaticity using flow fields, New Journal of Physics 19, 125009 (2017).
  * Plata  _et al._ [2020] C. A. Plata, D. Guéry-Odelin, E. Trizac, and A. Prados, Finite-time adiabatic processes: Derivation and speed limit, Physical Review E 101, 032129 (2020).
  * Zhang [2020] Y. Zhang, Work needed to drive a thermodynamic system between two distributions, EPL (Europhysics Letters) 128, 30002 (2020).
  * Funo  _et al._ [2020] K. Funo, N. Lambert, F. Nori, and C. Flindt, Shortcuts to Adiabatic Pumping in Classical Stochastic Systems, Physical Review Letters 124, 150603 (2020).
  * Plata  _et al._ [2021] C. A. Plata, A. Prados, E. Trizac, and D. Guéry-Odelin, Taming the Time Evolution in Overdamped Systems: Shortcuts Elaborated from Fast-Forward and Time-Reversed Protocols, Physical Review Letters 127, 190605 (2021).
  * Ilker  _et al._ [2022] E. Ilker, O. Güngör, B. Kuznets-Speck, J. Chiel, S. Deffner, and M. Hinczewski, Shortcuts in Stochastic Systems and Control of Biophysical Processes, Physical Review X 12, 021048 (2022).
  * Glauber [1963] R. J. Glauber, Time-dependent statistics of the Ising model, Journal of Mathematical Physics 4, 294 (1963).
  * Lu and Raz [2017] Z. Lu and O. Raz, Nonequilibrium thermodynamics of the Markovian Mpemba effect and its inverse, Proceedings of the National Academy of Sciences 114, 5083 (2017).
  * Klich  _et al._ [2019] I. Klich, O. Raz, O. Hirschberg, and M. Vucelja, Mpemba Index and Anomalous Relaxation, Physical Review X 9, 021060 (2019).
  * Kumar and Bechhoefer [2020] A. Kumar and J. Bechhoefer, Exponentially faster cooling in a colloidal system, Nature 584, 64 (2020).
  * Busiello  _et al._ [2021] D. M. Busiello, D. Gupta, and A. Maritan, Inducing and optimizing Markovian Mpemba effect with stochastic reset, New Journal of Physics 23, 103012 (2021).
  * Kumar  _et al._ [2022] A. Kumar, R. Chétrite, and J. Bechhoefer, Anomalous heating in a colloidal system, Proceedings of the National Academy of Sciences 119, e2118484119 (2022).
  * Zhang and Hou [2022] S. Zhang and J.-X. Hou, Theoretical model for the Mpemba effect through the canonical first-order phase transition, Physical Review E 106, 034131 (2022).
  * Biswas  _et al._ [2023] A. Biswas, R. Rajesh, and A. Pal, Mpemba effect in a Langevin system: Population statistics, metastability, and other exact results, The Journal of Chemical Physics 159, 044120 (2023).
  * Patrón  _et al._ [2023] A. Patrón, B. Sánchez-Rey, C. A. Plata, and A. Prados, Non-equilibrium memory effects: Granular fluids and beyond, Europhysics Letters 143, 61002 (2023).
  * Gal and Raz [2020] A. Gal and O. Raz, Precooling Strategy Allows Exponentially Faster Heating, Physical Review Letters 124, 060602 (2020).
  * Chittari and Lu [2023] S. S. Chittari and Z. Lu, Geometric approach to nonequilibrium hasty shortcuts, The Journal of Chemical Physics 159, 084106 (2023).
  * Anderson  _et al._ [1999] E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen,  _LAPACK Users’ Guide_, 3rd ed. (Society for Industrial and Applied Mathematics, 1999).
  * Johnson [2007] S. G. Johnson, The NLopt nonlinear-optimization package, https://github.com/stevengj/nlopt (2007).