# A robust and memory-efficient transition state search method for complex
energy landscapes

Samuel J. Avis  Jack R. Panter<EMAIL_ADDRESS>Halim Kusumaatmaja
<EMAIL_ADDRESS>Department of Physics, Durham University,
South Road, Durham DH1 3LE, UK

###### Abstract

Locating transition states is crucial for investigating transition mechanisms
in wide-ranging phenomena, from atomistic to macroscale systems. Existing
methods, however, can struggle in problems with a large number of degrees of
freedom, on-the-fly adaptive remeshing and coarse-graining, and energy
landscapes that are locally flat or discontinuous. To resolve these
challenges, we introduce a new double-ended method, the Binary-Image
Transition State Search (BITSS). It uses just two states that converge to the
transition state, resulting in a fast, flexible, and memory-efficient method.
We also show it is more robust compared to existing bracketing methods that
use only two states. We demonstrate its versatility by applying BITSS to three
very different classes of problems: Lennard-Jones clusters, shell buckling,
and multiphase phase-field models.

##  I Introduction

Transition states are central to the description of reconfiguration mechanisms
for systems in chemistry, condensed matter physics, and engineering.
Historically, many computational methods for locating transition states have
grown from an atomistic or particulate perspective. These have proven to be
important tools for understanding, for example, protein folding Bryngelson
_et al._ (1995); Onuchic, Luthey-Schulten, and Wolynes (1997), biological and
industrial catalysis Boehr  _et al._ (2006); Kerns  _et al._ (2015); Guo  _et
al._ (2018), quantum tunnelling Richardson  _et al._ (2016); Vaillant,
Althorpe, and Wales (2019), crystallisation Richard and Speck (2018), and
cluster formation Wales, Miller, and Walsh (1998); Wales (2012).

More recently, it is increasingly being recognised that transition states are
useful in mesoscale or macroscale systems. Here, the minimum energy barriers
provide important lower bounds to the energy input required for transitions to
occur. This has been used to understand failure in structural engineering
applications Panter  _et al._ (2019); Hutchinson and Thompson (2018), for the
development of super liquid-repellent surfaces Zhang and Ren (2014); Panter,
Gizaw, and Kusumaatmaja (2019); Amabili  _et al._ (2017), and investigating
locomotion through complex terrain for robotics Othayoth, Thoms, and Li
(2020). Moreover, it is becoming desirable to tailor elastic deformation
transitions to enable technologies such as advanced deployable structures
Filipov  _et al._ (2015); Zhai, Wang, and Jiang (2018), mechanical sensors and
actuators Bertoldi  _et al._ (2017); Truby and Lewis (2016); Chi  _et al._
(2022); Bonfanti  _et al._ (2020), and energy absorbers Shan  _et al._ (2015);
Giri and Mailen (2021) to name but a few.

Transition state search methods generally fall into two categories, single-
and double-ended methods. Single-ended methods are initialised at a single
state and attempt to climb to a nearby saddle point. Examples include
eigenvector following Cerjan and Miller (1981), the dimer method Heyden, Bell,
and Keil (2005); Kästner and Sherwood (2008); Zhang, Du, and Zheng (2016), and
climbing image methods E, Ren, and Vanden-Eijnden (2007); Ren and Vanden-
Eijnden (2013). Double-ended methods can be further subdivided into two
groups. The first utilise a chain of states between two minima which is then
minimised to provide an estimate for the full transition pathway in addition
to the transition state. Examples are the string method E, Ren, and Vanden-
Eijnden (2002, 2007) and doubly-nudged elastic band (DNEB) Trygubenko and
Wales (2004). These methods require an appropriate initial interpolation,
which can sometimes be challenging to obtain Wales and Carr (2012). The second
group are bracketing methods, which involve two states converging to the
transition state from either side. These include the Dewar-Healy-Stewart (DHS)
algorithm Dewar, Healy, and Stewart (1984), ridge method Ionova and Carter
(1993), the step and slide method Miron and Fichthorn (2001), and the double-
ended surface walking method Zhang, Shang, and Liu (2013).

A large range of landscapes, however, prove challenging or impossible to
explore via these methods. One key problem arises from the push towards larger
and more complex systems Trefethen and Thiyagalingam (2013); Shalf (2020);
Alexander  _et al._ (2020), resulting in the need to develop algorithms that
are more computationally and memory efficient, and can incorporate
optimisation strategies such as on-the-fly adaptive remeshing and coarse-
graining. These typically involve changing the resolution or discretisation of
the systems to focus the computational time on important regions, such as
using a higher resolution mesh in regions of high stress in finite element
simulations Lee and Bathe (1994). However, chain-of-states methods involve a
coupling between the configurations of each state, and so there is an issue if
they have different discretisations and numbers of degrees of freedom.
Meanwhile, single-ended methods can be inefficient because they are not well
suited for identifying specific pathways and can spend a large amount of time
searching for undesired transition states. Another major challenge in studying
complex energy landscapes relates to the presence of locally flat or
discontinuous regions, such as when considering patchy McMullen  _et al._
(2018); Nguemaha and Zhou (2018); Chen, Zhang, and Torquato (2018) and hard-
body Richard and Speck (2018); Santra, Singh, and Bagchi (2018) interactions
in atomistic simulations, systems of polymer chains Mökkönen, Ala-Nissila, and
Jónsson (2016), or collision constraints for macroscopic objects Wriggers
(2006). Flat zero-modes in the landscape pose issues for single-ended search
methods and current bracketing methods that rely only upon local information.
Specialist treatment can sometimes be used such as in the case of global
rotation and translation Page and McIver (1988), but they are thwarted by
local zero-modes. Finally, current methods cannot typically be applied in the
case of discontinuous potentials, or if the gradient is prohibitively
expensive to compute, because continuous, differentiable optimisation
functions are required.

In this work we introduce a new double-ended bracketing method, the Binary-
Image Transition State Search (BITSS). Using a range of different
applications, we demonstrate that it successfully addresses each of the above
challenges. In addition, we show that BITSS is superior compared to existing
bracketing methods, allowing us to access the transition states when other
methods fail.

##  II BITSS method

Figure 1:  Schematics of the BITSS method on a simple 2D potential with two
minima. (a) The orange line shows the trajectories of the two states from the
minima (blue) to the transition states (red) under the BITSS method. The
minimum energy pathway is shown by the dashed line. (b) A snapshot of the
BITSS minimisation showing the driving forces on each state due to the energy
constraint, $\bm{F}_{\mathrm{E}}$, and distance constraint,
$\bm{F}_{\mathrm{D}}$, with $E_{1}<E_{2}$ and
$d(\bm{x}_{1},\bm{x}_{2})<d_{i}$. (c) The final configuration of the BITSS
method showing the two states in orange, the transition state in red, and the
negative curvature eigenvector, $\bm{\hat{\tau}}$.

The method begins by first initialising the states, $\bm{x}_{1}$ and
$\bm{x}_{2}$, in the basins of attraction of different local minima, such as
the two blue spots in the 2d potential in fig. 1a. These can be set to the
minima, but this is not a necessary requirement. The energies of these two
states are then minimised, while constraining their separation. This is
iteratively reduced to zero, such that, at iteration $i$, their separation is

$d_{i}=(1-f)d_{i-1},$ (1)

with $d_{0}$ taking the value of the separation between the two initial
states. A reduction factor of $f=0.5$ is successful for most applications, but
this can be made smaller to ensure that the states do not slide off the ridge
between the two basins of attraction. Different metrics may be used to compute
this distance, although in this work we simply use the Euclidean distance,

$d(\bm{x}_{1},\bm{x}_{2})=\sqrt{\sum_{i}(x_{1,i}-x_{2,i})^{2}}.$ (2)

To further ensure that neither state is pulled over the ridge, a secondary
constraint enforces equal energies for the two states. Using this strategy,
the two states will meet at the lowest point on the ridge, the transition
state.

The two constraints are implemented using energy penalty terms, which result
in driving forces on the two states if the constraints are not met, such as in
fig. 1b. Including these energy penalty terms gives the total BITSS energy for
the pair of states,

$E_{\mathrm{BITSS}}(\bm{x}_{1},\bm{x}_{2})=E_{1}+E_{2}+\kappa_{e}\left(E_{1}-E_{2}\right)^{2}\\\
+\kappa_{d}\left(d(\bm{x}_{1},\bm{x}_{2})-d_{i}\right)^{2},$ (3)

where $E_{1}$ and $E_{2}$ are the single-state energies, and $\kappa_{e}$ and
$\kappa_{d}$ parametrise the strengths of the energy and distance constraints.

In this work the L-BFGS algorithm is chosen to minimise this energy, owing to
its fast convergence and low memory requirement for large numbers of degrees
of freedom Liu and Nocedal (1989). However, any other minimisation method can
be used instead.

To ensure that the transition state is located successfully, the constraint
strengths $\kappa_{d}$ and $\kappa_{e}$ are updated as the algorithm proceeds
using information from the system. These are set such that the driving forces
due to the constraints and single-state energies are of similar size. This
prevents the constraints from dominating the underlying potential or causing
large jumps that make a state pass over the ridge. This results in the
following equations (see Supplementary Note I for the derivation),

$\displaystyle\kappa_{e}=\frac{\alpha}{2E_{\mathrm{B}}},$ (4)
$\displaystyle\kappa_{d}=\max\left(\frac{\sqrt{\left|\bm{\nabla}E_{1}\right|^{2}+\left|\bm{\nabla}E_{2}\right|^{2}}}{2\sqrt{2}\beta
d_{i}}\;,\;\frac{E_{\mathrm{B}}}{\beta d_{i}^{2}}\right),$ (5)

where $\bm{\nabla}E_{1}$ and $\bm{\nabla}E_{2}$ are the gradients of the
energies of the two states, and $\alpha$ and $\beta$ are parameters with
recommended values of $\alpha=10$ and $\beta=0.1$. Here, $E_{\mathrm{B}}$ is
an estimation for the current energy barrier, evaluated using the difference
between the highest energy along a linear interpolation between the two states
and the average energy of the two states. These constraints are initially
calculated at the start of each minimisation, and regularly recalculated
throughout (once per 100 iterations is used in this work).

In practice, when numerically minimising, the states will jump about slightly
which can result in large gradients perpendicular to the optimal movement
direction. To reduce this effect, the gradients used in eq. 5 are projected in
the direction of the separation between the two states:

$\left|\bm{\nabla}E_{n}\right|\approx\frac{\left|(\bm{x}_{1}-\bm{x}_{2})\cdot\bm{\nabla}E_{n}\right|}{\left|\bm{x}_{1}-\bm{x}_{2}\right|}.$
(6)

In summary, the method involves iteratively performing the following three
steps:

  1. 1.

Reduce the constrained separation, according to eq. 1.

  2. 2.

Minimise the potential of the pair of states, eq. 3.

  3. 3.

Recompute the constraint coefficients, $\kappa_{e}$ and $\kappa_{d}$, at
regular intervals using eqs. 4 and 5.

This process is completed once a suitable convergence criterion is reached.
This can either be based upon the separation between the states, the size of
gradient at the midpoint between them, or the change in the position of the
midpoint.

Using the BITSS approach, the typical trajectories of the states are
demonstrated for a simple 2D potential in fig. 1a. Initially, the lower energy
state jumps up to satisfy the equal energy constraint and then moves to
minimise the separation without increasing its energy. Then, the two states
converge directly towards one another, before being deflected towards the
saddle in the ridge. Consequently, if there are multiple possible pathways
between two states, BITSS will be biased towards identifying those that are
more direct or with lower energy. Furthermore, the final two states are
positioned either side of the transition state in the direction of the
negative curvature eigenvector, $\bm{\hat{\tau}}$ (fig. 1c). So, BITSS
automatically identifies the ‘reactive mode’ and associated eigenvalue in
addition to the transition state. Once the transition state has been
identified, it is possible to find the full minimum energy pathway by tracing
the trajectory of downhill minimisations from the two final states, which are
either side of the saddle.

In the event that there are intermediate stable states, there will be a chain
of multiple transition states between the two minima. In this case, the equal-
energy constraint will not prevent the states from passing over the lower
energy transition states, so BITSS should converge to the transition state
with the highest energy. This enables the identification of the overall energy
barrier, providing estimates about the overall ease of the transition, or the
rate for chemical processes. However, as demonstrated in Supplementary Note
II, if multiple transition states have very similar energies then a smaller
distance reduction factor, $f$, may be necessary to ensure that it does indeed
converge to the highest transition state. Furthermore, if all of transition
states or the full pathway are desired, BITSS can be continually repeated from
one of the minima downhill from the located transition state and one of the
initial minima until the initial minima are piecewise connected by a full
pathway.

##  III Results and Discussion

###  III.1 Comparison with other bracketing methods

Figure 2:  Trajectories of the bracketing methods on a hooked potential with a
single saddle point. The equation for this potential is provided in section
C.1. The minimum energy pathway is shown by the dashed line.

The BITSS potential in eq. 3 and iterative steps above offer key advantages
over existing bracketing methods that also use two states to locate the
transition state. For instance, in the ridge method Ionova and Carter (1993),
the two images are initially chosen to bracket the largest energy point on an
interpolated path between the two endpoints. However, this is not guaranteed
to be on the ridge containing the transition state, and specialist methods are
required to avoid high-energy local maxima, or when the initial path contains
multiple candidate maxima. In another example, the double ended surface
walking method Zhang, Shang, and Liu (2013) requires Gaussian bias potentials
to be added at each iteration to force two dimers to climb uphill in the
landscape. For high numbers of degrees of freedom and many iterations, this
becomes very computationally expensive.

The two methods most similar to BITSS are the DHS Dewar, Healy, and Stewart
(1984) and step and slide Miron and Fichthorn (2001) methods. In the step and
slide method, the separation between two images is minimised while their
energy is fixed (iteratively increasing the energy up to the transition
state). Conversely in the DHS method, the energy of an image is minimised
while the image separation is fixed (iteratively decreasing the separation and
changing the frozen image up to the transition state). To illustrate how BITSS
is superior compared to these methods, we consider the hooked 2d potential in
fig. 2. For this potential the energies of the images ascend higher than that
of the transition state and consequently both of these methods fail to
converge to the saddle point regardless of the parameters that are used. The
step and slide method fails in this situation because it always expects that
the energy of the two states is below the saddle point if they have not
converged, so it has no means of descending down the ridge. For DHS, the
images reach a certain point at which one state can pass over the ridge by
minimising its energy. At this point DHS will fail even if the distance is
reduced very slowly. In contrast, the BITSS method is successful for this
potential. This is because the combination of distance and energy constraints
allow BITSS to approach a transition state from both below and above (by
sliding down a ridge).

Furthermore, using both an energy and distance constraint with BITSS provides
improved efficiency over these methods which each use just one of the
constraints. In the case of DHS, fixing one state in place and optimising the
other means that the amount that the separation is reduced must be much
smaller than BITSS to ensure that it does not pass over the ridge. Meanwhile,
in step and slide, it is difficult to obtain a reasonable energy increment
when the two states are far from the transition state, leading to a larger
than necessary number of iterations. Also, restricting the minimisation to a
constant energy surface can result in a considerably more complex method, as
the states must be constantly projected back onto this surface.

###  III.2 Comparison with chain-of-states methods

In contrast to BITSS, chain-of-states methods do not typically find transition
states directly. Instead, they are designed to find the full pathway (or an
approximation thereof), and a secondary method can then be used to refine to
the transition state. As we will demonstrate, this strategy is successful for
simple, linear pathways, but faces two key challenges when the pathway is
highly non-linear. Firstly, for such complex pathways, a large number of
states are required to sufficiently approximate the minimum energy pathway.
The second is that choosing a suitable initial interpolation can be
problematic to achieve. BITSS can be advantageous in both these regards, as
only two states are evolved, regardless of the pathway complexity, and no
initial interpolation is required.

Here we compare the speeds of convergence of BITSS to two widely used
approaches for finding transition states that employ chain-of-states methods:
climbing image nudged elastic band (CINEB) Henkelman, Uberuaga, and Jónsson
(2000), and DNEB with hybrid eigenvector following (DNEB-HEVF) Cerjan and
Miller (1981). The core of these methods involves minimising the total energy
of a chain of states, connected by elastic springs to keep them equally spaced
along the transition pathway. We fix the two end-points at the minima, so the
number of states that are minimised is two fewer than the number of states in
the chain. CINEB modifies the method by altering the behaviour of the state
with the highest energy. The direction of minimisation on this state is
inverted along the pathway direction, effectively converting the saddle point
into a local minimum. Alternatively, DNEB-HEVF involves first minimising the
chain of states until a convergence criterion is met, and then performing
hybrid eigenvector following from the highest energy state, moving uphill
along the smallest eigenvector of the Hessian until it reaches the transition
state. For completeness, we also combine hybrid eigenvector following with
BITSS and include the results in the convergence comparison. Additional
implementation details for these methods are included in appendix D.

It is also possible to use the string method with a climbing image E, Ren, and
Vanden-Eijnden (2007) or eigenvector following Zimmerman (2013); although, in
this case, the results are expected to be similar to the nudged elastic band
methods. We note that our aim in this section is to observe how the BITSS
method behaves for different systems, rather than providing a comprehensive
comparison of the currently available methods, which has been performed in
other works Koslover and Wales (2007); Sheppard, Terrell, and Henkelman
(2008).

Figure 3:  The three systems used for comparison with chain-of-states methods.
These are: (a) a Lennard-Jones seven-particle cluster, (b) cylindrical shell
buckling, and (c) wetting of a chemically-striped surface. The configurations
shown correspond to the two minimum energy states and the transition state,
marked by an asterisk.

Three diverse systems are used for this comparison, exhibiting a broad range
of energy landscapes. The first system is a two-dimensional, seven-particle
cluster, interacting via a Lennard-Jones pair potential. This is a frequently
used test system for studying transition rates Wales (2002); Passerone and
Parrinello (2001). Here, the 14 degrees of freedom are the particle
coordinates. The characteristic transition shown in fig. 3a sees a particle
rearrangement between two close-packed clusters.

The second system is an elastic cylindrical shell, modelled by a triangulated
mesh of nodes, which interact via extensional and angular springs. The $35400$
degrees of freedom are the node coordinates in three-dimensional space. The
characteristic transition in fig. 3b shows the formation of a stable dimple
from an initially unbuckled cylinder. This transition is essential to capture
and predict mechanical failure under strain Panter  _et al._ (2019); Virot
_et al._ (2017).

System | BITSS |  BITSS- HEVF | DNEB-HEVF | CINEB
---|---|---|---|---
3 | 5 | 10 | 20 | 3 | 5 | 10 | 20
LJ-7 | 148 | 144 | 135 | 138 | 208 | 361 | 30 | 153 | 2000 | 1692
Buckling | $12866$ | $14446$ | $73694$ | $7315$ | $6242$ | $12094$ | — | $8352$ | $19360$ | $77904$
Wetting | $12100$ | $17262$ | $17840$ | $17721$ | $20221$ | $48286$ | — | — | — | —
Table 1:  Number of potential gradient calculations required to reach the
transition state for the three comparison examples. The climbing image nudged
elastic band (CINEB) and DNEB with hybrid eigenvector following (DNEB-HEVF)
methods have been run for different numbers of images. Convergence is
determined to be when the root-mean-square of the gradient at the estimate for
the transition state is less than $10^{-4}$. The fields left blank indicate
that the method has not converged to the correct transition state.

The final system involves a droplet situated on a chemically striped surface
with both hydrophilic and hydrophobic regions. Droplet transitions on
patterned surfaces such as this are vital to understand as powerful bio-
inspired liquid manipulation strategies Kusumaatmaja  _et al._ (2006); Brown
and Bhushan (2016). In the example shown in fig. 3c, a droplet transitions
from two hydrophilic patches to one patch. Here, the system is represented by
a diffuse-interface model, in which the $40000$ degrees of freedom are the
local fluid compositions at each site of the discretised domain. The pathway
for this example is highly non-linear in the coordinate space because each
degree of freedom only varies when it is at the interface of the droplet. As a
result, the initial pathway for the chain-of-states methods cannot be a simple
linear interpolation. Instead, the position of a semi-circular droplet is
interpolated between the two final positions.

The results for the three systems are shown in table 1. First, we note that
for all three systems, using hybrid eigenvector following does not
significantly improve the speed of BITSS. Indeed, for buckling and especially
wetting, HEVF is detrimental to performance. Next, it is interesting to
compare each method’s performance between simple and complex pathways. In
contrast to the wetting transition’s highly non-linear pathway, the pathways
of the LJ-7 rearrangement and the buckling system can be simply tracked
following a gradual variation in the order parameters. For LJ-7, this is the
translation of atoms 5 and 6, and for buckling, this is the radial
displacement of the centre of the dimple Panter  _et al._ (2019). For the
simpler pathways, BITSS is generally slower, but for the complex pathways,
BITSS is faster. Moreover, we see that for the wetting example, CINEB does not
converge to the transition state because the estimated tangent vector is
highly inaccurate due to the non-linearity of the pathway.

For situations where memory is limited, it is important to minimise the number
of images used. However, efficiently finding the TS is challenging for both
CINEB and DNEB-HEVF if too few images are used, as observed for the
cylindrical buckling with three images. BITSS, on the other hand, converges
using only two images.

###  III.3 Adaptive discretisation

Figure 4:  Demonstration of BITSS addressing the challenges associated with
adaptive remeshing. (a) Snapshots of the BITSS method for the buckling of a
cylinder with a changing mesh. The radial displacement relative to the
unbuckled cylinder is shown, as well as the underlying triangular mesh. (b)
Snapshots for the striped wetting example with different resolutions for the
two states. Each grid cell denotes 50x50 lattice nodes. The zoomed axis shows
the difference in the fluid interface between the two final states, as well as
the approximated transition state (solid black line). This is compared to the
transition state found using a high resolution (dashed line).

Adaptive remeshing and coarse-graining are widely used techniques that we can
utilise to further increase the efficiency of BITSS. These techniques cause
issues for most existing double-ended methods because the coupled states may
end up with different degrees of freedom. However, in BITSS the only direct
coupling is in the distance measure, $d(\bm{x}_{1},\bm{x}_{2})$, which is
relatively easy to adapt. Here we demonstrate the use of adaptive remeshing by
considering two separate issues.

Firstly, we show in fig. 4a that BITSS is able to handle the discretisation
adapting, and the number of degrees of freedom changing, as the method runs.
For this we use the cylindrical buckling example with the resolution
increasing from 40 to 100 triangles around the cylinder, corresponding to an
increase from $1760$ to $11000$ degrees of freedom. This demonstrates that
BITSS is able to converge to the transition state so long as the remeshing is
not so significant as to shift a state into the basin of attraction of the
other minimum.

In the second test, shown in fig. 4b, we demonstrate the use of different
meshes for the two states in the striped wetting example. In this case, the
distance measure is adapted by interpolating one state onto the other mesh and
computing the Euclidean distance. However, for some applications a simpler
measure may be sufficient, such as the difference between average values of
the system. Using this approach, BITSS is able to closely approach the
transition state. The precision of this convergence is now limited by the
transition state energy differing slightly on each grid, but this effect will
be reduced when using an adaptive method or a higher resolution.

###  III.4 Complex landscapes

Figure 5:  Demonstration of BITSS applied to flat and discontinuous
potentials. (a) Energy profile of the BITSS pathway on a 2D potential with
flat regions. Blue and red dots denote the minima and transition state,
respectively. Points of interest are labelled by i–iii (see text). Top inset:
A zoomed in view around the transition state. Bottom inset: The pathway taken,
with the edges of the flat regions marked by dashed lines. (b) The
discontinuous hard-core pair potential used in the seven-particle cluster
(orange). The standard Lennard-Jones potential is also shown in grey. (c)
Disconnectivity graphs of the energy landscapes for the seven-particle cluster
with the two potentials. The two graphs are offset for visibility.

The final challenges we will address are those related to complex landscapes
that prove challenging for previous algorithms. The first is the presence of
flat regions in the landscape. Fig. 5a shows BITSS applied to a 2D landscape
with two such regions (i & ii) that are flat in the $x$-direction. We see
BITSS is able to successfully converge past these flat regions, even with one
very close to the transition state (ii). In these regions there are no driving
forces due to the potential and the energy constraint, which use purely local
information about the gradient. However, the distance constraint continues to
pull the states together, preventing them from getting stuck. When only a
single state has a zero-gradient mode then the other is likely to slide down
the potential slightly (iii), but the two states still remain either side of
the dividing ridge and so the result is unaffected.

An additional consideration is the case where the potential energy surface is
flat at the top of the pathway. There are two possibilities here, one is that
the potential is flat in a direction perpendicular to the tangent of the
pathway, such that the ridge is level. In this case BITSS is unaffected and it
will be able to converge to some point along the ridge. An example of this is
the free global rotation and translation of the Lennard-Jones cluster in fig.
3a. The other possibility is that the flat mode is in the direction of the
pathway. In this case there is no single transition state along the pathway,
but instead a region. BITSS would be ill suited in this situation because the
equal-energy constraint would not prevent the images from passing over the
saddle and falling to a minimum.

Finally, we investigate the application of BITSS to systems with undefined
gradients, such as when the landscape is discontinuous. To account for this,
the equations for the coefficients must be adapted to not depend upon the
gradients, and a gradient-free minimiser (simulated annealing) is used. These
changes are detailed in the methods section. This has been tested using a
7-particle cluster with a hard-core Lennard-Jones pair-potential, shown in
fig. 5b, which results in a discontinuous landscape. Using the gradient-free
approach, BITSS is able to successfully find the transition states, allowing
us to plot the disconnectivity graph of the system, shown in fig. 5c. Compared
with the results for the standard Lennard-Jones cluster, the energies of the
minima are largely unchanged, but the energies of the transition states are
found to be slightly higher. This indicates that the particles in the Lennard-
Jones cluster cut the corner slightly as they transition, whereas this is not
possible using the discontinuous potential, resulting in higher energies.
Despite this gradient-free method being feasible, it is worth noting that a
gradient-based approach is significantly more efficient, and so should be
preferred if possible.

##  IV Conclusion

Overall, we have developed the binary image transition state search (BITSS)
algorithm for the efficient location of transition states in traditionally
challenging landscapes. This has distinct advantages for complex pathways
owing to the lack of a required initial pathway estimate, as well as the
identification of the transition state that provides the overall energy
barrier in multi-step pathways. From the speed analysis, we find that the
combination of chain-of-states methods with single-ended transition state
search methods provides good performance for near-linear pathways, such as for
the Lennard-Jones cluster and cylindrical buckling. However, for highly
complex and non-linear pathways, as exhibited by the striped wetting example,
BITSS is superior. Indeed, the demonstrated speed and memory-efficiency will
be key as we move towards studying larger and more complex systems using
BITSS.

A second source of efficiency in the BITSS method comes from the ability to
adaptively change the degrees of freedom as the algorithm proceeds. We
demonstrated how transition states could be found by both increasing the
resolution upon convergence, and coupling systems with different
discretisations. The ease of coupling two copies of a system and adaptive
remeshing, now leads to the possibility of incorporating BITSS into existing
open-source optimisation methods, such as surface evolver Brakke (1992) or
finite element methods Kolev  _et al._ (2021), to provide important energy
barrier functionality.

Finally, we showed how BITSS can be used to survey discontinuous energy
landscapes, demonstrated for a system of attractive hard-core particles. This
opens up possibilities for studying a broad range of systems previously out of
reach of conventional landscape methods, but where transition information is
valuable. These include systems with very short range interactions, such as in
colloidal clusters, or hard contact forces, such as in the folding of elastic
materials, or locomotion and environmental interaction in robotics.

The distance metric between the two BITSS images is interesting to analyse
further. One question that emerges is whether transition states can be located
by coupling two images through a small number of collective properties, rather
than the total distance between all degrees of freedom in the system. A second
question concerns landscapes with multiple competing pathways between states.
In such cases, it may be possible to access transition states different from
the most direct one by using a biased distance metric. A further investigation
that is now open to pursue is when discontinuities in the landscape occur at
‘stationary points’ (now properly referred to as critical points). In this
case, a transition state can no longer be defined by its Hessian eigenvalues,
but instead is more broadly defined as a region of locally minimal energy that
separates two basins of attraction to minima. Overall, it will be interesting
to explore how BITSS enables access to even more challenging landscapes, and
those not yet amenable to traditional landscape exploration techniques.

## Supplementary material

See the supplementary material for derivations of the expressions of the
constraint coefficients in eqs. 5 and 4, and the demonstration of the BITSS
method applied to a path with multiple transition states.

###### Acknowledgements.

S. J. A. is supported by a studentship from the Engineering and Physical
Sciences Research Council [Grant No. EP/R513039/1]. H. K. and J. R. P.
acknowledge funding from the Engineering and Physical Sciences Research
Council [Grant No. EP/V034154/1].

## Conflict of Interest

The authors have no conflicts to disclose.

## Author Contributions

All authors conceived and planned the project. S. J. A. performed the
numerical experiments and analysed the data. J. R. P. and H. K. supervised the
study. S. J. A. wrote the paper, with review and input from and J. R. P. and
H. K.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A BITSS: Changes for undefined gradients

A couple of alterations to the method must be made to account for situations
where the gradients are unknown. First, the calculation of $\kappa_{d}$ in eq.
5 must be adapted to avoid the use of gradients. This can be done by simply
removing the first term and just using the second term in the equation.
Secondly, L-BFGS can no longer be used because it requires knowledge of the
gradients. We must instead use a minimiser that does not require a
differentiable optimisation function, for which we use simulated annealing
Kirkpatrick, Gelatt, and Vecchi (1983). This has a chance of randomly jumping
one state over the dividing barrier, but we can reduce this probability by
limiting the initial temperature and maximum random displacement. We typically
employ $T_{0}=E_{\mathrm{B}}/10$, and
$d_{\mathrm{max}}=d(\bm{x}_{1},\bm{x}_{2})/100$.

##  Appendix B Adaptive discretisation test details

Here we provide the details for the interpolations and mapping involved in the
two examples demonstrating the feasibility of using an adaptive discretisation
method. For the cylindrical buckling example with a changing mesh, the
resolution is refined each time the separation between the two states has
halved, and is performed at the end of each iteration of the BITSS method.
This involves the number of triangles around the circumference of the cylinder
increasing along the sequence: $40\rightarrow 60\rightarrow 80\rightarrow
100$; with the number of degrees of freedom increasing by:
$$1760$\rightarrow$3960$\rightarrow$7040$\rightarrow$11000$$. The positions of
the nodes on the new grid, $\\{\bm{n}_{i}\\}$, are determined by linear
interpolation from the previous grid $\\{\bm{p}_{i}\\}$, using the positions
of the unbuckled meshes, $\\{\bm{n^{\prime}}_{i}\\}$ and
$\\{\bm{p^{\prime}}_{i}\\}$. For each node of the new grid,
$\bm{n^{\prime}}_{i}$, the triangle that contains it is first identified,
which we will denote
$\\{\bm{p^{\prime}}_{1},\bm{p^{\prime}}_{2},\bm{p^{\prime}}_{3}\\}$, and the
barycentric coordinates of the point are computed,
$\\{\lambda_{1},\lambda_{2},\lambda_{3}\\}$. The new position is then given by
$\bm{n}_{i}=\lambda_{1}\bm{p}_{1}+\lambda_{2}\bm{p}_{2}+\lambda_{3}\bm{p}_{3}$.

In the wetting example with different resolutions for the two states, the
distance is obtained by first mapping the phase field from the higher
resolution grid, $\\{\phi_{k,l}|k,l\in\\{0,1,\cdots,399\\}\\}$, to the low
resolution grid, $\\{\phi^{\prime}_{i,j}|i,j\in\\{0,1,\cdots,199\\}\\}$.
Because a square grid is used with a resolution ratio of two, the mapping
involves averaging each 2x2 block to a single point:

$\phi^{\prime}_{i,j}=\frac{1}{4}\left(\phi_{2i,2j}+\phi_{2i+1,2j}+\phi_{2i,2j+1}+\phi_{2i+1,2j+1}\right).$
(7)

Then the separation from the other state, $\\{\widetilde{\phi}_{i,j}\\}$, is
computed using the 2-norm,

$d=\sqrt{\sum_{i,j}\left(\phi^{\prime}_{i,j}-\widetilde{\phi}_{i,j}\right)^{2}}.$
(8)

Finally, the gradient of the distance with respect to each point must be
mapped back to the higher-resolution grid, which is done by assigning a
quarter of each component back to its the original four points:

$\frac{\partial
d}{\partial\phi_{k,l}}=\frac{\partial\phi^{\prime}_{i,j}}{\partial\phi_{k,l}}\frac{\partial
d}{\partial\phi^{\prime}_{i,j}}=\frac{1}{4}\frac{\partial
d}{\partial\phi^{\prime}_{i,j}}=\frac{1}{4}\frac{\phi^{\prime}_{i,j}-\widetilde{\phi}_{i,j}}{d},$
(9)

where $k\in\\{2i,2i+1\\}$, and $l\in\\{2j,2j+1\\}$.

##  Appendix C Energy and gradient expressions for the example systems

###  C.1 2D potentials

The 2D potentials in figs. 1 and 2 use a sum of Gaussian potentials, eq. 10,
the parameters for which are provided in tables 2 and 3.

$E(x,y)=\sum_{i}a_{i}\exp\\!\left(-\frac{(x-b_{x,i})^{2}}{c_{x,i}}-\frac{(y-b_{y,i})^{2}}{c_{y,i}}\right)$
(10)

$a$ | $b_{x}$ | $b_{y}$ | $c_{x}$ | $c_{y}$
---|---|---|---|---
-3 | -1.4 | 0 | 1 | 1
-2 | 1.4 | 0 | 1 | 1
-1 | 0.07 | 1 | 1 | 1
Table 2:  Parameters for the Gaussians to produce the potential used in fig. 1.  $a$ | $b_{x}$ | $b_{y}$ | $c_{x}$ | $c_{y}$
---|---|---|---|---
-1 | 0 | 0 | 10 | 10
1 | 0 | 0 | 1 | 1
5 | 2 | 0 | 1 | 0.1
-1 | 1 | 1 | 0.1 | 0.1
-1 | 1 | -1 | 0.1 | 0.1
0.01 | 0 | 0 | 1 | 1
0.5 | -2 | 0 | 1 | 1
Table 3:  Parameters for the Gaussians to produce the potential used in fig.
2.

###  C.2 Particle cluster system

In the particle cluster example system, the Lennard-Jones potential is used
for the interaction between each pair of particles. Therefore, the potential
and its gradient for each pair of particles are

$\displaystyle
E=4\epsilon\left[\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right],$
(11) $\displaystyle\frac{\partial E}{\partial\bm{x}_{1}}=-\frac{\partial
E}{\partial\bm{x}_{2}}=\frac{24\epsilon(\bm{x}_{2}-\bm{x}_{1})}{r^{2}}\left[2\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}\right],$
(12)

where $\bm{x}_{1}$ and $\bm{x}_{2}$ are the positions of the two particles,
$r$ is their separation, $\epsilon$ is the interaction strength, and $\sigma$
is the particle radius.

###  C.3 Cylindrical buckling system

A 2D triangular mesh is used to model the cylindrical buckling system, with
the ends of the cylinder fixed in place to apply an axial compression of
0.14%. The energy of the system is evaluated by treating all bonds in the mesh
as an elastic spring to obtain the stretching energy, and all pairs of
adjacent triangles to be connected by elastic hinges, providing the bending
energy. Their expressions are given by

$E=\sum_{i}k^{\mathrm{S}}_{i}(r_{i}-r^{0}_{i})^{2}+\sum_{j}k^{\mathrm{B}}_{j}[1+\cos(\theta_{j}-\theta^{0}_{j})].$
(13)

The first term is the stretching energy, where $r_{i}$ is the length,
$r^{0}_{i}$ is the equilibrium length, and $k^{\mathrm{S}}_{i}$ is the
stretching rigidity of bond $i$. The second term provides the bending energy,
where $\theta_{j}$ is the dihedral angle, $\theta^{0}_{j}$ is the equilibrium
angle, and $k^{\mathrm{B}}_{j}$ is the bending rigidity of hinge $j$.

Figure 6:  Schematic of the bar and hinge model showing the relevant
parameters for a single hinge element. $h_{i}$ and $\bm{\hat{n}}_{i}$ denote
the height and unit normal of each triangle respectively.

The gradient of the energy can be obtained by individually considering the
stretching and bending energies of a single bond and hinge. For simplicity, we
will ignore the index for the bond and hinge. Using the variables shown in the
schematic in fig. 6, the gradient of the stretching energy of the bond between
$\bm{x}_{2}$ and $\bm{x}_{3}$ is given by,

$\frac{\partial E^{\mathrm{S}}}{\partial\bm{x}_{2}}=-\frac{\partial
E^{\mathrm{S}}}{\partial\bm{x}_{3}}=2k^{\mathrm{S}}(r-r^{0})(\bm{x_{2}}-\bm{x_{3}}).$
(14)

The gradients of the bending energy of the hinge are

$\displaystyle\frac{\partial E^{\mathrm{B}}}{\partial\bm{x}_{1}}$
$\displaystyle=k^{\mathrm{B}}\sin(\theta-\theta^{0})\frac{\bm{\hat{n}}_{a}}{h_{a}},$
(15) $\displaystyle\frac{\partial E^{\mathrm{B}}}{\partial\bm{x}_{2}}$
$\displaystyle=-k^{\mathrm{B}}\sin(\theta-\theta^{0})\left[\frac{\bm{\hat{n}}_{a}}{h_{a}}+\frac{\bm{\hat{n}}_{a}+\bm{\hat{n}}_{b}}{r}\right],$
(16) $\displaystyle\frac{\partial E^{\mathrm{B}}}{\partial\bm{x}_{3}}$
$\displaystyle=-k^{\mathrm{B}}\sin(\theta-\theta^{0})\left[\frac{\bm{\hat{n}}_{b}}{h_{b}}+\frac{\bm{\hat{n}}_{a}+\bm{\hat{n}}_{b}}{r}\right],$
(17) $\displaystyle\frac{\partial E^{\mathrm{B}}}{\partial\bm{x}_{4}}$
$\displaystyle=k^{\mathrm{B}}\sin(\theta-\theta^{0})\frac{\bm{\hat{n}}_{b}}{h_{b}}.$
(18)

###  C.4 Striped wetting system

This is modelled on a 200x200 2D grid (and 400x400 in section III.3) using a
phase-field model Panter, Gizaw, and Kusumaatmaja (2019), which has an order
parameter, $\phi(\bm{r})$, representing the phase of the liquid ($\phi=-1$ for
gas, $\phi=1$ for liquid). The energy functional has four separate terms,

$E[\phi]=E^{\mathrm{B}}[\phi]+E^{\mathrm{I}}[\phi]+E^{\mathrm{S}}[\phi]+E^{\mathrm{V}}[\phi].$
(19)

The first term, $E^{\mathrm{B}}$, uses a double well potential to set values
of $\phi=\pm 1$ in the bulk. The second term then provides the interfacial
energy between the liquid and gas by imposing an energy penalty to gradients
in $\phi$. $E^{\mathrm{S}}$ is the solid-liquid interaction energy, which sets
the contact angles of the hydrophilic and hydrophobic regions to
60$\mathrm{\SIUnitSymbolDegree}$ and 110$\mathrm{\SIUnitSymbolDegree}$,
respectively. Finally, $E^{\mathrm{V}}$ constrains the volume of the liquid
drop by penalising any variation from the target volume.

These four sections of the phase-field model are obtained using the following
equations,

$\displaystyle E^{\mathrm{B}}[\phi]$
$\displaystyle=\sum_{i}\frac{1}{\epsilon}\left(\frac{{\phi_{i}}^{4}}{4}-\frac{{\phi_{i}}^{2}}{2}+\frac{1}{4}\right)\mathop{}\\!\Delta
V,$ (20) $\displaystyle E^{\mathrm{I}}[\phi]$
$\displaystyle=\sum_{i}\frac{\epsilon}{2}\left|\bm{\nabla}\phi_{i}\right|^{2}\mathop{}\\!\Delta
V,$ (21) $\displaystyle E^{\mathrm{S}}[\phi]$
$\displaystyle=\sum_{j}\sqrt{2}\cos\theta_{j}\left(\frac{{\phi_{j}}^{3}}{6}-\frac{\phi_{j}}{2}-\frac{1}{3}\right)\mathop{}\\!\Delta
S,$ (22) $\displaystyle E^{\mathrm{V}}[\phi]$
$\displaystyle=k_{\mathrm{V}}\left[\sum_{i}\frac{\phi_{i}+1}{2}\mathop{}\\!\Delta
V-V_{0}\right]^{2},$ (23)

where the index $i$ includes all of the nodes, while $j$ represents the nodes
along the solid surface. $\mathop{}\\!\Delta V$ and $\mathop{}\\!\Delta S$
are, respectively, the volume and solid surface areas contained by each
individual node. $\epsilon$ is the liquid-gas interface width (set to 2.5
lattice units), $\theta_{j}$ is the contact-angle with the solid surface, and
$V_{0}$ is the constrained volume of the liquid drop. The strength of the
volume constraint is parametrised by $k_{\mathrm{V}}$ for which we use a value
of $10^{4}$. The gradients of these terms are given by,

$\displaystyle\frac{\partial E^{\mathrm{B}}}{\partial\phi_{i}}$
$\displaystyle=\frac{1}{\epsilon}\left({\phi_{i}}^{3}-\phi_{i}\right)\mathop{}\\!\Delta
V,$ (24) $\displaystyle\frac{\partial E^{\mathrm{I}}}{\partial\phi_{i}}$
$\displaystyle=\frac{\epsilon}{2}\left[\frac{\partial\left|\bm{\nabla}\phi_{i}\right|^{2}}{\partial\phi_{i}}+\sum_{k}\frac{\partial\left|\bm{\nabla}\phi_{k}\right|^{2}}{\partial\phi_{i}}\right]\mathop{}\\!\Delta
V,$ (25) $\displaystyle\frac{\partial E^{\mathrm{S}}}{\partial\phi_{j}}$
$\displaystyle=\sqrt{2}\cos\theta_{j}\left(\frac{{\phi_{j}}^{2}}{2}-\frac{1}{2}\right)\mathop{}\\!\Delta
S,$ (26) $\displaystyle\frac{\partial E^{\mathrm{V}}}{\partial\phi_{i}}$
$\displaystyle=k_{\mathrm{V}}\left[\sum_{i^{\prime}}\frac{\phi_{i^{\prime}}+1}{2}\mathop{}\\!\Delta
V-V_{0}\right]\mathop{}\\!\Delta V,$ (27)

where index $k$ denotes neighbouring nodes at which the evaluation of the
gradient uses $\phi_{i}$.

##  Appendix D Implementation details for the CINEB, DNEB-HEVF, and BITSS-HEVF
methods

The spring constants used to connect the states in the CINEB and DNEB methods
are system dependent. They are chosen such that they keep the states
equidistant without overwhelming the gradients arising from the potential
energy landscapes under consideration. We employ $10^{-1}$ for the Lennard-
Jones particle cluster, $10^{-2}$ for the cylindrical buckling example, and
$10^{-6}$ for the striped wetting system.

When hybrid eigenvector following is combined with BITSS or DNEB, we need to
set out criteria to determine when the double-ended method has sufficiently
converged, at which point the hybrid eigenvector following method begins. For
BITSS, the criterion is when the average of the two states changes by less
than a tenth of the reduction in the separation given by eq. 1 during a BITSS
step. For DNEB, the root-mean-square of the total energy gradient of the chain
of states is used with a convergence criterion of $10^{-3}$ for the particle
cluster and cylindrical buckling systems, and $10^{-5}$ for the striped
wetting. We employ the hybrid eigenvector following method implemented in the
program OPTIM Wales .

## References

  * Bryngelson  _et al._ (1995) J. D. Bryngelson, J. N. Onuchic, N. D. Socci, and P. G. Wolynes, Proteins: Structure, Function, and Bioinformatics 21, 167 (1995).
  * Onuchic, Luthey-Schulten, and Wolynes (1997) J. N. Onuchic, Z. Luthey-Schulten, and P. G. Wolynes, Annual Review of Physical Chemistry 48, 545 (1997).
  * Boehr  _et al._ (2006) D. D. Boehr, D. McElheny, H. J. Dyson, and P. E. Wrightt, Science 313, 1638 (2006).
  * Kerns  _et al._ (2015) S. J. Kerns, R. V. Agafonov, Y. J. Cho, F. Pontiggia, R. Otten, D. V. Pachov, S. Kutter, L. A. Phung, P. N. Murphy, V. Thai, T. Alber, M. F. Hagan, and D. Kern, Nature Structural and Molecular Biology 22, 124 (2015).
  * Guo  _et al._ (2018) C. Guo, Z. Wang, D. Wang, H. F. Wang, and P. Hu, Journal of Physical Chemistry C 122, 21478 (2018).
  * Richardson  _et al._ (2016) J. O. Richardson, C. Pérez, S. Lobsiger, A. A. Reid, B. Temelso, G. C. Shields, Z. Kisiel, D. J. Wales, B. H. Pate, and S. C. Althorpe, Science 351, 1310 (2016).
  * Vaillant, Althorpe, and Wales (2019) C. L. Vaillant, S. C. Althorpe, and D. J. Wales, Journal of Chemical Theory and Computation 15, 33 (2019).
  * Richard and Speck (2018) D. Richard and T. Speck, The Journal of Chemical Physics 148, 124110 (2018).
  * Wales, Miller, and Walsh (1998) D. J. Wales, M. A. Miller, and T. R. Walsh, Nature 394, 758 (1998).
  * Wales (2012) D. J. Wales, Philosophical Transactions of the Royal Society A 370, 2877 (2012).
  * Panter  _et al._ (2019) J. R. Panter, J. Chen, T. Zhang, and H. Kusumaatmaja, Communications Physics 2, 151 (2019).
  * Hutchinson and Thompson (2018) J. W. Hutchinson and J. M. T. Thompson, International Journal of Solids and Structures 148-149, 157 (2018).
  * Zhang and Ren (2014) Y. Zhang and W. Ren, The Journal of Chemical Physics 141, 244705 (2014).
  * Panter, Gizaw, and Kusumaatmaja (2019) J. R. Panter, Y. Gizaw, and H. Kusumaatmaja, Science Advances 5, eaav7328 (2019).
  * Amabili  _et al._ (2017) M. Amabili, A. Giacomello, S. Meloni, and C. M. Casciola, Physical Review Fluids 2, 034202 (2017).
  * Othayoth, Thoms, and Li (2020) R. Othayoth, G. Thoms, and C. Li, Proceedings of the National Academy of Sciences of the United States of America 117, 14987 (2020).
  * Filipov  _et al._ (2015) E. T. Filipov, T. Tachi, G. H. Paulino, and D. A. Weitz, Proceedings of the National Academy of Sciences of the United States of America 112, 12321 (2015).
  * Zhai, Wang, and Jiang (2018) Z. Zhai, Y. Wang, and H. Jiang, Proceedings of the National Academy of Sciences of the United States of America 115, 2032 (2018).
  * Bertoldi  _et al._ (2017) K. Bertoldi, V. Vitelli, J. Christensen, and M. Van Hecke, Nature Reviews Materials 2, 17066 (2017).
  * Truby and Lewis (2016) R. L. Truby and J. A. Lewis, Nature 540, 371 (2016).
  * Chi  _et al._ (2022) Y. Chi, Y. Li, Y. Zhao, Y. Hong, Y. Tang, and J. Yin, Advanced Materials 34, 2110384 (2022).
  * Bonfanti  _et al._ (2020) S. Bonfanti, R. Guerra, F. Font-Clos, D. Rayneau-Kirkhope, and S. Zapperi, Nature Communications 11, 4162 (2020).
  * Shan  _et al._ (2015) S. Shan, S. H. Kang, J. R. Raney, P. Wang, L. Fang, F. Candido, J. A. Lewis, and K. Bertoldi, Advanced Materials 27, 4296 (2015).
  * Giri and Mailen (2021) T. R. Giri and R. Mailen, International Journal of Mechanical Sciences 204, 106541 (2021).
  * Cerjan and Miller (1981) C. J. Cerjan and W. H. Miller, The Journal of Chemical Physics 75, 2800 (1981).
  * Heyden, Bell, and Keil (2005) A. Heyden, A. T. Bell, and F. J. Keil, Journal of Chemical Physics 123, 224101 (2005).
  * Kästner and Sherwood (2008) J. Kästner and P. Sherwood, Journal of Chemical Physics 128, 014106 (2008).
  * Zhang, Du, and Zheng (2016) L. Zhang, Q. Du, and Z. Zheng, SIAM Journal on Scientific Computing 38, A528 (2016).
  * E, Ren, and Vanden-Eijnden (2007) W. E, W. Ren, and E. Vanden-Eijnden, The Journal of Chemical Physics 126, 164103 (2007).
  * Ren and Vanden-Eijnden (2013) W. Ren and E. Vanden-Eijnden, Journal of Chemical Physics 138, 134105 (2013).
  * E, Ren, and Vanden-Eijnden (2002) W. E, W. Ren, and E. Vanden-Eijnden, Physical Review B 66, 052301 (2002).
  * Trygubenko and Wales (2004) S. A. Trygubenko and D. J. Wales, Journal of Chemical Physics 120, 2082 (2004).
  * Wales and Carr (2012) D. J. Wales and J. M. Carr, Journal of Chemical Theory and Computation 8, 5020 (2012).
  * Dewar, Healy, and Stewart (1984) M. J. S. Dewar, E. F. Healy, and J. J. P. Stewart, Journal of the Chemical Society, Faraday Transactions 2 80, 227 (1984).
  * Ionova and Carter (1993) I. V. Ionova and E. A. Carter, The Journal of Chemical Physics 98, 6377 (1993).
  * Miron and Fichthorn (2001) R. A. Miron and K. A. Fichthorn, The Journal of Chemical Physics 115, 8742 (2001).
  * Zhang, Shang, and Liu (2013) X.-J. Zhang, C. Shang, and Z.-P. Liu, Journal of Chemical Theory and Computation 9, 5745 (2013).
  * Trefethen and Thiyagalingam (2013) A. E. Trefethen and J. Thiyagalingam, Journal of Computational Science 4, 444 (2013).
  * Shalf (2020) J. Shalf, Philosophical Transactions of the Royal Society A 378, 20190061 (2020).
  * Alexander  _et al._ (2020) F. Alexander, A. Almgren, J. Bell, A. Bhattacharjee, J. Chen, P. Colella, D. Daniel, J. DeSlippe, L. Diachin, E. Draeger, A. Dubey, T. Dunning, T. Evans, I. Foster, M. Francois, T. Germann, M. Gordon, S. Habib, M. Halappanavar, S. Hamilton, W. Hart, Z. Huang, A. Hungerford, D. Kasen, P. R. Kent, T. Kolev, D. B. Kothe, A. Kronfeld, Y. Luo, P. Mackenzie, D. McCallen, B. Messer, S. Mniszewski, C. Oehmen, A. Perazzo, D. Perez, D. Richards, W. J. Rider, R. Rieben, K. Roche, A. Siegel, M. Sprague, C. Steefel, R. Stevens, M. Syamlal, M. Taylor, J. Turner, J. L. Vay, A. F. Voter, T. L. Windus, and K. Yelick, Philosophical Transactions of the Royal Society A 378, 20190056 (2020).
  * Lee and Bathe (1994) N.-S. Lee and K.-J. Bathe, Finite Elements in Analysis and Design 16, 99 (1994).
  * McMullen  _et al._ (2018) A. McMullen, M. Holmes-Cerfon, F. Sciortino, A. Y. Grosberg, and J. Brujic, Physical Review Letters 121, 138002 (2018).
  * Nguemaha and Zhou (2018) V. Nguemaha and H. X. Zhou, Scientific Reports 8, 1 (2018).
  * Chen, Zhang, and Torquato (2018) D. Chen, G. Zhang, and S. Torquato, Journal of Physical Chemistry B 122, 8462 (2018).
  * Santra, Singh, and Bagchi (2018) M. Santra, R. S. Singh, and B. Bagchi, Physical Review E 98, 032606 (2018).
  * Mökkönen, Ala-Nissila, and Jónsson (2016) H. Mökkönen, T. Ala-Nissila, and H. Jónsson, The Journal of Chemical Physics 145, 094901 (2016).
  * Wriggers (2006) P. Wriggers,  _Computational Contact Mechanics_ (Springer Berlin Heidelberg, Berlin, Heidelberg, 2006).
  * Page and McIver (1988) M. Page and J. W. McIver, The Journal of Chemical Physics 88, 922 (1988).
  * Liu and Nocedal (1989) D. C. Liu and J. Nocedal, Mathematical Programming 45, 503 (1989).
  * Henkelman, Uberuaga, and Jónsson (2000) G. Henkelman, B. P. Uberuaga, and H. Jónsson, The Journal of Chemical Physics 113, 9901 (2000).
  * Zimmerman (2013) P. Zimmerman, Journal of Chemical Theory and Computation 9, 3043 (2013).
  * Koslover and Wales (2007) E. F. Koslover and D. J. Wales, The Journal of Chemical Physics 127, 134102 (2007).
  * Sheppard, Terrell, and Henkelman (2008) D. Sheppard, R. Terrell, and G. Henkelman, The Journal of Chemical Physics 128, 134106 (2008).
  * Wales (2002) D. J. Wales, Molecular Physics 100, 3285 (2002).
  * Passerone and Parrinello (2001) D. Passerone and M. Parrinello, Physical Review Letters 87, 108302 (2001).
  * Virot  _et al._ (2017) E. Virot, T. Kreilos, T. M. Schneider, and S. M. Rubinstein, Physical Review Letters 119, 224101 (2017).
  * Kusumaatmaja  _et al._ (2006) H. Kusumaatmaja, J. Léopoldès, A. Dupuis, and J. M. Yeomans, Europhysics Letters 73, 740 (2006).
  * Brown and Bhushan (2016) P. S. Brown and B. Bhushan, Philosophical Transactions of the Royal Society A 374, 20160135 (2016).
  * Brakke (1992) K. A. Brakke, Experimental Mathematics 1, 141 (1992).
  * Kolev  _et al._ (2021) T. Kolev, P. Fischer, M. Min, J. Dongarra, J. Brown, V. Dobrev, T. Warburton, S. Tomov, M. S. Shephard, A. Abdelfattah, V. Barra, N. Beams, J.-S. Camier, N. Chalmers, Y. Dudouit, A. Karakus, I. Karlin, S. Kerkemeier, Y.-H. Lan, D. Medina, E. Merzari, A. Obabko, W. Pazner, T. Rathnayake, C. W. Smith, L. Spies, K. Swirydowicz, J. Thompson, A. Tomboulides, and V. Tomov, The International Journal of High Performance Computing Applications 35, 527 (2021).
  * Kirkpatrick, Gelatt, and Vecchi (1983) S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, Science 220, 671 (1983).
  * (62) D. J. Wales, “Optim: A program for optimizing geometries and calculating reaction pathways,” https://www-wales.ch.cam.ac.uk/software.html, accessed: 2022-08-31.


# Path integrals for fractional Brownian motion and fractional Gaussian noise

Baruch Meerson<EMAIL_ADDRESS>Racah Institute of Physics, Hebrew
University of Jerusalem, Jerusalem 91904, Israel  Olivier Bénichou
<EMAIL_ADDRESS>Sorbonne Université, CNRS, Laboratoire
de Physique Théorique de la Matière Condensée (UMR CNRS 7600), 4 Place
Jussieu, 75252 Paris Cedex 05, France  Gleb Oshanin  gleb.oshanin@sorbonne-
universite.fr Sorbonne Université, CNRS, Laboratoire de Physique Théorique de
la Matière Condensée (UMR CNRS 7600), 4 Place Jussieu, 75252 Paris Cedex 05,
France  Dipartimento di Scienze Matematiche, Politecnico di Torino, Corso Duca
degli Abruzzi 24, 10129 Torino, Italy

###### Abstract

The Wiener’s path integral plays a central role in the studies of Brownian
motion. Here we derive exact path-integral representations for the more
general _fractional_ Brownian motion (fBm) and for its time derivative process
– the fractional Gaussian noise (fGn). These paradigmatic non-Markovian
stochastic processes, introduced by Kolmogorov, Mandelbrot and van Ness, found
numerous applications across the disciplines, ranging from anomalous diffusion
in cellular environments to mathematical finance. Still, their exact path-
integral representations were previously unknown. Our formalism exploits the
Gaussianity of the fBm and fGn, relies on theory of singular integral
equations and overcomes some technical difficulties by representing the action
functional for the fBm in terms of the fGn for the sub-diffusive fBm, and in
terms of the derivative of the fGn for the super-diffusive fBm. We also extend
the formalism to include external forcing. The exact and explicit path-
integral representations open new inroads into the studies of the fBm and fGn.

Introduction. – The importance of path integrals in theoretical physics is
broadly recognized. Their application proved to be rewarding not only as a
computational tool, both analytical and numerical, but also as a powerful and
versatile conceptual framework. The notion of path integrals was introduced in
the 1920s by Wiener Wiener  for the Brownian motion (Bm). Since then it helped
uncover many nontrivial statistical properties of the Bm kac ; satya ; greg ;
vin ; dean . Feynman reinvented path integrals in the 1940s within his
reformulation of quantum mechanics Feynman ; hibbs ; Feynman2 . He is also
credited for making path integrals an intrinsic part of physicist’s toolbox
book0 ; wiegel ; book1 ; zinn ; book3 ; kamenevbook ; book4 ; leticia .

Path-integral representations of stochastic processes and fields are
especially useful in the studies of large-deviation statistics of physical
quantities. Performing a saddle-point evaluation of the pertinent path
integral (which relies on a problem-specific large parameter), one can
determine the optimal, that is the most likely, history of the system which
dominates the statistics in question. This method of large deviation analysis
appears in different areas of physics under different names: the optimal
fluctuation method, the instanton method, the weak-noise theory, the
macroscopic fluctuation theory, the dissipative WKB approximation, etc. A full
list of references on different applications of this method would exceed a
hundred.

The key object of a path-integral representation of Bm and its functionals is
the probability density $P[x(t)]$ of a given realization of a Brownian
trajectory $x(t)$, $P[x(t)]\sim\exp(-S[x(t)])$, where the action functional
$S[x(t)]$ is given by the Wiener’s formula Wiener

$S[x(t)]=\frac{1}{2}\int dt\,\dot{x}^{2}(t)$ (1)

(the dot here and henceforth denotes the time derivative, and we set the
diffusion coefficient to $1/2$ for brevity). The local-in-time Wiener’s action
(1) reflects the Markovian nature of the Bm. The last two decades have
witnessed a great interest in the _fractional_ Brownian motion (fBm),
introduced by Mandelbrot and van Ness Mandelbrot , and earlier by Kolmogorov
Kolmogorov . The Mandelbrot-van Ness (MvN) fBm is a non-Markovian
generalization of the Brownian motion which keeps the important properties of
Gaussianity, stationarity of the increment, and dynamical scale invariance.
For the two-sided (that is, pre-thermalized) fBm, time $t$ is defined on the
entire axis $|t|<\infty$. For the one-sided fBm $0\leq t<\infty$. Here the
process starts at $t=0$, and there is no past. Both versions of the fBm are
zero-mean Gaussian processes (for convenience we set $x(0)=0$), and they are
completely defined by their covariance functions

$\begin{split}\\!\kappa_{2}(t,t^{\prime})\\!=\\!\langle
x(t)x(t^{\prime})\rangle\\!&=\\!\frac{1}{2}\left(|t|^{2H}\\!+\\!|t^{\prime}|^{2H}\\!-\\!|t-t^{\prime}|^{2H}\right)\,,\\\
\kappa_{1}(t,t^{\prime})\\!=\\!\langle
x(t)x(t^{\prime})\rangle\\!&=\\!\frac{1}{2}\left(t^{2H}+t^{\prime
2H}-|t-t^{\prime}|^{2H}\right)\,.\end{split}$ (2)

Here the subscript 1 and 2 stand for the one- and two-sided processes,
respectively, the angle brackets denote ensemble averaging, and $0<H<1$ is the
Hurst index which quantifies the dynamical scale-invariance of the process
Stanley  and its ruggedness. For $H<1/2$ the fBm is sub-diffusive, i.e. the
mean-squared displacement $\\!\langle x^{2}(t)\rangle\\!=t^{2H}$ grows sub-
linearly with time. For $H>1/2$ the fBm is super-diffusive. In the borderline
case $H=1/2$ one recovers the standard Bm. Figure 1 presents examples of
numerical stochastic realizations of fBm for $H=1/4$, $1/2$ (standard Bm) and
$3/4$.

The fractional Gaussian noise (fGn) was introduced by Mandelbrot and van Ness
Mandelbrot  as the time-derivative of $x(t)$. That is, by definition, the fBm
$x(t)$ obeys the Langevin equation $\dot{x}(t)=y(t)$, where $y(t)$ is the fGn.
For $H<1/2$ the fGn is anti-persistent (that is, it has negative
autocorrelations). For $H>1/2$ it is positively correlated. For $H=1/2$ the
delta-correlated white Gaussian noise is recovered. The subsequent analysis
covers both sub- and super-diffusive cases.

Multiple physical processes have been successfully modeled as fBm. These
include fluctuating interfaces Krug  dynamics in crowded fluids weiss ; weiss2
, sub-diffusive dynamics of bacterial loci in a cytoplasm weber , telomere
diffusion in the cell nucleus garini ; prx , modeling of conformations of
serotonergic axons vojta , diffusion of a tagged bead of a polymer Walter ;
Amitai , translocation of a polymer through a pore Amitai ; Zoia ; Dubbeldam ;
Palyulin , single-file diffusion in ion channels Kukla ; Wei ; chanel , etc. A
review can be found in ralf . In its turn, the fGn Mandelbrot  is used to
model anti-persistent or persistent dependency structures in observed time
series in many applications including hydrology Molz , information theory
Barton , climate data analysis climate  and physiology physiol , to mention a
few.

By now the MvN fBm has become a standard model of anomalous diffusion in
systems with memory. Still, a satisfactory path-integral representation of
this process is unavailable without . This is in spite of the fact that, for
non-Markovian but Gaussian processes, such as the MvN fBm, there is a
straightforward path zinn  to constructing an analog of Eq. (1). It involves
the determination of a (highly singular) nonlocal kernel, inverse to the
covariance function (2), via solving a singular integral equation [such as Eq.
(5) below]. For MvN fBm this equation is hard to deal with analytically, which
explains the scarcity of results on path-integral representations of the fBm
RLfBm .

These technical obstacles were circumvented in the early work bead , where a
nonlocal analog of the Wiener’s action (1) was derived, by a different method,
in the particular case of the dynamics of a tagged bead in an infinitely long
pre-thermalized Rouse polymer bead . Under some natural assumptions this non-
Markovian system is equivalent to a fluctuating interface in one dimension,
and the latter is known to be describable by the MvN fBm with the Hurst
exponent $H=1/4$ Krug . The action, calculated in Ref. bead , is given, up to
a constant factor, by the expression

$S_{\text{bead}}[x(t)]\sim\iint\frac{dt_{1}dt_{2}}{|t_{1}-t_{2}|^{1/2}}\,\dot{x}(t_{1})\dot{x}(t_{2})\,,$
(3)

see also Ref. nech . We should also mention a series of works 19 ; 20 ; 21 ;
22 ; 23 ; khadem  aimed at determining $S[x(t)]$ for the one-sided MvN fBm in
the form of a perturbation expansion around the Wiener’s action (1). By
construction, such an expansion, based on the small parameter $|H-1/2|\ll 1$,
is quite limited in its validity.

In this work we find exact and explicit non-local analogs of the Wiener’s
action (1) for the MvN fBm: for arbitrary $0<H<1$ and for both two-sided and
one-sided versions of the fBm. We also extend the path integrals to include
overdamped motion of the particle under external force. We achieve these goals
by seeking, from the start, the action functional for the fBm in terms of its
time derivative processes: the first-derivative process (that is, the fGn) for
the sub-diffusive fBM, and the second-derivative process for the super-
diffusive fBm motivation . Our formalism fully exploits the Gaussianity of the
fBm and relies on the well-established theory of singular integral equations,
see e.g. Refs. kilbas ; polyanin . The resulting path integrals are convenient
to work with, as they involve only mildly singular kernels. Finally, we extend
the formalism to include external forcing.

Figure 1: Stochastic realizations of one-sided MvN fBm $x(t)$ for $H=1/4$ (top
left) and $H=3/4$ (bottom). $H=1/2$ (top right) corresponds to the standard
Bm.

General expressions and main results. – Quite generally, the action functional
$S[X(t)]$ of a Gaussian process $X(t)$ on a time interval $\Omega$ can be
represented as zinn

$S[X(t)]=\frac{1}{2}\int_{\Omega}dt_{1}\int_{\Omega}dt_{2}\,\mathcal{K}(t_{1},t_{2})\,X(t_{1})\,X(t_{2})\,.$
(4)

The kernel $\mathcal{K}(t_{1},t_{2})$ (a symmetric function of $t_{1}$ and
$t_{2}$) is the inverse of the covariance function $\kappa(t_{1},t_{2})$ of
the process $X(t)$:

$\int_{\Omega}dt_{1}\,\kappa\left(t_{1},t_{3}\right)\,\mathcal{K}(t_{1},t_{2})=\delta\left(t_{2}-t_{3}\right)\,.$
(5)

Once $\mathcal{K}(t_{1},t_{2})$ is known, the action functional (4) is
completely defined, giving the probability density $P[X(t)]$ of a given
realization of the process $X(t)$. Now we present our main results for the
action functionals of the MvN fBm $x(t)$. They have different forms for the
sub-diffusive and super-diffusive fBm, and for the two- and one-sided
processes.

We start with the sub-diffusion. For the two-sided sub-diffusive ($0<H<1/2$)
fBm $x(t)$, the action $S=S[x(t)]$ is given by

$S=\frac{{\rm cot}(\pi H)}{4\pi
H}\int^{\infty}_{-\infty}\int^{\infty}_{-\infty}\frac{dt_{1}dt_{2}}{|t_{1}-t_{2}|^{2H}}\,\dot{x}(t_{1})\dot{x}(t_{2})\,.$
(6)

For the one-sided sub-diffusive fBm we obtain

$S\\!=\\!\frac{{\rm cot}(\pi H)}{4\pi
H}\\!\int^{\infty}_{0}\int^{\infty}_{0}\\!\\!\frac{dt_{1}dt_{2}\,I_{z}\left(\frac{1}{2}-H,H\right)}{|t_{1}-t_{2}|^{2H}}\dot{x}(t_{1})\dot{x}(t_{2}),$
(7)

where $I_{z}(a,b)$ is the regularized incomplete beta function

$I_{z}(a,b)=\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\int^{z}_{0}x^{a-1}(1-x)^{b-1}dx\,,$
(8)

$\Gamma(\dots)$ is the gamma function, and $z=4t_{1}t_{2}(t_{1}+t_{2})^{-2}$.

As one can see, the action functionals (6) and (7) are non-local in time and
written in terms of the fGn $\dot{x}(t)$, rather than in terms of $x(t)$
itself. The expression (7) for the one-sided case is more complicated, than
that for the two-sided one, Eq. (6). In particular, the two-sided kernel in
Eq. (6) is a difference kernel, which reflects the stationarity in time of the
two-sided derivative process, the fGn. The one-sided kernel (7) is not a
difference kernel in spite of the stationarity of the fGn. The non-
stationarity, however, is temporary, as it is caused by a transient created by
the initial condition $x(t=0)=0$. Indeed, in the limit of
$t_{1},t_{2}\to\infty$, and $t_{1}-t_{2}=\text{const}$, $z$ tends to $1$, the
one-sided kernel coincides with the two-sided one, and the stationarity is
restored.

In the limiting case $H=1/2$, the kernels in Eqs. (6) and (7) become delta-
functions and yield the classical Wiener’s formula (1), as we show in Ref. SM
.

For $H=1/4$ Eq. (6) has the same functional form as the two-sided expression
(3), as to be expected in view of the pre-thermalization of the Rouse polymer
bead . We also remark that Eq. (6) was postulated in Ref. nech  as an
effective Hamiltonian of topologically stabilized polymers in melts,
permitting to cover various conformations ranging from ideal Gaussian coils to
crumpled globules. Our derivation validates their approach.

Now we present our results for the super-diffusive fBm, $1/2<H<1$. In the two-
sided case we obtain

$S\\!=\\!\frac{\sigma(H)}{2}\,\\!\int^{\infty}_{-\infty}\\!\int^{\infty}_{-\infty}dt_{1}\,dt_{2}\,|t_{1}-t_{2}|^{2-2H}\,\ddot{x}(t_{1})\,\ddot{x}(t_{2})\,,$
(9)

where

$\sigma(H)=-\frac{{\rm cot}(\pi H)}{4\pi H(1-H)(2H-1)}\,,$ (10)

a positive function. For the one-sided case

$\displaystyle S$ $\displaystyle=$
$\displaystyle\frac{\sigma(H)}{2}\,\int^{\infty}_{0}\int^{\infty}_{0}dt_{1}\,dt_{2}\,|t_{1}-t_{2}|^{2-2H}$
(11) $\displaystyle\times$ $\displaystyle
I_{z^{\prime}}\left(\frac{3}{2}-H,2H-2\right)\ddot{x}(t_{1})\ddot{x}(t_{2})\,,$

where $z^{\prime}={\rm min}(t_{1},t_{2})/{\rm max}(t_{1},t_{2})$. Again, the
expressions in Eqs. (9) and (11) are non-local in time, but now they are
written in terms of $\ddot{x}(t)$, that is in terms of the first derivative of
the fGn. The two-sided kernel is a difference kernel. The one-sided kernel is
not, but it approaches the difference form following an initial transient.
Also, the classical Wiener’s form (1) is recovered in the limit $H\to 1/2$ SM
.

Expressions (6)-(11), alongside with Eqs. (25) and (27) below, represent the
main results of this work. Here we present derivations of Eqs. (6) and (9) for
the two-sided sub-diffusive and super-diffusive fBm, respectively. The
derivation of the (a bit more bulky) one-sided expressions in Eqs. (7) and
(11) is relegated to the SM SM .

Sub-diffusion – Here we work directly with the fGn. Its covariance function
$c(t_{1},t_{2})$ can be readily calculated:

$\displaystyle c(t_{1},t_{2})=\langle
y(t_{1})y(t_{2})\rangle=\langle\dot{x}(t_{1})\dot{x}(t_{2})\rangle$
$\displaystyle=\frac{\partial^{2}}{\partial t_{1}\partial t_{2}}\langle
x(t_{1})x(t_{2})\rangle=\frac{d}{d\tau}\left(H|\tau|^{2H-1}\text{sgn}\,\tau\right)\,,$
(12)

where $\tau=t_{1}-t_{2}$, and we used Eq. (2). Equation (12) holds both for
the two-sided and the one-sided process and for all $0<H<1$. Notably, the fGn
is a stationary process. For $H=1/2$ Eq. (12) gives
$c(\tau)=(1/2)(d/d\tau)\,\text{sgn}\,\tau=\delta(\tau)$, as to be expected for
the white noise.

Let us denote by $\mathcal{C}(\tau)$ the kernel inverse to $c(\tau)$. For the
two-sided process, $\mathcal{C}(\tau)$ is defined by the equation
$\int_{-\infty}^{\infty}d\tau\,c\left(\tau-t\right)\,\mathcal{C}(\tau)=\delta(t)$
or, in the explicit form,

$\int_{-\infty}^{\infty}d\tau\,\mathcal{C}(\tau)\frac{d}{d\tau}\left[|\tau-t|^{2H-1}\text{sgn}(\tau-t)\right]=\frac{1}{H}\,\delta(t).$
(13)

Integrating by part and assuming that the boundary terms are zero (as can be
verified a posteriori), we arrive at the integral equation

$\int_{-\infty}^{\infty}d\tau\,\frac{\text{sgn}(\tau-t)}{|\tau-t|^{1-2H}}\,\mathcal{D}(\tau)=-\frac{1}{H}\,\delta(t)$
(14)

for the unknown function $\mathcal{D}(\tau)=d\mathcal{C}(\tau)/d\tau$. The
solution can be found in Ref. kilbas :

$\mathcal{D}(\tau)=\frac{d\mathcal{C}(\tau)}{d\tau}=\frac{\text{cot}\,(\pi
H)}{2\pi H}\frac{d}{d\tau}\frac{1}{|\tau|^{2H}}\,.$ (15)

Getting rid of the $\tau$-derivative and using the fact that the kernel must
vanish at $|\tau|\to\infty$, we obtain

$\mathcal{C}(\tau)=\frac{\text{cot}\,(\pi H)}{2\pi H}\frac{1}{|\tau|^{2H}}\,.$
(16)

The ensuing Gaussian action functional (4), written in terms of
$X(t)=\dot{x}(t)$, yields the announced equation (6).

Super-diffusion – Here we work with the second-derivative process
$z(t)=\ddot{x}(t)$. Its covariance is

$q(t_{1},t_{2})=\frac{d^{3}}{d\tau^{3}}\left[H|\tau|^{2H-1}\text{sgn}(\tau)\right]\,.$
(17)

For the two-sided process the inverse kernel $\mathcal{Q}(t_{1},t_{2})$ is
defied by the equation
$\int_{-\infty}^{\infty}d\tau\,q\left(\tau-t\right)\,\mathcal{Q}(\tau)=\delta(t)$,
or, in the explicit form,

$\int_{-\infty}^{\infty}d\tau\,\mathcal{Q}(\tau)\frac{d^{3}}{d\tau^{3}}\left[|\tau-t|^{2H-1}\text{sgn}(\tau-t)\right]=\frac{1}{H}\,\delta(t).$
(18)

Integrating three times by part and assuming that the boundary terms are zero
(as verified a posteriori), we arrive at the equation

$\int_{-\infty}^{\infty}d\tau\,\frac{\text{sgn}(\tau-t)}{|\tau-t|^{1-2H}}\,\mathcal{Z}(\tau)=-\frac{1}{H}\,\delta(t),$
(19)

where $\mathcal{Z}(\tau)=d^{3}\mathcal{Q}(\tau)/d\tau^{3}$. This is exactly
the same equation as Eq. (14), but now $1-2H<0$. It is convenient to rewrite
this equation as

$\displaystyle-\int_{-\infty}^{t}d\tau\,\left(t-\tau\right)^{2H-1}\mathcal{Z}(\tau)$
$\displaystyle+$
$\displaystyle\int_{t}^{\infty}d\tau\,\left(\tau-t\right)^{2H-1}\mathcal{Z}(\tau)$
(20) $\displaystyle=$ $\displaystyle-\frac{1}{H}\,\delta(t)$

and differentiate both sides of Eq. (20) with respect to $t$. The resulting
equation,

$\int_{-\infty}^{\infty}d\tau\,\frac{\mathcal{Z}(\tau)}{|\tau-t|^{2-2H}}=\frac{1}{H(2H-1)}\,\delta^{\prime}(t),$
(21)

is solvable kilbas , and we obtain

$\mathcal{Z}(\tau)=-\frac{\text{cot}(\pi H)}{2\pi
H(2H-1)}\,\frac{d^{2}}{d\tau^{2}}\frac{\text{sgn}\,\tau}{|\tau|^{2H-1}}\,.$
(22)

Integrating this expression over $\tau$ three times and using account the fact
that the kernel must vanish at $|\tau|\to\infty$, we obtain the desired
inverse kernel:

$\mathcal{Q}(\tau)=\sigma(H)\,|\tau|^{2(1-H)}\,,$ (23)

where $\sigma(H)$ is defined in Eq. (10). The resulting Gaussian action
functional (4), written in terms of $X(t)=z(t)\equiv\ddot{x}(t)$, yields the
announced Eq. (9).

External force – An important extension of this formalism deals with
situations where the fBm of a particle is accompanied by its overdamped motion
under external force $f(x)$. A natural approach to modelling this situation
employs the non-Markovian Langevin equation Metzler2021

$\dot{x}(t)=f[x(t)]+y(t)\,,$ (24)

where the noise term $y(t)$ describes fGn. When the external force $f(x)$ is
confining, the $x$-distribution approaches a steady state. This steady state,
however, is non-Boltzmann. Therefore, not surprisingly, it violates the
fluctuation-dissipation theorem Metzler2021 . As the fGn $y(t)$ is a Gaussian
process, a natural path-integral representation for Eq. (24) is provided by
the action functional

$\displaystyle S[x(t)]$ $\displaystyle=$
$\displaystyle\frac{1}{2}\int_{-\infty}^{\infty}dt_{1}\int_{-\infty}^{\infty}dt_{2}\,\mathcal{C}(t_{1}-t_{2})\,\\{\dot{x}(t_{1})-f[x(t_{1})]\\}$
(25) $\displaystyle\times$ $\displaystyle\\{\dot{x}(t_{2})-f[x(t_{2})]\\}\,,$

where $\mathcal{C}(\tau)$ is the inverse kernel for the fGn, given by Eq.
(16). Here we assumed a two-sided sub-diffusive fBm.

For a super-diffusive fBm a suitable non-Markovian Langevin equation can be
obtained by a formal differentiation of Eq. (24) with respect to time, leading
to

$\ddot{x}(t)=f^{\prime}[x(t)]\dot{x}(t)+z(t)\,,$ (26)

where $f^{\prime}(x)\equiv df(x)/dx$, and the noise term $z(t)$ is the time
derivative of the fGn. The corresponding path integral for the two-sided
process is given by the action functional

$\displaystyle
S[x(t)]=\frac{1}{2}\int_{-\infty}^{\infty}dt_{1}\int_{-\infty}^{\infty}dt_{2}\,\mathcal{Q}(t_{1}-t_{2})$
(27) $\displaystyle\times$
$\displaystyle\\!\\!\left[\ddot{x}(t_{1})-f^{\prime}[x(t_{1})]\,\dot{x}(t_{1})\right]\left[\ddot{x}(t_{2})-f^{\prime}[x(t_{2})]\,\dot{x}(t_{2})\right],$

where $\mathcal{Q}(\tau)$ is given by Eq. (23). Expressions similar to Eqs.
(25) and (27), but with the one-side kernels as in Eqs. (7) and (11), hold for
the one-sided sub- and super-diffusive fBm, respectively.

Summary. – We generalized the classical Wiener’s path integral for the Bm and
found exact path-integral representations for the two-sided and one-sided MvN
fBm for the whole range $0<H<1$ of the Hurst exponent. We also extended the
formalism to include external forcing. The exact and explicit path-integral
representations open new inroads into analytical and numerical studies of fBm
– an important paradigm of scale-invariant stochastic processes with memory –
in a multitude of applications in natural sciences, technology and finance.

Acknowledgments. – We are grateful to P. Chigansky, D. S. Dean, S. N. Majumdar
and K. L. Sebastian for useful discussions. B. M. was supported by the Israel
Science Foundation (Grant No. 1499/20).

##  Appendix A Supplemental Material

Here we present some details of derivations of the results obtained for the
one-sided case and also show that in the limit $H\to 1/2$ the expressions (6),
(7), (9) and (11) converge to the Wiener result in Eq. (1).

##  Appendix B Two-sided fBm

###  B.1 Sub-diffusion, $0<H<1/2$. The limit of $H\to 1/2$

In order to take the limit $H\to 1/2$, we take advantage of the identity

$\displaystyle\frac{1}{|t_{1}-t_{2}|^{2H}}=\frac{1}{(1-2H)}\frac{d}{dt_{1}}\frac{{\rm
sgn}(t_{1}-t_{2})}{|t_{1}-t_{2}|^{2H-1}}\,,$ (S1)

which permits us to formally rewrite the kernel in Eq. (3) as

$\displaystyle\frac{{\rm cot}(\pi H)}{2\pi
H}\frac{1}{|t_{1}-t_{2}|^{2H}}=\frac{{\rm cot}(\pi H)}{2\pi
H(1-2H)}\frac{d}{dt_{1}}\frac{{\rm
sgn}(t_{1}-t_{2})}{|t_{1}-t_{2}|^{2H-1}}\,.$ (S2)

Taking the limit $H\to 1/2$ in the both sides of the latter equality and
noticing that

$\displaystyle\lim_{H\to 1/2}\frac{{\rm cot}(\pi H)}{(1-2H)}=\frac{\pi}{2}\,,$
(S3)

we get

$\displaystyle\lim_{H\to 1/2}\frac{{\rm cot}(\pi H)}{2\pi
H}\frac{1}{|t_{1}-t_{2}|^{2H}}$ $\displaystyle=\lim_{H\to 1/2}\frac{{\rm
cot}(\pi H)}{2\pi H(1-2H)}\frac{d}{dt_{1}}\frac{{\rm
sgn}(t_{1}-t_{2})}{|t_{1}-t_{2}|^{2H-1}}$ (S4)
$\displaystyle=\frac{1}{2}\frac{d}{dt_{1}}{\rm
sgn}(t_{1}-t_{2})=\delta(t_{1}-t_{2})\,,$ (S5)

which yields the Wiener expression (1).

###  B.2 Super-diffusion, $1/2<H<1$. The limit of $H\to 1/2$

We turn to the limit $H\to 1/2$ directly in Eq. (22) in the main text to get

$\displaystyle\lim_{H\to
1/2}\mathcal{Z}(\tau)=\mathcal{Z}_{H=1/2}(\tau)=\frac{1}{2}\frac{d^{2}}{d\tau^{2}}{\rm
sgn}(\tau)=\frac{d}{d\tau}\delta(\tau)$ (S6)

The action $S$ written in terms of the derivative of the fractional Gaussian
noise involves the kernel function $\mathcal{Q}$, which is given by a triple
integral of $\mathcal{Z}_{H=1/2}(\tau)$ [see Eq. (23)]. Consequently, the
action has the form

$\displaystyle
S=\frac{1}{2}\int^{\infty}_{-\infty}\int^{\infty}_{-\infty}dt_{1}dt_{2}\ddot{x}(t_{1})\ddot{x}(t_{2})\int^{t_{1}}_{-\infty}\int^{t_{2}}_{-\infty}d\tau_{1}d\tau_{2}\delta(\tau_{1}-\tau_{2})\,.$
(S7)

Integrating the latter expression by parts, we arrive at the Wiener’s result
in Eq. (1).

##  Appendix C One-sided fBm

Here we present brief derivations of our expressions (7) and (11) of the main
text.

###  C.1 Sub-diffusion, $0<H<1/2$

For the one-sided sub-diffusive fBm the inverse kernel
$\mathcal{C}(t_{1},t_{2})$ is not a difference kernel, and Eq. (13) gives way
to the equation

$\displaystyle\int_{0}^{\infty}dt_{1}\,\mathcal{C}(t_{1},t_{2})\frac{d}{dt_{1}}\left[|t_{1}-t_{3}|^{2H-1}\text{sgn}(t_{1}-t_{3})\right]=\frac{1}{H}\,\delta(t_{2}-t_{3})\,.$
(S8)

Integrating by part, we arrive at

$\int_{0}^{\infty}dt_{1}\,\frac{\text{sgn}(t_{1}-t_{3})}{|t_{1}-t_{3}|^{1-2H}}\,\mathcal{D}(t_{1},t_{2})=-\frac{1}{H}\,\delta(t_{2}-t_{3}),$
(S9)

where $\mathcal{D}(t_{1},t_{2})=(\partial/\partial
t_{1})\,\mathcal{C}(t_{1},t_{2})$. The solution can be found in Ref. [47].
Getting rid of the $t_{1}$-derivative, we obtain $\mathcal{C}(t_{1},t_{2})$ in
the following three alternative (but equivalent) forms

$\displaystyle\mathcal{C}(t_{1},t_{2})$
$\displaystyle=\frac{(t_{1}t_{2})^{-H}\sqrt{z^{\prime}}}{H(1-2H)\sin(\pi
H)\Gamma(2H)\Gamma^{2}\left(\frac{1}{2}-H\right)}\,_{2}F_{1}\left(1,H+\frac{1}{2};\frac{3}{2}-H;z^{\prime}\right)\,,\quad
z^{\prime}=\frac{{\rm min}(t_{1},t_{2})}{{\rm max}(t_{1},t_{2})}$ (S10)
$\displaystyle=\frac{{\rm cot}(\pi H)(t_{1}t_{2})^{-H}\sqrt{z}}{2\pi
H(1-2H)B(1/2-H,2H)}\,_{2}F_{1}\left(\frac{1}{2},1;\frac{3}{2}-H;z\right)\,,\quad
z=\frac{4t_{1}t_{2}}{(t_{1}+t_{2})^{2}}$ (S11) $\displaystyle=\frac{{\rm
cot}(\pi
H)}{2\sqrt{\pi}\Gamma(1/2-H)\Gamma(1+H)}\frac{B_{z}(1/2-H,H)}{|t_{1}-t_{2}|^{2H}}\,,$
(S12)

where $B(a,b)$ and $B_{z}(a,b)$ are the complete and incomplete beta-
functions, respectively, and ${}_{2}F_{1}(\dots)$ is the hypergeometric
function. Recalling next the definition of the regularized incomplete beta-
function [see Eq. (8)], we obtain our result in Eq. (7).

Limit $H\to 1/2$. To take the limit $H\to 1/2$ in Eq. (7) in the main text it
is expedient to use the representations of the kernel
$\mathcal{C}(t_{1},t_{2})$ given in the second line in Eq. (S10). We formally
rewrite the Gauss hypergeometric function entering this representation as

$\,{}_{2}F_{1}\left(\frac{1}{2},1;\frac{3}{2}-H;z\right)=\frac{\Gamma(H)\Gamma(3/2-H)}{\sqrt{\pi}(1-z)^{H}z^{1/2-H}}-\frac{(1/2-H)}{H}\,_{2}F_{1}\left(\frac{1}{2},1;1+H;1-z\right)$
(S13)

such that, after some algebra, the kernel $\mathcal{C}(t_{1},t_{2})$ can be
cast into the form

$\displaystyle\mathcal{C}(t_{1},t_{2})=\frac{{\rm cot}(\pi H)}{2\pi
H}\frac{1}{|t_{1}-t_{2}|^{2H}}-\frac{{\rm cot}(\pi H)}{2\pi
H^{2}B(1/2-H,2H)}\frac{(t_{1}t_{2})^{1/2-H}}{(t_{1}+t_{2})}\,_{2}F_{1}\left(\frac{1}{2},1;1+H;\left(\frac{t_{1}-t_{2}}{t_{1}+t_{2}}\right)^{2}\right)$
(S14)

Further on, we observe that the numerical $H$-dependent amplitude in the
second term in the latter expression vanishes in the limit $H\to 1/2$ as

$\displaystyle\frac{{\rm cot}(\pi H)}{2\pi H^{2}B(1/2-H,2H)}\simeq
2\left(\frac{1}{2}-H\right)^{2}\,,$ (S15)

which signifies that this term does not contribute in this limit. On the
contrary, as demonstrated above, the first term converges to the delta-
function which ensures that the action in Eq. (7) in the main text converges
to the Wiener’s result in Eq. (1).

###  C.2 Super-diffusion, $1/2<H<1$

Here we present a derivation of Eq. (11) and also check the limit of $H\to
1/2$. In particular, this derivation highlights the reason why the
representation of the action in terms of $\ddot{x}$ is advantageous.

We start the derivation by representing the action in terms of the first
derivative of the fBm, that is in terms of the fGn:

$\displaystyle
S=\frac{1}{2}\int^{\infty}_{0}\int^{\infty}_{0}dt_{1}dt_{2}\mathcal{C}(t_{1},t_{2})\dot{x}(t_{1})\dot{x}(t_{2})\,.$
(S16)

The inverse kernel $\mathcal{C}(t_{1},t_{2})$ obeys the integral equation

$\displaystyle\int^{\infty}_{0}\frac{dt_{1}\,\mathcal{C}(t_{1},t_{2})}{|t_{1}-t_{3}|^{2-2H}}=\frac{\delta(t_{2}-t_{3})}{H(2H-1)}\,.$
(S17)

The explicit solution of this equation can be found in Ref. Lundgren . After
straightforward transformations, it reads:

$\displaystyle\mathcal{C}(t_{1},t_{2})$ $\displaystyle=\frac{1}{2H\sin(\pi
H)\Gamma(2H)\Gamma^{2}(3/2-H)}\frac{d^{2}}{dt_{1}dt_{2}}\int_{0}^{{\rm
min}(t_{1},t_{2})}\frac{d\tau}{(t_{1}-\tau)^{H-1/2}(t_{2}-\tau)^{H-1/2}}$
(S18) $\displaystyle=\frac{1}{2H\sin(\pi
H)\Gamma(2H)\Gamma^{2}(3/2-H)}\frac{d^{2}}{dt_{1}dt_{2}}|t_{1}-t_{2}|^{2-2H}B_{z^{\prime}}(3/2-H,2H-2)\,,\quad
z^{\prime}=\frac{{\rm min}(t_{1},t_{2})}{{\rm max}(t_{1},t_{2})}\,.$ (S19)

We observe that the function $\mathcal{C}(t_{1},t_{2})$, defined in the right-
hand-side of Eq. (S18), contains a non-integrable singularity at
$t_{1}=t_{2}$, which shows why the representation of the action in Eq. (S16)
in terms of the first derivatives $\dot{x}(t)$ is problematic. To get a
regular result, we integrate Eq. (S16) by part, i.e. express it in terms of
the second derivatives $\ddot{x}(t_{1})$ and $\ddot{x}(t_{2})$. Then, using
the definition of the regularized incomplete beta-function in Eq. (8), we
arrive at the final result in Eq. (11).

Limit $H\to 1/2$. The limit $H\to 1/2$ can be conveniently taken in the
expression given in the first line of Eq. (S18). This yields

$\displaystyle\lim_{H\to
1/2}\mathcal{C}(t_{1},t_{2})=\frac{d^{2}}{dt_{1}dt_{2}}{\rm
min}(t_{1},t_{2})=\delta(t_{1}-t_{2})\,,$ (S20)

and we recover the Wiener’s expression in Eq. (1).

## References

  * (1) N. Wiener, Proc. Natl. Acad. Sci. USA 7, 253 (1921); Proc. Natl. Acad. Sci. USA 7, 294 (1921); J. Math. Phys. Sci. 2, 132 (1923); Proc. London Math. Soc. 22, 454 (1924); Acta Math. 55, 117 (1930).
  * (2) M. Kac, Trans. Amer. Math. Soc. 65, 1 (1949).
  * (3) S. N. Majumdar, Curr. Sci. 89, 2076 (2005).
  * (4) A. J. Bray, S. N. Majumdar, and G. Schehr, Adv. Phys. 62, 225 (2013).
  * (5) V. Démery and D. S. Dean, Phys. Rev. E 84, 011148 (2011).
  * (6) D. S. Dean and R. Horgan, Phys. Rev. E 76, 041102 (2007).
  * (7) R. P. Feynman, Rev. Mod. Phys. 20, 367 (1948).
  * (8) R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals (McGraw-Hill, New York, 1965).
  * (9) Feynman’s Thesis – A New Approach to Quantum Theory, edited by L. M. Brown (World Scientific, Singapor, 2005).
  * (10) L. S. Schulman, Techniques and Applications of Path Integration (Wiley, New York, 1981).
  * (11) F.W. Wiegel, Introduction to Path-Integral Methods in Physics and Polymer Science, (World Scientific, Philadelphia, 1986).
  * (12) M. Chaichian and A. P. Demichev, Stochastic Processes and Quantum Mechanics (Path Integrals in Physics., vol. 1.) (Institute of Physics Publishing, Bristol, 2001).
  * (13) J. Zinn-Justin. Quantum Field Theory and Critical Phenomena (International Series of Monographs on Physics, vol. 113) (Clarendon, Oxford, 2002).
  * (14) H. Kleinert, Path Integrals in Quantum Mechanics, Statistics, Polymer Physics, and Financial Markets (World Scientific, Singapore, 2009).
  * (15) A. Kamenev, Field Theory of Non-Equilibrium Systems (Cambridge University Press, Cambridge, UK, 2011).
  * (16) H. S. Wio, Path Integrals for Stochastic Processes. An Introduction (World Scientific, Singapore, 2013).
  * (17) L. F. Cugliandolo, V. Lecomte, and F. van Wijland, J. Phys. A: Math. Theor. 52 50LT01 (2019).
  * (18) B. B. Mandelbrot and J. W. van Ness, SIAM Review 10, 422 (1968).
  * (19) A. N. Kolmogorov, CR (Dokl.) Acad. Sci. URSS 26, 115 (1940).
  * (20) A.-L. Barabási and H. E. Stanley, Fractal Concepts in Surface Growth (Cambridge University Press, Cambridge, England, 1995).
  * (21) J. Krug, H. Kallabis, S. N. Majumdar, S. J. Cornell, A. J. Bray,and C. Sire, Phys. Rev. E 56, 2702 (1997).
  * (22) M. Weiss, Phys. Rev. E 88, 010101(R) (2013).
  * (23) D. Ernst, M. Hellmann, J. Köhler, and M. Weiss, Soft Matter 8, 4886 (2012).
  * (24) S. C. Weber, A. J. Spakowitz, and J. A. Theriot, Phys. Rev. Lett. 104, 238102 (2010).
  * (25) I. Bronshtein, E. Kepten, I. Kanter, S. Berezin, M. Lindner, A. B. Redwood, S. Mai, S. Gonzalo, R. Foisner, Y. Shav-Tal, and Y. Garini, Nat. Commun. 6, 8044 (2015).
  * (26) D. Krapf, N. Lukat, E. Marinari, R. Metzler, G. Oshanin, C. Selhuber-Unkel, A. Squarcini, L. Stadler, M. Weiss, and X. Xu, Phys. Rev. X 9, 011019 (2019).
  * (27) S. Janusonis, N. Detering, R. Metzler, and T. Vojta, Frontiers Comp. Neurosci. 14, 56 (2020).
  * (28) J.-C. Walter, A. Ferrantini, E. Carlon, and C. Vanderzande, Phys. Rev. E 85, 031120 (2012).
  * (29) A. Amitai, Y. Kantor, and M. Kardar, Phys. Rev. E 81, 011107 (2010).
  * (30) A. Zoia, A. Rosso, and S. N. Majumdar, Phys. Rev. Lett. 102, 120602 (2009).
  * (31) J. L. A. Dubbeldam, V. G. Rostiashvili, A. Milchev, and T. A. Vilgis, Phys. Rev. E 83, 011802 (2011).
  * (32) V. Palyulin, T. Ala-Nissila, and R. Metzler, Soft Matter 10, 9016 (2014).
  * (33) V. Kukla, J. Kornatowski, D. Demuth, I. Girnus, H. Pfeifer, L. V. C. Rees, S. Schunk, K. K. Unger, and J. K¨arger, Science 272, 702 (1996).
  * (34) Q.-H. Wei, C. Bechinger, and P. Leiderer, Science 287, 625 (2000).
  * (35) O. Bénichou, P. Illian, G. Oshanin, A. Sarracino, and R. Voituriez, J. Phys.: Condens. Matter 30, 443001 (2018).
  * (36) R. Metzler, J.-H. Jeon, A. G. Cherstvy, and E. Barkai, Phys. Chem. Chem. Phys. 16, 24128 (2014).
  * (37) F. J. Molz, H. H. Liu, and J. Szulga, Water Resources Res. 33, 2273 (1997).
  * (38) R.J. Barton and H.V. Poor, IEEE Trans. Inform. Theory 34, 943 (1988).
  * (39) E. Myrvoll-Nilsen, H.-B. Fredriksen, S. H. Sørbye, and M. Rypdal, Front. Earth Sci. 7, 214 (2019).
  * (40) V. Maxim, L. Sendur, J. Fadili, J. Suckling, R. Gould, R. Howard, and E. Bullmore, Neuroimage 25, 141 (2005).
  * (41) In simple large-deviation problems, analyzed with the optimal fluctuation method, the explicit knowledge of the inverse kernel can be unnecessary. This happens when the ensuing nonlocal Euler-Lagrange equation (a linear integral equation) can be transformed into a form containing only the covariance function of the process baruch ; longtimeLD ; escape . For more complicated problems, however, no such transformation exists, and the knowledge of the inverse kernel is indispensable.
  * (42) For the Riemann-Liouville fBm (a different mathematical model of anomalous diffusion, introduced by Levy Levy ), the action $S[x(t)]$ can be determined by employing fractional calculus, as was done in Refs. seb ; san  and in Chapter 10 in Ref. book4 . Also, in Refs. Calvo ; Wiopaper  a path-integral representation was found for the fractional Levy motion Calvo , still another mathematical model. All these results, however, are irrelevant to the MvN fBm studied in our work.
  * (43) S. F. Burlatskii and G. Oshanin, Theor. Math. Phys. 75(3), 659 (1988).
  * (44) K. Polovnikov, S. Nechaev, and M. V. Tamm, Soft Matter 14, 6561 (2018).
  * (45) K. J. Wiese, S. N. Majumdar, and A. Rosso, Phys. Rev. E 83, 061141 (2011).
  * (46) M. Delorme and K. J. Wiese, Phys. Rev. Lett. 115, 210601 (2015).
  * (47) M. Delorme and K. J. Wiese, Phys. Rev. E 94, 012134 (2016).
  * (48) M. Delorme and K. J. Wiese, Phys. Rev. E 94, 052105 (2016).
  * (49) T. Sadhu, M. Delorme and K. J. Wiese, Phys. Rev. Lett. 120, 040603 (2018).
  * (50) S. M. J. Khadem, R. Klages and S. H. L. Klapp, arXiv:2205.15791.
  * (51) Indeed, although $x(t)$ is a Gaussian process, the Wiener’s action (1) is expressed through $\dot{x}(t)$, rather than through $x(t)$. Furthermore, one can rewrite Eq. (1) as $S[x(t)]=(1/2)\iint dt_{1}dt_{2}\,\delta(t_{1}-t_{2})\dot{x}(t_{1})\dot{x}(t_{2})$, with a delta-function kernel. In the $x$-representation the kernel, $\partial^{2}_{t_{1}t_{2}}\,\delta(t_{1}-t_{2})$, is more singular and less convenient to work with.
  * (52) S. G. Samko, A. A. Kilbas and A. I. Marichev, Fractional integrals and derivatives: Theory and Applications, (Gordon and Breach, Paris, 1993).
  * (53) A. D. Polyanin and A. V. Manzhirov, Handbook of Integral Equations (Taylor & Francis Group, Boca Raton, FL, USA, 2008).
  * (54) See Supplemental Material at … for some derivations and technical details.
  * (55) T. Guggenberger, A. Chechkin, and R. Metzler, J. Phys. A: Math. Theor. 54, 29 (2021).
  * (56) B. Meerson and G. Oshanin, Phys. Rev. E 105, 064137 (2022).
  * (57) B. Meerson, Phys. Rev. E 100, 042135 (2019).
  * (58) B. Meerson, Phys. Rev. E 105, 034106 (2022).
  * (59) P. Lévy, University of California Publications in Statistics 1, 331 (1953).
  * (60) K. L. Sebastian, J. Phys. A: Math. Gen. 28, 4305 (1995).
  * (61) I. Calvo and R. Sanchez, J. Phys. A: Math. Theor. 41, 282002 (2008).
  * (62) I. Calvo, R. Sánchez and B. A. Carreras, J. Phys. A: Math. Theor. 42, 055003 (2009).
  * (63) H. S. Wio, J. Phys. A: Math. Theor. 46, 115005 (2013).
  * (64) T. Lundgren and D. Chiang, Quart. Appl. Math. 24, 303 (1967).