Integrable extensions of Adler’s map via Grassmann algebras

P. Adamopoulou†, S. Konstantinou-Rizos‡ and G. Papamikos⋆

†School of Mathematical & Computer Sciences, Heriot-Watt University, UK

‡ Centre of Integrable Systems, P.G. Demidov Yaroslavl State University,
Russia

⋆Department of Mathematical Sciences, University of Essex, UK

E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>

We study certain extensions of the Adler map on Grassmann algebras $\Gamma(n)$
of order $n$. We consider a known Grassmann-extended Adler map, and assuming
that $n=1$ we obtain a commutative extension of Adler’s map in six dimensions.
We show that the map satisfies the Yang–Baxter equation, admits three
invariants and is Liouville integrable. We solve the map explicitly, viewed as
a discrete dynamical system.

PACS numbers: 02.30.Ik

Mathematics Subject Classification: 15A75, 35Q53, 16T25, 17B80, 37J70

Keywords: Yang–Baxter maps, Grassmann algebras, Liouville integrability,

solutions of discrete dynamical systems, symplectic structures

##  1 Introduction

Since the late 80s, when the study of set-theoretical solutions to the
Yang–Baxter equation started [4, 24] (formally proposed by Drinfeld in 1990
[8]), a lot of important results have been obtained on the relation of the
Yang–Baxter equation with the field of integrable systems and various
algebraic structures. For example, the relation of the Yang–Baxter equation to
matrix refactorisation problems [18, 19, 25], the classification of
quadrirational Yang–Baxter maps [2, 22], the relation of Yang–Baxter maps with
Darboux transformations [15, 17, 21], with reflection maps in soliton theory
[5, 6], and more recently with the theory of braces [7, 23].

The increasing popularity of noncommutative structures over the past few
decades and their applications in the field of mathematical physics and of
integrable systems motivated the noncommutative (Grassmann) extension of
Darboux and Bäcklund transformations [12, 27, 28, 29]. This, in turn, gave
rise to the construction of Grassmann extended Yang–Baxter maps. In
particular, the first examples of such maps appeared in [11] in relation to
the nonlinear Schrödinger (NLS) equation and the derivative NLS equation, in
other works [14, 16] in relation to the KdV equation, as well as for the
Boussinesq equation in [13], where a Grassmann-extended Yang–Baxter together
with associated quad-graph systems was presented.

A definition for the complete integrability in the Liouville sense of
Yang–Baxter maps over Grassmann algebras is not yet available. However,
starting from a Grassmann Yang–Baxter map one can obtain a hierarchy of
commutative maps by fixing the number of generators of the underlying
Grassmann algebra. In this paper we study the Liouville integrability of such
a commutative map associated to a parametric Grassmann Adler map [16] by
restricting to the Grassmann algebra with one canonical generator.

This paper is organised as follows: In Section 2 we give all the necessary
definitions and notations that make this text self-contained. In particular,
we give the definition of Grassmann algebras, their grading structure, and
provide the definitions of superdeterminant and supertrace for matrices with
Grassmann elements. Also, we give the parametric set-theoretic Yang–Baxter
equation and define its solutions over sets with Grassmann variables. We call
such solutions parametric Grassmann (or supersymmetric) Yang–Baxter maps. We
further explain the relation of such maps to matrix refactorisation problems.
In Section 8 we provide a Grassmann-extended Adler map and invariant and anti-
invariant quantities. In Section 4 we give the definition of the Liouville
integrability for commutative maps, we present a six-dimensional commutative
extension of the Adler map, and we prove its Liouville integrability. Finally,
we discuss how the iterations of this map, seen as a discrete dynamical
system, can be solved exactly. In Section 5 we close we some concluding
remarks and ideas for future work.

##  2 Yang-Baxter equation and matrix refactorisation problems on Grassmann
algebras

A Grassmann algebra of order $n$ denoted by $\Gamma(n)$, over a field
$\mathbb{F}$ of characteristic zero, is an associative algebra with unit
$\mathbb{1}$ and $n$ generators $\theta_{i}$, $i=1,\ldots,n$, satisfying

$\theta_{i}\theta_{j}+\theta_{j}\theta_{i}=0\,.$ (1)

From (1), it follows that $\theta_{i}^{2}=0$ for all $i$ and that a generic
element of the Grassmann algebra $\Gamma(n)$ is of the form

$f=\sum_{s\geq 0}\sum_{i_{1}<\cdots<i_{s}}f_{i_{1}\ldots
i_{s}}\theta_{i_{1}}\cdots\theta_{i_{s}}.$ (2)

We say that those elements $\eqref{eq:graselem}$ of $\Gamma(n)$ that contain
sums of products of only even number of $\theta_{i}$’s are even, while those
that contain sums of products of only odd number of $\theta_{i}$’s are called
odd. We denote the subset of all even and odd elements by $\Gamma(n)_{0}$ and
$\Gamma(n)_{1}$, respectively. It follows that $\Gamma(n)$ admits a
$\mathbb{Z}_{2}$-gradation, i.e. it can be written as
$\Gamma(n)=\Gamma(n)_{0}\oplus\Gamma(n)_{1}$ with
$\Gamma(n)_{i}\Gamma(n)_{j}\subseteq\Gamma(n)_{(i+j)\mod 2}$. From the
$\mathbb{Z}_{2}-$grading it follows that $\Gamma(n)_{0}$ is a subalgebra of
$\Gamma(n)$ and that the even elements commute with all elements of
$\Gamma(n)$, while the odd elements anticommute with each other. Hence
$\Gamma(n)$ has the structure of a superalgebra with $\Gamma(n)_{0}$ being the
set of bosonic elements while $\Gamma(n)_{1}$ the set of fermionic elements.
In what follows, unless stated otherwise, elements of $\Gamma(n)_{0}$ will be
denoted by Latin letters, while elements of $\Gamma(n)_{1}$ by Greek letters,
with the exception of $\lambda$ which plays the role of a spectral
parameter111Objects used in this paper arise in the study of spectral problems
in soliton theory..

The notions of the determinant and the trace of a matrix with elements in
$\Gamma(n)$ are defined for square matrices, $M$, of the block-form
$M=\left(\begin{matrix}P&\Pi\\\ \Lambda&L\end{matrix}\right)$. The elements of
matrices $P$ and $L$ belong in $\Gamma(n)_{0}$, while those of $\Pi$ and
$\Lambda$ in $\Gamma(n)_{1}$. Matrices $\Pi$ and $\Lambda$ are not necessarily
square. The superdeterminant of $M$, denoted by $\operatorname{sdet}(M)$, is
defined by:

$\operatorname{sdet}(M)=\det(P-\Pi
L^{-1}\Lambda)\det(L^{-1})=\det(P^{-1})\det(L-\Lambda P^{-1}\Pi),$

where $\det(\cdot)$ is the usual determinant of a matrix, while the
supertrace,

$\operatorname{str}(M)=\operatorname{tr}(P)-\operatorname{tr}(L),$

where $\operatorname{tr}(\cdot)$ is the usual trace of a matrix. For more
details on Grassmann algebras, see [3, 10].

We are interested in solutions $S_{a,b}$ to the parametric Yang–Baxter
equation

$S^{12}_{a,b}\circ S^{13}_{a,c}\circ S^{23}_{b,c}=S^{23}_{b,c}\circ
S^{13}_{a,c}\circ S^{12}_{a,b}\,,$ (3)

where $a,b,c\in\mathbb{F}$ and $S_{a,b}:V_{n}^{k,l}\times
V_{n}^{k,l}\rightarrow V_{n}^{k,l}\times V_{n}^{k,l}$, with

$V_{n}^{k,l}:=\\{(\bm{x},\bm{\chi})\,|\,\bm{x}\in\Gamma(n)^{k}_{0},~{}\bm{\chi}\in\Gamma(n)^{l}_{1}\\}\,,$
(4)

and the indices $k,\,l$ are the number of even and odd Grassmann variables,
respectively. The maps $S_{a,b}^{ij}:V_{n}^{k,l}\times V_{n}^{k,l}\times
V_{n}^{k,l}\rightarrow V_{n}^{k,l}\times V_{n}^{k,l}\times V_{n}^{k,l}$,
$i,j=1,2,3$, $i\neq j$, are defined by the following relations

$S_{a,b}^{12}=S_{a,b}\times id\,,\quad S_{b,c}^{23}=id\times S_{b,c}\,,\quad
S_{a,c}^{13}=\pi^{12}S_{a,c}^{23}\pi^{12}\,,$

where $\pi^{12}$ is the permutation defined by

$\pi^{12}((\bm{x},\bm{\chi}),(\bm{y},\bm{\psi}),(\bm{z},\bm{\zeta}))=((\bm{y},\bm{\psi}),(\bm{x},\bm{\chi}),(\bm{z},\bm{\zeta}))\,.$

We call such solutions $S_{a,b}$ to the Yang–Baxter equation (3) parametric
Grassmann Yang–Baxter maps or parametric supersymmetric Yang–Baxter maps.

Following [25], we define a Lax matrix of the parametric Yang–Baxter map

$S_{a,b}((\bm{x},\bm{\chi}),(\bm{y},\bm{\psi}))=((\bm{u},\bm{\xi}),(\bm{v},\bm{\eta}))$

to be a matrix $\mathcal{L}_{a}((\bm{x},\bm{\chi});\lambda)$, depending on the
point $(\bm{x},\bm{\chi})\in V_{n}^{k,l}$, a parameter $a\in\mathbb{F}$ and a
spectral parameter $\lambda$, such that

$\mathcal{L}_{a}(\bm{u},\bm{\xi})\mathcal{L}_{b}(\bm{v},\bm{\eta})=\mathcal{L}_{b}(\bm{y},\bm{\psi})\mathcal{L}_{a}(\bm{x},\bm{\chi})\,,$
(5)

where we have suppressed the dependence on $\lambda$ for convenience.

The quantity
$\operatorname{str}\left(\mathcal{L}_{b}(\bm{y},\bm{\psi})\mathcal{L}_{a}(\bm{x},\bm{\chi})\right)$
constitutes a generating function of invariants for the map $S_{a,b}$. By an
invariant we mean a function $I$ of both even and odd variables such that

$I(\bm{u},\bm{\xi},\bm{v},\bm{\eta})=I(\bm{x},\bm{\chi},\bm{y},\bm{\psi}).$
(6)

It is also possible that the map admits anti-invariants, that is functions $J$
such that

$J(\bm{u},\bm{\xi},\bm{v},\bm{\eta})=-J(\bm{x},\bm{\chi},\bm{y},\bm{\psi}).$
(7)

The product of two different anti-invariants or the square of an anti-
invariant are invariants.

##  3 A Grassmann-extended Adler map

A non-commutative extension on Grassmann algebras of the Adler map [1]

$R_{a,b}:(x,y)\mapsto\left(y+\frac{a-b}{x+y},x+\frac{b-a}{x+y}\right)$ (8)

first appeared in [16], where the authors derived the map
$S_{a,b}:V_{n}^{1,2}\times V_{n}^{1,2}\rightarrow V_{n}^{1,2}\times
V_{n}^{1,2}$

$S_{a,b}:\left((x,\chi_{1},\chi_{2}),(y,\psi_{1},\psi_{2})\right)\mapsto\left((u,\xi_{1},\xi_{2}),(v,\eta_{1},\eta_{2})\right)$
(9)

with components given by

$\displaystyle x$ $\displaystyle\mapsto
u=y+\frac{a-b}{x+y-\chi_{1}\psi_{2}}\,,$ (10a) $\displaystyle\chi_{1}$
$\displaystyle\mapsto\xi_{1}=\psi_{1}-\frac{a-b}{x+y}\chi_{1}\,,$ (10b)
$\displaystyle\chi_{2}$ $\displaystyle\mapsto\xi_{2}=\psi_{2}\,,$ (10c)
$\displaystyle y$ $\displaystyle\mapsto
v=x-\frac{a-b}{x+y-\chi_{1}\psi_{2}}\,,$ (10d) $\displaystyle\psi_{1}$
$\displaystyle\mapsto\eta_{1}=\chi_{1}\,,$ (10e) $\displaystyle\psi_{2}$
$\displaystyle\mapsto\eta_{2}=\chi_{2}+\frac{a-b}{x+y}\psi_{2}\,.$ (10f)

The map $S_{a,b}$ given in (9) arises from the re-factorisation problem

$\mathcal{L}_{a}(u,\bm{\xi})\mathcal{L}_{b}(v,\bm{\eta})=\mathcal{L}_{b}(y,\bm{\psi})\mathcal{L}_{a}(x,\bm{\chi})$
(11)

of the following Lax matrix

$\mathcal{L}_{a}(x,\bm{\chi})=\begin{pmatrix}x&1&0\\\
x^{2}+\chi_{1}\chi_{2}-a+\lambda&x&\chi_{1}\\\ \chi_{2}&0&1\end{pmatrix},$
(12)

with $(x,\bm{\chi})=(x,\chi_{1},\chi_{2})\in V_{n}^{1,2}$. Matrix
$\mathcal{L}_{a}(x,\bm{\chi})$ was first derived in [28] in relation to the
Darboux transformation for a generalised super KdV system. We note that
$\mathcal{L}_{a}(x,\bm{\chi})$ has constant superdeterminant. Here we find
certain invariants and anti-invariants of the map (9)-(10) which was shown to
be a Yang–Baxter map in [16].

######  Proposition 3.0.1.

The map $S_{a,b}$ given in (9)-(10) admits the following invariants

$\displaystyle I_{1}$ $\displaystyle=x+y\,,$ (13) $\displaystyle I_{2}$
$\displaystyle=\chi_{1}\chi_{2}+\psi_{1}\psi_{2}\,,$ (14) $\displaystyle
I_{3}$
$\displaystyle=b\chi_{1}\chi_{2}+a\psi_{1}\psi_{2}+(x+y)(\chi_{1}\psi_{2}+\psi_{1}\chi_{2})\,,$
(15) $\displaystyle I_{4}$ $\displaystyle=\chi_{1}\psi_{1}\chi_{2}\psi_{2}\,.$
(16)

###### Proof.

Invariants $I_{i}$, $i=1,2$, follow from
$\operatorname{str}(\mathcal{L}_{b}(y,\psi_{1},\psi_{2})\mathcal{L}_{a}(x,\chi_{1},\chi_{2}))$
and were already found in [16]. The invariance of $I_{3}$ can be verified by
straightforward calculations. We notice that the quantities

$J_{1}=\chi_{1}\psi_{1}\quad\mbox{and}\quad J_{2}=\chi_{2}\psi_{2}$

are anti-invariants. Indeed,

$J_{1}\circ
S_{a,b}=\xi_{1}\eta_{1}=\left(\psi_{1}-\frac{a-b}{x+y}\chi_{1}\right)\chi_{1}=\psi_{1}\chi_{1}=-J_{1}\,,$

and similarly for $J_{2}$. The anti-invariants $J_{1}$, $J_{2}$ lead to the
invariant $I_{4}=J_{1}J_{2}=\chi_{1}\psi_{1}\chi_{2}\psi_{2}$. ∎

######  Remark 3.0.2.

In the bosonic limit, where all odd variables are zero and all even variables
are in $\mathbb{F}$, map (10) becomes the Adler map (8).

##  4 Adler map over $\Gamma(1)$

In this section we restrict the parametric Grassmann Yang–Baxter map $S_{a,b}$
given in (9)-(10) in the case where all variables are in $\Gamma(1)$. Hence,
we investigate the dynamical and integrability properties of the map
$S_{a,b}:V_{1}^{1,2}\times V_{1}^{1,2}\rightarrow V_{1}^{1,2}\times
V_{1}^{1,2}$. We first provide the definition of Liouville integrability for
maps in the commutative setting.

######  Definition 4.0.1.

Let $Y:(x_{1},...,x_{2N+M})\mapsto(u_{1},...,u_{2N+M})$,
$u_{i}=u_{i}(x_{1},...,x_{2N+M})$, $i=1,...,2N+M$, be a $(2N+M)$-dimensional
map. If

  * i.

there is a Poisson bracket $\\{\cdot,\cdot\\}$ such that
$\left\\{x_{i},x_{j}\right\\}$ has rank $2N$ and is invariant under the map
$Y$;

  * ii.

map $Y$ admits $N$ functionally independent invariants, $I_{i}$, such that
$\left\\{I_{i},I_{j}\right\\}=0$, $i,j=1,\ldots,N$;

  * iii.

there are $M$ Casimir functions, $C_{i}$, $i=1,\ldots,M$, i.e.
$\left\\{C_{i},f\right\\}=0$, for any arbitrary function
$f=f(x_{1},...,x_{2N+M})$, which are invariants of the map;

then $Y$ is said to be completely integrable or Liouville integrable [9, 20,
26].

In the case of the Grassmann algebra $\Gamma(1)$ with unit $\mathbb{1}$ and
generator $\theta$, such that $\theta^{2}=0$, even elements can be expressed
as $p\,\mathbb{1}$ with $p\in\mathbb{F}$, while for odd elements we have
$\rho\,\theta$, $\rho\in\mathbb{F}$.

######  Proposition 4.0.1.

The map (9)–(10) over $\Gamma(1)$ becomes a commutative map
$\mathbb{F}^{6}\rightarrow\mathbb{F}^{6}$

$S_{a,b}(x,\chi_{1},\chi_{2},y,\psi_{1},\psi_{2})=(u,\xi_{1},\xi_{2},v,\eta_{1},\eta_{2}),$
(17)

where

$\displaystyle x\mapsto u$ $\displaystyle=y+\frac{a-b}{x+y}\,,$ (18a)
$\displaystyle\chi_{1}\mapsto\xi_{1}$
$\displaystyle=\psi_{1}-\frac{a-b}{x+y}\chi_{1}\,,$ (18b)
$\displaystyle\chi_{2}\mapsto\xi_{2}$ $\displaystyle=\psi_{2}\,,$ (18c)
$\displaystyle y\mapsto v$ $\displaystyle=x-\frac{a-b}{x+y}\,,$ (18d)
$\displaystyle\psi_{1}\mapsto\eta_{1}$ $\displaystyle=\chi_{1}\,,$ (18e)
$\displaystyle\psi_{2}\mapsto\eta_{2}$
$\displaystyle=\chi_{2}+\frac{a-b}{x+y}\psi_{2}\,.$ (18f)

Map (17)-(18) is a Yang–Baxter map and admits the following invariants:

$\displaystyle I_{1}$ $\displaystyle=x+y\,,$ (19a) $\displaystyle I_{2}$
$\displaystyle=\chi_{1}\chi_{2}+\psi_{1}\psi_{2},$ (19b) $\displaystyle I_{3}$
$\displaystyle=b\chi_{1}\chi_{2}+a\psi_{1}\psi_{2}+(x+y)(\chi_{1}\psi_{2}+\psi_{1}\chi_{2}).$
(19c)

###### Proof.

We obtain the map in (18) by restricting all the variables in the map (10) in
the algebra $\Gamma(1)$. Then we obtain the expressions (18) by equating the
coefficients of $\mathbb{1}$ and $\theta$. The Yang–Baxter property can be
readily verified by straightforward substitution to the classical
(commutative) Yang–Baxter equation. Finally, the invariance of
$I_{1},I_{2},I_{3}$ under the map (18) can be directly verified. ∎

######  Remark 4.0.2.

It can also be verified that $J_{1}=\chi_{1}\psi_{1}$ and
$J_{2}=\chi_{2}\psi_{2}$ are not anti-invariants of the map (18). It also
follows that $I_{4}=J_{1}J_{2}$ given in Proposition (4.0.1) is not an
invariant.

######  Theorem 4.0.3.

Map (17)-(18) is completely integrable in the Liouville sense.

###### Proof.

The gradients $\nabla I_{i}$, $i=2,3$, are linearly independent, thus the
invariants $I_{i}$, $i=2,3$ are functionally independent. Moreover, the latter
are in involution $\\{I_{i},I_{j}\\}=0$ with respect to the following Poisson
bracket:

$\\{\chi_{1},\chi_{2}\\}=\\{\psi_{1},\psi_{2}\\}=1$ (20)

and all the rest fundamental brackets vanish. The rank of the associated
$6\times 6$ Poisson matrix $P$ is four. Additionally, for the following
quantities

$C_{1}=I_{1}=x+y,\quad C_{2}=ax+by+xy(x+y)$

we have that $C_{i}\circ S_{a,b}=C_{i}$, $i=1,2$, i.e. they are invariants
themselves, and also $(\nabla C_{i})^{t}P=0$, from which follows that the
latter are Casimir functions. Moreover, $\\{C_{i},I_{j}\\}=0$ for all possible
combinations. The Casimir $C_{2}$ can be constructed using the involutivity of
the Adler map (8) and averaging of the function
$(a+b)y^{2}+2(b-a)xy-(a+b)x^{2}$ over an orbit of the map. Finally, map
$S_{a,b}$ preserves the Poisson bracket, which completes the proof. ∎

######  Remark 4.0.4.

We observe that the commutative map (18) can be decomposed into Adler’s map
(8) and a linear map with matrix coefficients depending on the variables that
appear on Adler’s map. In particular the matrix coefficients are functions of
the Casimir $C_{1}$ and hence constant for a given orbit. This implies the
solvability of the map (18) if seen as discrete dynamical system.

##  5 Conclusions

In this paper we presented a commutative six-dimensional extension of Adler’s
map, and studied its integrability and dynamical properties. The map was
derived by restricting the parametric Grassmann Yang–Baxter map given in
(9)-(10) in the case of the $\Gamma(1)$ Grassmann algebra. We conjecture that
for each Grassmann algebra $\Gamma(n)$ the commutative map obtained, in a
similar way as in the case of $\Gamma(1)$, will be integrable, leading in such
a way to a hierarchy of ($3\times 2^{n}$)-dimensional commutative integrable
Yang–Baxter maps for $n=1,2,\ldots$. It is interesting to define complete
integrability in the Grassmann setting in such a way that the integrability of
the hierarchy of these commutative maps would directly follow. We aim to
extend the ideas presented in this paper to the case of entwining Yang–Baxter
maps.

## Acknowledgements

This paper is submitted for the proceedings of the second international
conference on Integrable Systems and Nonlinear Dynamics in Yaroslavl 19–23
October 2020.

This work started while S.K.-R. visited Heriot-Watt University in January
2019, and the University of Essex as an International Visiting Fellowship in
November 2019. S.K.-R.’s work was funded by the Russian Science Foundation
(project number 20-71-10110).

## References

  * [1] V. E. Adler, Recuttings of polygons, Funct. Anal. Its Appl. 27 141–143 (1993)
  * [2] V. E. Adler, A. I. Bobenko and Yu. B. Suris, Geometry of Yang–Baxter maps: pencils of conics and quadrirational mappings, Comm. Anal. Geom. 12 967–1007 (2004)
  * [3] F. A. Berezin, Introduction to Superanalysis, Springer, Dordrecht (1987)
  * [4] V. M. Bukhshtaber, Yang–Baxter mappings, Uspekhi Mat. Nauk 53 241–242 (1998) (English translation: Russian Math. Surveys 53 1343–1345)
  * [5] V. Caudrelier, N. Crampé and Q. C. Zhang, Set-theoretical reflection equation: classification of reflection maps, J. Phys. A: Math. Theor. 46 095203 (2013)
  * [6] V. Caudrelier and Q. C. Zhang, Yang–Baxter and reflection maps from vector solitons with a boundary, Nonlinearity 27 1081 (2014)
  * [7] A. Doikou, A. Smoktunowicz, From Braces to Hecke algebras & Quantum Groups, arXiv:1912.03091v2
  * [8] V. G. Drinfeld, On some unsolved problems in quantum group theory, in Quantum Groups, Springer, Leningrad (1990)
  * [9] A. P. Fordy, Periodic cluster mutations and related integrable maps, J. Phys. A 47 474003 (2014)
  * [10] L. Frappat, A. Sciarrino and P. Sorba, Dictionary on Lie Algebras and Superalgebras, Academic Press, Cornwall (2000)
  * [11] G. G. Grahovski, S. Konstantinou-Rizos and A. V. Mikhailov, Grassmann extensions of Yang–Baxter maps, J. Phys. A: Math. Theor.  49 145202 (2016)
  * [12] G. G. Grahovski and A. V. Mikhailov, Integrable discretisations for a class of nonlinear Schrödinger equations on Grassmann algebras, Phys. Lett. A 377 3254–3259 (2013)
  * [13] S. Konstantinou-Rizos, On the 3D consistency of a Grassmann extended lattice Boussinesq system, Nucl. Phys. B 951 114878 (2020)
  * [14] S. Konstantinou-Rizos and T. E. Kouloukas, A noncommutative discrete potential KdV lift, J. Math. Phys. 59 063506 (2018)
  * [15] S. Konstantinou-Rizos and A. V. Mikhailov, Darboux transformations, finite reduction groups and related Yang–Baxter maps, J. Phys. A: Math. Theor. 46 425201 (2013)
  * [16] S. Konstantinou-Rizos and A. V. Mikhailov, Anticommutative extension of the Adler map, J. Phys. A: Math. Theor. 49 30LT03 (2016)
  * [17] S. Konstantinou-Rizos and G. Papamikos, Entwining Yang–Baxter maps related to NLS type equations, J. Phys. A: Math. Theor. 52 485201 (2019)
  * [18] T. E. Kouloukas and V. G. Papageorgiou, Entwining Yang–Baxter maps and integrable lattices, Banach Center Publications 93 163–175 (2011)
  * [19] T. E. Kouloukas and V. G. Papageorgiou, Yang–Baxter maps with first-degree-polynomial $2\times 2$ Lax matrices J. Phys. A: Math. Theor. 42 404012 (2009)
  * [20] S. Maeda, Completely integrable symplectic mapping, Proc. Japan Acad. 63 198–200 (1987)
  * [21] A. V. Mikhailov, G. Papamikos and J. P. Wang, Darboux Transformation for the Vector sine-Gordon Equation and Integrable Equations on a Sphere, Lett. Math. Phys. 106 973–996 (2016)
  * [22] V. G. Papageorgiou, A. G. Tongas and A. P. Veselov, Yang–Baxter maps and symmetries of integrable equations on quad–graphs, J. Math. Phys. 47 083502 (2006)
  * [23] W. Rump, Modules over braces, Algebra Discrete Math., 2 127–137 (2006)
  * [24] E. K. Sklyanin, Classical limits of the $SU(2)$-invariant solutions of the Yang–Baxter equation, J. Math. Sci. 40 93–107 (1988)
  * [25] Yu. B. Suris and A. P. Veselov, Lax matrices for Yang–Baxter maps, J. Nonlinear Math. Phys. 10 223–230 (2003)
  * [26] A. P. Veselov, Integrable maps Russ. Math. Surv., 46 1–51 (1991)
  * [27] L.-L. Xue, D. Levi and Q.P. Liu, Supersymmetric KdV equation: Darboux transformation and discrete systems, J. Phys. A: Math. Theor. 46 502001 (2013)
  * [28] L.-L. Xue and Q. P. Liu, Bäcklund–Darboux transformations and discretizations of super KdV equation, Symmetry, Integrability and Geometry: Methods and Applications, 10 (2014)
  * [29] L.-L. Xue and Q. P. Liu, A supersymmetric AKNS problem and its Darboux-Bäcklund transformations and discrete systems, Stud. Appl. Math. 135 35–62 (2015)


# Convergence of the Discrete Minimum Energy Path††thanks:  This work was
supported by the National Key R&D Program of China (No. 2020YFA0712900). HC’s
work was also supported by National Natural Science Foundation of China (No.
NSFC11971066).

Xuanyu Liu111<EMAIL_ADDRESS>School of Mathematical Sciences,
Beijing Normal University, Beijing 100875, China., Huajie Chen222
<EMAIL_ADDRESS>School of Mathematical Sciences, Beijing Normal
University, Beijing 100875, China. and Christoph Ortner333<EMAIL_ADDRESS>Department of Mathematics, University of British Columbia, 1984 Mathematics
Road, Vancouver, BC, Canada V6T 1Z2.

###### Abstract

The minimum energy path (MEP) describes the mechanism of reaction, and the
energy barrier along the path can be used to calculate the reaction rate in
thermal systems. The nudged elastic band (NEB) method is one of the most
commonly used schemes to compute MEPs numerically. It approximates an MEP by a
discrete set of configuration images, where the discretization size determines
both computational cost and accuracy of the simulations. In this paper, we
consider a discrete MEP to be a stationary state of the NEB method and prove
an optimal convergence rate of the discrete MEP with respect to the number of
images. Numerical simulations for the transitions of some several proto-
typical model systems are performed to support the theory.

##  1 Introduction

In computational chemistry, the reaction mechanism is often represented by a
reaction path on the potential energy landscape, which is of great importance
to determine the thermodynamic and kinetic properties of chemical reactions.
The minimum energy path (MEP) has been most widely used to represent the
reaction path, giving the route that needs the least amount of work for the
system to undergo the transition. The MEP connects two local minimizers and
goes through one or more transition state (saddle), along which the energy
barrier (at the saddle) can be used to calculate the transition rate (see e.g.
[1, 6, 16]) of the reaction.

The most widely used techniques for finding the MEP are the nudged elastic
band (NEB) method [7, 8] and the string method [4, 5, 15]. They both
iteratively evolve a discrete path of images in projected steepest descent
directions, while keeping a smooth distribution of the images along the path.
The number of images along the path determines the accuracy and efficiency of
the simulation, especially when the energy functional is expensive to evaluate
(for example, when the system is described by quantum mechanical models). This
makes it particuarly important to obtain sharp error bounds of the MEP
approximation with respect to the number of images.

Notwithstanding the importance of such a convergence analysis, there is in
general very limited work devoted to the analysis of the MEP. In [3], the
authors give a comprehensive qualitative analysis of the evolution of the path
by the string method from a dynamical systems perspective, characterizing the
set of limiting curves and establishing conditions under which the limiting
curve is indeed a MEP. The techniques of this work are not suitable to perform
a numerical discretization analysis.

In [9], the authors investigated and analyzed the convergence of the string
method. Under some assumptions related to the critical points, they proposed a
uniform asymptotic stability of the MEP in the sense that any curve near the
MEP can be arbitrarily close to it in the Hausdorff distance under the
gradient decent dynamics. With this stability concept they showed that the
string method initialized in a neighborhood of the MEP can converge to an
arbitrarily small neighborhood of the limiting MEP as the number of images is
increased. However, for a curve near the MEP, this notion of stability does
not give an explicit relationship between its deviation to the MEP and its
force, hence it may not be able to provide an error bound with respect to the
number of images.

To prepare for a comprehensive and sharp convergence analysis of the NEB and
string methods, we began in [11] to establish a new framework for studying the
stability of MEPs. We showed that the displacement of a curve to the MEP can
be directly bounded by its force under appropriate weighted function space
norms. This result in particular directly implies the stability of MEPs under
perturbations of the potential energy landscape. Our result also suggests a
path towards studying the convergence of the discrete MEP, to which end we
still need an analogous stability for the discrete MEP (see Section 4.2).

The purpose of the present work is to obtain such a discrete stability result
and hence deduce sharp convergence rates of discrete MEPs to their continuous
limits. Under mild and natural assumptions on the limit MEP, we will indeed
obtain an optimal convergence rate. More precisely, we will show that (see
Theorem 2.1)

$\big{\|}\bar{\varphi}-\bar{\phi}_{M}\big{\|}\leq CM^{-1},$

where $\bar{\varphi}$ and $\bar{\phi}_{M}$ are the MEP and the discrete MEP
respectively, $\|\cdot\|$ denotes a discrete $C^{1}$-type-norm while $M$ is
the number of images used to discretize the path. To our best knowledge, this
is the first result that gives an explicit convergence rate for the MEP with
respect to the discretization parematers.

Outline. The rest of this paper is organized as follows. In Section 2, we
present the main results of this paper, including the convergence result and
numerical experiments. In Section 3, we give some conclusions. In Section 4,
we present detailed proofs for the convergence results.

Notation. Denote the Euclidean norm and $\ell^{\infty}$-norm of a vector by
$\lvert\cdot\rvert$ and $|\cdot|_{\infty}$, respectively. And
$\|\cdot\|_{\infty}$ is used to denote the matrix norm induced by
$|\cdot|_{\infty}$. For $b>a$, we denote by
$C\left([a,b];\mathbb{R}^{N}\right)$ the space of continuous curves in the
configuration space $\mathbb{R}^{N}$ with the norm
$\lVert\varphi\rVert_{C([a,b];\mathbb{R}^{N})}:=\sup_{\alpha\in[a,b]}|\varphi(\alpha)|_{\infty}$;
and $C^{1}\left([a,b];\mathbb{R}^{N}\right)$ the space of continuously
differentiable curves with the norm
$\lVert\varphi\rVert_{C^{1}([a,b];\mathbb{R}^{N})}:=\lVert\varphi\rVert_{C([a,b];\mathbb{R}^{N})}+\lVert\varphi^{\prime}\rVert_{C([a,b];\mathbb{R}^{N})}$.
Let $X$ and $Y$ be Banach spaces with the norm $\|\cdot\|_{X}$ and
$\|\cdot\|_{Y}$ respectively. We will denote by $\mathscr{L}(X,Y)$ the Banach
space of all linear bounded operators from $X$ to $Y$ with the operator norm
$\lVert\cdot\rVert_{\mathscr{L}(X,Y)}$. For a given functional $\mathscr{F}\in
C^{1}(X)$, we will denote its first variation by $\delta\mathscr{F}(x)v$ with
$v\in X$. We will use $C$ to denote a generic positive constant that may
change from one line to the next. In the error estimates, $C$ will always be
independent on the discretization parameters or the choice of test functions.
The dependencies of $C$ on model parameters (in our context, the energy
landscape) will normally be clear from the context or stated explicitly.

##  2 Main results

###  2.1 Minimum energy path

Let $E:\mathbb{R}^{N}\rightarrow\mathbb{R}$ be a potential energy functional
with $N\in\mathbb{N}$ the dimensionality of the configuration space, which
could encode atomic positions, the structure of crystal lattices, a
discretization of a function space, and many other examples. Throughout this
paper we assume that $E\in C^{4}(\mathbb{R}^{N})$. This regularity assumption
is required because we will need to control perturbations of the eigenvalues
and eigenvectors of the Hessian $\nabla^{2}E$ to be $C^{2}$ in our analysis.

Given an energy function $E$, we call $y\in\mathbb{R}^{N}$ a critical point if
$\nabla E(y)=0$. We call a critical point $y$ a strong local minimizer if the
Hessian $\nabla^{2}E(y)\in\mathbb{R}^{N\times N}$ is positive definite, and an
index-1 saddle point if $\nabla^{2}E(y)$ has exactly one negative eigenvalue
while all the other eigenvalues are positive. For the sake of brevity we will
omit the qualifiers “strong” and “index-1” and simply say “local minimizer”
and “saddle point”. We assume throughout that $E$ has at least two local
minimizers on the energy landscape denoted, respectively, by
$y^{A}_{M}\in\mathbb{R}^{N}$ and $y^{B}_{M}\in\mathbb{R}^{N}$.

A minimum energy path (MEP) is a curve $\varphi\in
C^{1}\big{(}[0,1];\mathbb{R}^{N}\big{)}$ connecting $y^{A}_{M}$ and
$y^{B}_{M}$ whose tangent is everywhere parallel to the gradient except at the
critical points. To give a rigorous definition of the MEP, we first define two
projection operators
$P_{v}~{},P_{v}^{\perp}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ with a given
a vector $v\in\mathbb{R}^{N}\setminus\\{0\\}$, by $\displaystyle
P_{v}:=vv^{T}/|v|^{2}$ and $P^{\perp}_{v}:=I-P_{v}$, where
$I\in\mathbb{R}^{N\times N}$ is the identity matrix. Define the admissible
class for “regular” curves connecting the minimizers $y^{A}_{M}$ and
$y^{B}_{M}$ as

$\displaystyle\mathscr{A}:=\Big{\\{}\varphi\in
C^{1}\big{(}[0,1];\mathbb{R}^{N}\big{)}:~{}\varphi(0)=y^{A}_{M},~{}\varphi(1)=y^{B}_{M},~{}\varphi^{\prime}(\alpha)\neq
0~{}~{}\forall~{}\alpha\in[0,1]\Big{\\}}.$

Then a MEP connecting the minimizers $y^{A}_{M}$ and $y^{B}_{M}$ is a solution
of the following problem: Find $\varphi\in\mathscr{A}$ such that

$\displaystyle\hskip 5.0ptP^{\perp}_{\varphi^{\prime}(\alpha)}\nabla
E\big{(}\varphi(\alpha)\big{)}={0}\hskip
21.52771pt~{}\forall~{}\alpha\in[0,1],$ (2.1a) $\displaystyle\hskip
5.0pt|\varphi^{\prime}(\alpha)|-L(\varphi)=0\hskip
30.1388pt~{}\forall~{}\alpha\in[0,1],$ (2.1b)

where $L:C^{1}([0,1];\mathbb{R}^{N})\rightarrow\mathbb{R}$ is the length
operator given by
$L(\varphi):=\int_{0}^{1}\lvert\varphi^{\prime}(s)\rvert~{}{\rm d}s$. Equation
(2.1b) enforces the curve to be parameterized by normalized arc length, which
removes the redundancy due to re-parameterization. With this additional
constraint one intuitively expects that the MEP is locally unique, at least
under some further natural assumptions that we now formalize:

If $\bar{\varphi}$ is a solution of (2.1), since $y^{A}_{M}$ and $y^{B}_{M}$
are local minimizers, there exists an $\bar{s}\in(0,1)$ with
$y_{S}=\bar{\varphi}(\bar{s})\in\mathbb{R}^{N}$ such that the energy
$E(y_{S})$ reaches the maximum along the MEP. This implies that $\nabla
E(y_{S})$ vanishes in the tangent direction $\bar{\varphi}^{\prime}(\bar{s})$
and thus it is a critical point. Since the energy $E(y_{S})$ is a maximum
along the tangent $\bar{\varphi}^{\prime}(\bar{s})$, we may generically expect
that the Hessian $\nabla^{2}E(y_{S})$ has at least one negative eigenvalue.
For the sake of simplicity of the analysis, we will assume throughout this
paper that

(A) $y_{S}$ is the only critical point between $y^{A}_{M}$ and $y^{B}_{M}$
along the MEP $\bar{\varphi}$. Moreover, $y^{A}_{M},y^{B}_{M}$ are strong
minimizers, while $y_{S}=\bar{\varphi}(\bar{s})$ is an index-1 saddle.

Although ${\bf(A)}$ is natural and will be satisfied by many (if not most)
MEPs one encounters in practice, there are also cases where this assumption
fails. For example, in [10] examples are given where there is more than one
index-1 saddle along an MEP. Our theory can be generalized to these cases by
adjusting the formulations, provided that all critical points along the MEP
satisfy certain stability conditions.

We observe by a direct calculation (see Lemma 4.1 for details) that, if
$\bar{\varphi}\in C^{2}\big{(}[0,1];\mathbb{R}^{N}\big{)}$ solves (2.1), then
$\bar{\varphi}^{\prime}(0)$, $\bar{\varphi}^{\prime}(\bar{s})$ and
$\bar{\varphi}^{\prime}(1)$ are eigenvectors of the Hessians
$\nabla^{2}E(y^{A}_{M})$, $\nabla^{2}E(y_{S})$ and $\nabla^{2}E(y^{B}_{M})$,
respectively. This implies that the MEP has to go through the critical points
in the direction of some eigenvector of the corresponding Hessian. The
following assumption formalizes the requirement that, to highest order, there
is a unique optimal path to exit the energy minimizer.

(B) Let $\sigma_{A}$ and $\sigma_{B}$ denote the eigenvalues associated,
respectively, with the eigenvectors
$\bar{\varphi}^{\prime}(0),\bar{\varphi}^{\prime}(1)$. We assume (i) that they
are the lowest eigenvalues of $\nabla^{2}E(y^{A}_{M})$ and
$\nabla^{2}E(y^{B}_{M})$; and (ii) that they are simple.

###  2.2 The discrete MEP and its convergence

We now describe the discretization of the MEP. A natural way to approximate
the continuous MEP curve is to use a discrete path of images connecting the
two minimizers. The images can best be thought of as interpolated by a spline.
The most widely used methods using this idea are the nudged elastic band (NEB)
method and the string method. In this paper, the will employ the NEB
formulation to define the discrete MEP, but our results can also be applied to
the string method provided that the same finite difference scheme as in the
NEB method is used to approximate the tangent (see Remark 2.2).

Let $M\in\mathbb{Z}_{+}$ be the number of images and $h:=\frac{1}{M}$ be the
corresponding meshsize. Let

$\displaystyle\mathscr{A}_{h}:=\Big{\\{}\phi_{h}:=\\{\phi_{h,k}\\}_{k=0}^{M}\in(\mathbb{R}^{N})^{M+1}:~{}\phi_{h,0}=y^{A}_{M},~{}\phi_{h,M}=y^{B}_{M},~{}{\rm
and}~{}\phi_{h,j}\neq\phi_{h,k}~{}\forall~{}j\neq k\Big{\\}}$

be the space for the discrete path connecting the two minimizers $y^{A}_{M}$
and $y^{B}_{M}$ with $M+1$ images.

Given an initial discrete path $\phi^{0}_{h}\in\mathscr{A}_{h}$, the the NEB
evolution equation is to find a $\phi_{h}\in
C^{1}\big{(}[0,+\infty);\mathscr{A}_{h}\big{)}$ such that

$\left\\{~{}\begin{aligned}
\dot{\phi}_{h}(t)&=\mathcal{F}\big{(}\phi_{h}(t)\big{)}\qquad\forall~{}t>0,\\\\[4.30554pt]
\phi_{h}(0)&=\phi^{0}_{h}.\end{aligned}\right.$ (2.2)

Here the driving force
$\mathcal{F}:\mathscr{A}_{h}\rightarrow(\mathbb{R}^{N})^{M+1}$ is defined by

$\big{(}\mathcal{F}(\phi_{h})\big{)}_{k}:=\left\\{~{}\begin{aligned} &0\hskip
71.13188pt{\rm for}~{}k=0,M,\\\\[4.30554pt] &-P_{\hat{\tau}_{k}}^{\perp}\nabla
E(\phi_{h,k})+ch^{-1}\Big{(}|\phi_{h,k}-\phi_{h,k+1}|-|\phi_{h,k}-\phi_{h,k-1}|\Big{)}\hat{\tau}_{k}\\\\[4.30554pt]
&~{}~{}\hskip 71.13188pt{\rm for}~{}k=1,\cdots,M-1,\end{aligned}\right.$ (2.3)

where $c$ is a parameter specifying the balance between the tangential force
and the projected downward pull from the curvature of the potential barrier
[8], and $\hat{\tau}_{k}$ is the normalized tangent direction at $\phi_{h,k}$
defined through an upwind scheme

$\displaystyle\hat{\tau}_{k}:=\frac{(D_{h}\phi_{h})_{k}}{|(D_{h}\phi_{h})_{k}|}\qquad{\rm
with}$
$\displaystyle(D_{h}\phi_{h})_{k}:=\left\\{\begin{array}[]{ll}\displaystyle\frac{\phi_{h,k+1}-\phi_{h,k}}{h}&{\rm
if}~{}k=0~{}{\rm
or}~{}E(\phi_{h,k+1})>E(\phi_{h,k})>E(\phi_{h,k-1}),\\\\[4.30554pt]
\displaystyle\frac{\phi_{h,k}-\phi_{h,k-1}}{h}&{\rm if}~{}k=M~{}{\rm
or}~{}E(\phi_{h,k-1})>E(\phi_{h,k})>E(\phi_{h,k+1}),\\\\[4.30554pt]
\displaystyle\frac{\phi_{h,k+1}-\phi_{h,k-1}}{2h}&{\rm
otherwise}.\end{array}\right.$ (2.7)

In (2.3), the force perpendicular to the tangent $\hat{\tau}_{k}$ is to push
the path towards the MEP, while the force parallel to the tangent
$\hat{\tau}_{k}$ is to enforce an approximate equidistribution of the images
along the path.

It is intuitive to expect that the NEB evolution path converges to the MEP in
the limit as $t\to\infty$ and $M\to\infty$. In the present paper, we focus on
the discretization error of the curve, i.e. the limit $M\to\infty$, while the
convergence of the evolution, i.e., the limit $t\to\infty$, will be considered
in our future work. Let $\bar{\phi}_{h}\in\mathscr{A}_{h}$ be a stationary
solution of the NEB evolution equation (2.2). Then $\bar{\phi}_{h}$ solves the
equation

$\displaystyle\mathcal{F}(\bar{\phi}_{h})={0},$ (2.8)

and we will call it a discrete MEP.

The following theorem is the main result of this paper, which states the
convergence and convergence rate of the discrete MEP to the continuous one.
The proof of this theorem is given in Section 4. We not only consider the
convergence of the path (under some discrete $C^{1}$-norm), but also the
convergence of the approximate energy barrier. For a given path
$\varphi\in\mathscr{A}$, its energy barrier can be written as

$\delta
E(\varphi):=\sup_{\alpha\in[0,1]}E\big{(}\varphi(\alpha)\big{)}-E\big{(}\varphi(0)\big{)}.$

Note that for a MEP $\bar{\varphi}$, the energy barrier is completely
determined by the energy at the saddle. Correspondingly, we can define the
energy barrier of a discrete path $\phi_{h}\in\mathscr{A}_{h}$ by

$\delta E(\phi_{h}):=\max_{0\leq k\leq M}E(\phi_{h,k})-E(\phi_{h,0}),$

with a bit of abuse of the notation $\delta E$.

######  Theorem 2.1.

Let $\bar{\varphi}\in C^{3}\big{(}[0,1];\mathbb{R}^{N}\big{)}$ be a solution
of (2.1). Assume that (A) and (B) are satisfied. Then for sufficiently small
$h$, there exists a solution $\bar{\phi}_{h}$ of (2.8) such that

$\displaystyle\max_{0\leq k\leq
M}\big{|}(D_{h}\bar{\phi}_{h})_{k}-\bar{\varphi}^{\prime}(kh)\big{|}+\max_{0\leq
k\leq M}\big{|}\bar{\phi}_{h,k}-\bar{\varphi}(kh)\big{|}\leq C_{\rm
p}h\qquad{\rm and}$ (2.9) $\displaystyle\quad\big{|}\delta
E(\bar{\varphi})-\delta E(\bar{\phi}_{h})\big{|}\leq C_{\rm e}h^{2},$ (2.10)

where $C_{\rm p}$ and $C_{\rm e}$ are positive constants depending only on $E$
and $\bar{\varphi}$,

######  Remark 2.1.

In our convergence analysis, we only consider numerical schemes that
approximate the tangent by first order finite differences (see the definition
of the upwind scheme in (2.2)). Therefore, the $\mathcal{O}(h)$ convergence
rate in Theorem 2.1 is optimal. Since the MEP itself is typically used only to
extract qualitative information this is sufficient for most applications. The
quantitative accuracy of the energy barrier is more important for the accurate
prediction of transition rates. If the accuracy provided by the MEP is
insufficient then one can easily improve on the prediction, e.g. using a dimer
method or a Krylov-Newton method, taking the maximum along the MEP as the
initial guess.

######  Remark 2.2.

Theorem 2.1 gives the convergence of the discrete MEP obtained by the NEB
method. For the other widely used method, the string method, the discrete path
$\phi_{h}$ is evolved by the driving force

$\big{(}\mathcal{F}^{\rm
string}(\phi_{h})\big{)}_{k}:=-P_{(D_{h}\phi_{h})_{k}}^{\perp}\nabla
E(\phi_{h,k})\qquad\quad{\rm for}~{}k=1,\cdots,M-1.$

Instead of using a force in tangent direction as that in the NEB method, a
redistribution of the images is performed at each step during the evolution
[5]. Then the discrete MEP defined by the stationary state of the string
method $\bar{\phi}_{h}^{\rm string}$ is such that the force $\mathcal{F}^{\rm
string}$ vanishes and the images are equally distributed along the path. This
implies $\mathcal{F}(\bar{\phi}_{h}^{\rm string})=0$, and thus
$\bar{\phi}_{h}^{\rm string}$ is exactly the same path as that obtained by the
NEB method. Therefore, the convergence result also holds for the discrete MEP
from the string method.

However, many implementations of the string method employ cubic spline
interpolatin to obtain the tangent and additional work would be required to
extend our analysis to those schemes.

###  2.3 Numerical experiments

We carry out some numerical experiments for the transitions of some typical
systems, to support and extend the theoretical results. All simulations are
implemented in open-source Julia packages JuLIP.jl [14], and performed on a PC
with Intel Core i7-CPU (2.2 GHz) with 32GB RAM. To test the decay of the
numerical errors, the results obtained by a very fine discretization (with up
to $1000$ images) are used as the reference MEP.

To support the theory developed in this paper, we test the convergence in the
discrete $C^{1}$-norm error (see (2.9)) and of the energy barrier error (see
(2.10)) of the discrete MEP. For completeness, we also test the convergence in
a discrete max-norm. These three error measures are, respectively, denoted by

$\displaystyle e_{h,C^{1}}:=\max_{0\leq k\leq
M}\big{|}(D_{h}\bar{\phi}_{h})_{k}-\bar{\varphi}^{\prime}(kh)\big{|}+\max_{0\leq
k\leq M}\big{|}\bar{\phi}_{h,k}-\bar{\varphi}(kh)\big{|},$ $\displaystyle
e_{h,E_{\rm b}}:=\big{|}\delta E(\bar{\varphi})-\delta
E(\bar{\phi}_{h})\big{|}\qquad{\rm and}$ $\displaystyle e_{h,C}:=\max_{0\leq
k\leq M}\big{|}\bar{\phi}_{h,k}-\bar{\varphi}(kh)\big{|}.$

Example 1. Consider a toy model with the potential
$E:\mathbb{R}\times\mathbb{R}_{+}\rightarrow\mathbb{R}$ with

$E(x,y)=(1-x^{2}-y^{2})^{2}+\frac{y^{2}}{x^{2}+y^{2}}\qquad y\geq 0.$

The contour lines of the energy landscape are shown in Figure 2.6. We see that
there are two minimizers located at $(-1,0)$ and $(1,0)$, respectively. We
present the decay of the numerical errors in Figure 2.6. We observe that the
convergence rate of the discrete $C^{1}$-norm error is $\mathcal{O}(h)$ and
that of the energy barrier error is $\mathcal{O}(h^{2})$, which agrees with
our theory. We further observe that the discrete $C$-norm error decays with
the same rate as the discrete $C^{1}$-norm error, i.e., there is no improved
rate when employing a weaker norm. This suggests that the spaces we employed
in the analysis are in fact natural.

In the above energy functional, the eigenvalues of $\nabla^{2}E(y^{A}_{M})$
are 2 and 8, and $\sigma_{A}=2$ is simple and the lowest eigenvalue. As
comparisons, we test two other energy functionals with similar form but not
satisfying the condition (B). First, let
$E_{1}:\mathbb{R}\times\mathbb{R}_{+}\rightarrow\mathbb{R}$ be given by

$E_{1}(x,y)=\frac{1}{4}(1-x^{2}-y^{2})^{2}+\frac{y^{2}}{x^{2}+y^{2}}\qquad
y\geq 0.$

The Hessian $\nabla^{2}E_{1}(y^{A}_{M})$ has an eigenvalue 2 with multiplicity
2, so $\sigma_{A}=2$ is the lowest but a degenerated eigenvalue of
$\nabla^{2}E_{1}(y^{A}_{M})$. Note that the exact MEP of this case is
explicitly known as the upper and lower branches of the unit circle, from
which we can calculate the exact numerical errors. We present the contour
lines of the energy landscape and the decay of numerical errors in Figure 2.6
and 2.6 respectively. Secondly, let
$E_{2}:\mathbb{R}\times\mathbb{R}_{+}\rightarrow\mathbb{R}$ be given by

$E_{2}(x,y)=\frac{1}{8}(1-x^{2}-y^{2})^{2}+\frac{y^{2}}{x^{2}+y^{2}}\qquad
y\geq 0.$

The eigenvalues of $\nabla^{2}E_{2}(y^{A}_{M})$ are 2 and 1, in which
$\sigma_{A}=2$ is simple but not the lowest eigenvalue. The exact MEP of this
case is also known as the upper and lower branches of the unit circle. We show
the contour lines of the energy landscape and the decay of numerical errors in
Figure 2.6 and 2.6 respectively. In the last two examples, we see that the
discrete MEP also converge in discrete $C^{1}$-norm, but the convergence rates
are significantly slower and it is only $O(h^{1/2})$ for the last case. This
indicates that the condition (B) is really necessary to obtain Theorem 2.1.
But surprisingly, we find that the discrete $C$-norm error and the energy
barrier error still decay as $O(h)$ and $O(h^{2})$, respectively. The methods
of our present work are insufficient to explain these effects; obtaining sharp
rates in the linearly unstable regime will require further work.

Figure 2.1: (Example 1) The contour lines of the energy landscape for $E$,
with two minimizers and the MEP (red line).

Figure 2.2: (Example 1) The convergence of the discrete MEP ($\sigma_{A}$ is
simple and lowest eigenvalue).

Figure 2.3: (Example 1) The contour lines of the energy landscape for $E_{1}$,
with two minimizers and the MEP (red line).

Figure 2.4: (Example 1) The convergence of the discrete MEP ($\sigma_{A}$ is a
degenerated eigenvalue).

Figure 2.5: (Example 1) The contour lines of the energy landscape for $E_{2}$,
with two minimizers and the MEP (red line).

Figure 2.6: (Example 1) The convergence of the discrete MEP ($\sigma_{A}$ is
not the lowest eigenvalue).

Example 2. Next, we consider the Muller potential [12]
$E:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ with

$E(x,y)=\sum_{k=1}^{4}T_{k}\exp\Big{(}a_{k}(x-\bar{x}_{k})^{2}+b_{k}(x-\bar{x}_{k})(y-\bar{y}_{k})+c_{k}(y-\bar{y}_{k})^{2}\Big{)},$

where the parameters are $T=(-200,-100,-170,-15)^{\rm T}$,
$a=(-1,-1,-6.5,0.7)^{\rm T}$, $b=(0,0,11,0.6)^{\rm T}$,
$c=(-10,-10,-6.5,0.7)^{\rm T}$, $\bar{x}=(1,0,-0.5,-1)^{\rm T}$ and
$\bar{y}=(0,0.5,1.5,1)^{\rm T}$. We present the contour lines of energy
landscape in Figure 2.8 and the decay of numerical errors in Figure 2.8. We
again observe an $\mathcal{O}(h)$ decay for the discrete $C^{1}$-norm and
$C$-norm errors, and $\mathcal{O}(h^{2})$ for the energy barrier error.

Figure 2.7: (Example 2) The contour lines of the energy landscape, with two
minimizers and the MEP (red line).

Figure 2.8: (Example 2) Convergence of the discrete MEP.

Figure 2.9: (Example 3) The configurations of two minimizers.

Figure 2.10: (Example 3) The contour lines of energy landscape, with two
minimizers and the MEP (red line). Since the energy explodes when the two
atoms collide, we only show the contour lines with energy less than $0$.

Figure 2.11: (Example 3) Energy path of MEP

Figure 2.12: (Example 3) Convergence of the discrete MEP.

Example 3. Consider a two dimensional periodic system with two atoms lying in
each unit cell. In the “first” cell containing the origin, we fix one atom at
the origin and put the other at the position $(x,y)\in\mathbb{R}^{2}$, and
then construct a corresponding lattice
$\Lambda:=\left(\begin{smallmatrix}2x&x\\\
0&y\end{smallmatrix}\right)\mathbb{Z}^{2}$ (see Figure 2.12). The interaction
between any two atoms in the lattice $\Lambda$ is modelled by the Lennard-
Jones pair potential

$\phi(r)=4\epsilon_{0}\left(\left(\frac{\sigma_{0}}{r}\right)^{12}-\left(\frac{\sigma_{0}}{r}\right)^{6}\right)\phi_{\rm
cut}(r)$

with parameters $\epsilon_{0}=1,~{}\sigma_{0}=1$ and a cubic spline cut-off
function [2] $\phi_{\rm cut}$ satisfying $\phi_{\rm cut}(r)=1$ for
$r\in[0,1.9]$ and $\phi_{\rm cut}(r)=0$ for $r\geq 2.7$. Then the “averaged”
energy of the system (i.e. energy per cell) is the sum of the site energy of
the two atoms located at origin and $(x,y)$,

$E(x,y)=\frac{1}{2}\sum_{(x^{\prime},y^{\prime})\in\Lambda}\phi\big{(}|(x^{\prime},y^{\prime})-(0,0)|\big{)}+\frac{1}{2}\sum_{(x^{\prime},y^{\prime})\in\Lambda}\phi\big{(}|(x^{\prime},y^{\prime})-(x,y)|\big{)}.$

We present the contour lines of the energy landscape in Figure 2.12, together
with the two minimizers and the MEP. We further present the energy path of the
MEP in Figure 2.12. The decay of the numerical errors are shown in Figure
2.12, from which we observe an $\mathcal{O}(h)$ convergence rate for the
discrete $C^{1}$-norm and discrete $C$-norm errors, and an
$\mathcal{O}(h^{2})$ convergence rate for the energy barrier error.

##  3 Conclusions

This paper studies the convergence of the discrete MEPs that are given by the
stationary states of the NEB methods. We establish an optimal convergence rate
of the discrete MEP with respect to the number of images along the path. For
numerical methods that find the MEP by evolving the path between local
minimizers (say NEB methods and string methods), an important remaining issue
is the convergence with respect to the evolution time, which will be studied
in our other work.

##  4 Proofs: Convergence of the discrete MEP

In this section, we will provide a proof of Theorem 2.1. We will first propose
abstract forms for the MEP equation (2.1) and the corresponding
discretization; then give a stability result of the discrete MEP; and finally
derive the convergence by connecting the abstract discrete MEP equation with
(2.8).

###  4.1 The abstract MEP and discrete MEP equations

Define an abstract MEP operator $\mathscr{F}:\mathscr{A}\rightarrow
C\big{(}[0,1];\mathbb{R}^{N}\big{)}$ with

$\displaystyle\hskip-8.5359pt\mathscr{F}(\varphi)(\alpha):=P^{\perp}_{\varphi^{\prime}(\alpha)}\nabla
E\big{(}\varphi(\alpha)\big{)}-\alpha(\alpha-1)(\alpha-\bar{s})\big{(}|\varphi^{\prime}(\alpha)|-L(\varphi)\big{)}\frac{\varphi^{\prime}(\alpha)}{\lvert\varphi^{\prime}(\alpha)\rvert}$
(4.1)

for $\alpha\in[0,1]$. Then we can bring in an abstract form of (2.1): Find
$\varphi\in\mathscr{A}$ such that

$\mathscr{F}(\varphi)=0.$ (4.2)

Comparing the components of $\mathscr{F}$ with (2.1), we see that
$\bar{\varphi}$ solves (2.1) if and only if it is a solution of (4.2).

Consider a discrete path with $M+1$ ($M\in\mathbb{Z}_{+}$) images. Let
$h=\frac{1}{M}$ be the corresponding mesh size. We denote by
$\alpha_{k}:=kh~{}(k=0,\cdots,M)$ the nodes of the discrete path. Since
$E(y_{S})$ reaches the maximum along the MEP, there exists a
$k_{\bar{s}}\in\mathbb{Z}_{+}$ such that
$\alpha_{k_{\bar{s}}}\in(\bar{s}-h,\bar{s}+h)$ and $\displaystyle
E\big{(}\bar{\varphi}(\alpha_{k_{\bar{s}}})\big{)}=\max_{0\leq k\leq
M}E(\bar{\varphi}(\alpha_{k}))$. Note that $\alpha_{k_{\bar{s}}}$ is to
approximate the index $\bar{s}$ and $\displaystyle\lim_{h\rightarrow
0^{+}}\alpha_{k_{\bar{s}}}=\bar{s}$.

To deduce a discretization of (4.2), we first introduce two first order
difference operators. Define
$\tilde{D}_{h},~{}\hat{D}_{h}:(\mathbb{R}^{N})^{M+1}\rightarrow(\mathbb{R}^{N})^{M+1}$
with

$\displaystyle(\tilde{D}_{h}\varphi_{h})_{k}$
$\displaystyle:=\left\\{\begin{array}[]{ll}\displaystyle\frac{\varphi_{h,k+1}-\varphi_{h,k}}{h}&{\rm
for}~{}k=0,\cdots,k_{\bar{s}}-1\\\\[4.30554pt]
\displaystyle\frac{\varphi_{h,k+1}-\varphi_{h,k-1}}{2h}&{\rm
for}~{}k=k_{\bar{s}}\\\\[4.30554pt]
\displaystyle\frac{\varphi_{h,k}-\varphi_{h,k-1}}{h}&{\rm
for}~{}k=k_{\bar{s}}+1,\cdots,M\end{array}\right.\qquad\qquad{\rm and}$ (4.6)
$\displaystyle\displaystyle(\hat{D}_{h}\varphi_{h})_{k}$
$\displaystyle:=\left\\{\begin{array}[]{ll}\displaystyle\frac{\varphi_{h,k+1}-\varphi_{h,k}}{h}&{\rm
for}~{}k=0\\\\[4.30554pt]
\displaystyle\frac{\varphi_{h,k}-\varphi_{h,k-1}}{h}&{\rm
for}~{}k=1,\cdots,M\end{array}\right..$ (4.9)

Let

$Y_{h}:=\big{\\{}f_{h}=\\{f_{h,k}\\}_{k=0}^{M}\in(\mathbb{R}^{N})^{M+1}:~{}f_{h,0}=f_{h,M}=0\big{\\}}$

equipped with the norm

$\|f_{h}\|_{Y_{h}}:=\max_{1\leq k\leq
M-1}\left|\frac{f_{h,k}}{\alpha_{k}(\alpha_{k}-1)}\right|_{\infty}+\max_{\begin{subarray}{c}0\leq
k\leq M,\\\\[1.42271pt] k\neq
k_{\bar{s}}\end{subarray}}\left|\frac{f_{h,k}-f_{h,k_{\bar{s}}}}{\alpha_{k}-\alpha_{k_{\bar{s}}}}\right|_{\infty}.$
(4.10)

Define $\mathscr{F}_{h}:\mathscr{A}_{h}\rightarrow Y_{h}$ by

$\displaystyle\big{(}\mathscr{F}_{h}(\varphi_{h})\big{)}_{k}$
$\displaystyle:=P_{(\tilde{D}_{h}\varphi_{h})_{k}}^{\perp}\nabla
E(\varphi_{h,k})-\alpha_{k}(\alpha_{k}-1)(\alpha_{k}-\alpha_{k_{\bar{s}}+\frac{1}{2}})\big{(}|(\hat{D}_{h}\varphi_{h})_{k}|-L_{h}(\varphi_{h})\big{)}\frac{(\tilde{D}_{h}\varphi_{h})_{k}}{|(\tilde{D}_{h}\varphi_{h})_{k}|}$
(4.11)

for $k=0,\cdots,M$, where
$\alpha_{k_{\bar{s}}+\frac{1}{2}}:=(k_{\bar{s}}+\frac{1}{2})h$ and
$L_{h}:\mathscr{A}_{h}\rightarrow\mathbb{R}$ is a discretization of the length
operator $L$

$\displaystyle L_{h}(\varphi_{h})$
$\displaystyle:=h\sum_{k=1}^{M}|(\hat{D}_{h}\varphi_{h})_{t}|=\sum_{k=1}^{M}|\varphi_{h,k}-\varphi_{h,k-1}|.$
(4.12)

We can then provide a discretization of (4.2): Find
$\varphi_{h}\in\mathscr{A}_{h}$ such that

$\mathscr{F}_{h}(\varphi_{h})=0.$ (4.13)

Then the solution $\bar{\varphi}_{h}$ of (4.13) will be taken as an
approximation of MEP $\bar{\varphi}$. We observe from (4.9), (4.12) and
(4.13)that

$\displaystyle\lvert\bar{\varphi}_{h,k}-\bar{\varphi}_{h,k-1}\rvert=\lvert\bar{\varphi}_{h,k}-\bar{\varphi}_{h,k+1}\rvert\qquad\forall~{}k=1,\cdots,M-1.$
(4.14)

This implies that the images on $\bar{\varphi}_{h}$ are equally distributed.

To give the first variation of $\mathscr{F}$ at $\bar{\varphi}$, we introduce
a function measuring the gradient of MEP, which will be heavily used in our
analysis. Let

$\bar{\lambda}(\alpha):=\nabla
E\big{(}\bar{\varphi}(\alpha)\big{)}^{T}\frac{\bar{\varphi}^{\prime}(\alpha)}{\lvert\bar{\varphi}^{\prime}(\alpha)\rvert^{2}},\qquad\alpha\in[0,1].$
(4.15)

We see from (2.1a) that

$\nabla
E\big{(}\bar{\varphi}(\alpha)\big{)}=\bar{\lambda}(\alpha)\bar{\varphi}^{\prime}(\alpha)\qquad{\rm
for}~{}\alpha\in[0,1],$ (4.16)

The following lemma states the properties of $\bar{\lambda}$ and it can be
seen form [11, Lemma 4.1]. We also observe from the following lemma that
$\sigma_{A}=\bar{\lambda}^{\prime}(0)$ and
$\sigma_{B}=\bar{\lambda}^{\prime}(1)$.

######  Lemma 4.1.

[11, Lemma 4.1] Let $\bar{\varphi}\in C^{2}\big{(}[0,1];\mathbb{R}^{N}\big{)}$
be the solution of (2.1). If (A) is satisfied, then

  * (i)

$\big{(}\bar{\lambda}^{\prime}(0),\bar{\varphi}^{\prime}(0)\big{)}$,
$\big{(}\bar{\lambda}^{\prime}(\bar{s}),\bar{\varphi}^{\prime}(\bar{s})\big{)}$
and $\big{(}\bar{\lambda}^{\prime}(1),\bar{\varphi}^{\prime}(1)\big{)}$ are
eigenpairs of the Hessians $\nabla^{2}E(y^{A}_{M})$, $\nabla^{2}E(y_{S})$ and
$\nabla^{2}E(y^{B}_{M})$, respectively;

  * (ii)

$\bar{\lambda}^{\prime}(0)>0,~{}\bar{\lambda}^{\prime}(1)>0$ and
$\bar{\lambda}^{\prime}(\bar{s})<0$; and

  * (iii)

there exist positive constants $\underline{c},\bar{c}$ depending only on
$\bar{\varphi}$, such that

$\underline{c}\leq\left\lvert\frac{\bar{\lambda}(\alpha)}{\alpha(\alpha-1)(\alpha-\bar{s})}\right\rvert\leq\bar{c}\qquad\forall~{}\alpha\in(0,\bar{s})\cup(\bar{s},1).$
(4.17)

Let

$\displaystyle X$ $\displaystyle:=\\{\psi\in
C^{1}\big{(}[0,1];\mathbb{R}^{N}\big{)}:~{}\psi(0)=\psi(1)=0\\}\hskip
56.9055pt{\rm and}$ $\displaystyle X_{h}$
$\displaystyle:=\Big{\\{}\psi_{h}=\\{\psi_{h,k}\\}_{k=0}^{M}\in(\mathbb{R}^{N})^{M+1}:~{}\psi_{h,0}=\psi_{h,M}=0\Big{\\}}$

equipped with the norms $\|\psi\|_{X}:=\|\psi\|_{C^{1}([0,1];\mathbb{R}^{N})}$
and
$\|\psi_{h}\|_{X_{h}}:=|\hat{D}_{h}\psi_{h}|_{\infty}+|\psi_{h}|_{\infty}$,
respectively. Note that $X_{h}$ is the same space as $Y_{h}$ but equipped with
a different norm. Using (4.16), we obtain by a direct computation that the
first variation $\delta\mathscr{F}(\bar{\varphi}):X\rightarrow
C\big{(}[0,1];\mathbb{R}^{N}\big{)}$ is given by

$\displaystyle\delta\mathscr{F}(\bar{\varphi})\psi(\alpha)$
$\displaystyle=P^{\perp}_{\bar{\varphi}^{\prime}(\alpha)}\Big{(}\nabla^{2}E\big{(}\bar{\varphi}(\alpha)\big{)}\psi(\alpha)-\bar{\lambda}(\alpha)\psi^{\prime}(\alpha)\Big{)}$
$\displaystyle\hskip-51.21504pt+\alpha(\alpha-1)(\alpha-\bar{s})\left(\frac{\bar{\varphi}^{\prime}(s)^{T}\psi^{\prime}(s)}{|\bar{\varphi}^{\prime}(s)|}-\int_{0}^{1}\frac{\bar{\varphi}^{\prime}(s)^{T}\psi^{\prime}(s)}{|\bar{\varphi}^{\prime}(s)|}~{}{\rm
d}s\right)\frac{\bar{\varphi}^{\prime}(\alpha)}{\lvert\bar{\varphi}^{\prime}(\alpha)\rvert}\qquad{\rm
for}~{}\alpha\in[0,1].$ (4.18)

For $\varphi_{h}\in\mathscr{A}_{h}$, we can obtain the first variation
$\delta\mathscr{F}_{h}(\varphi_{h}):X_{h}\rightarrow Y_{h}$ by a direct
calculation that for $k=0,\cdots,M$,

$\displaystyle\hskip-17.07182pt\big{(}\delta\mathscr{F}_{h}(\varphi_{h})\psi_{h}\big{)}_{k}=P_{(\tilde{D}_{h}\varphi_{h})_{k}}^{\perp}\left(\nabla^{2}E(\varphi_{h,k})\psi_{h,k}-\frac{\nabla
E(\varphi_{h,k})^{T}(\tilde{D}_{h}\varphi_{h})_{k}}{|(\tilde{D}_{h}\varphi_{h})_{k}|^{2}}(\tilde{D}_{h}\psi_{h})_{k}\right)$
$\displaystyle-\frac{\big{(}P_{(\tilde{D}_{h}\varphi_{h})_{k}}^{\perp}\nabla
E(\varphi_{h,k})\big{)}^{T}(\tilde{D}_{h}\psi_{h})_{k}}{|(\tilde{D}_{h}\varphi_{h})_{k}|}\frac{(\tilde{D}_{h}\varphi_{h})_{k}}{|(\tilde{D}_{h}\varphi_{h})_{k}|}$
$\displaystyle+\alpha_{k}(\alpha_{k}-1)(\alpha_{k}-\alpha_{k_{\bar{s}}+\frac{1}{2}})\left(\frac{(\hat{D}_{h}\psi_{h})_{k}^{T}(\hat{D}_{h}\varphi_{h})_{k}}{|(\hat{D}_{h}\varphi_{h})_{k}|}-h\sum_{t=1}^{M}\frac{(\hat{D}_{h}\psi_{h})_{t}^{T}(\hat{D}_{h}\varphi_{h})_{t}}{|(\hat{D}_{h}\varphi_{h})_{t}|}\right)\frac{(\tilde{D}_{h}\varphi_{h})_{k}}{|(\tilde{D}_{h}\varphi_{h})_{k}|}$
$\displaystyle+\alpha_{k}(\alpha_{k}-1)(\alpha_{k}-\alpha_{k_{\bar{s}}+\frac{1}{2}})\big{(}|(\hat{D}_{h}\varphi_{h})_{k}|-L_{h}(\varphi_{h})\big{)}\frac{P_{(\tilde{D}_{h}\varphi_{h})_{k}}^{\perp}(\tilde{D}_{h}\psi_{h})_{k}}{|(\tilde{D}_{h}\varphi_{h})_{k}|}.$
(4.19)

###  4.2 Stability of the discrete MEP

Define a projection operator
$\Pi_{h}:C([0,1];\mathbb{R}^{N})\rightarrow(\mathbb{R}^{N})^{M+1}$ by

$(\Pi_{h}\varphi)_{k}:=\varphi(kh),\qquad k=0,\cdots,M.$

The following lemma states the stability of the discrete MEP, which is a key
ingredient in the convergence analysis.

######  Lemma 4.2.

Let $\bar{\varphi}\in C^{3}\big{(}[0,1];\mathbb{R}^{N}\big{)}$ be a solution
of (2.1). Assume that (A) and (B) are satisfied. Then for sufficiently small
$h$, there exists a constant $C>0$ depending only on $E$ and $\bar{\varphi}$
such that

$\|\psi_{h}\|_{X_{h}}\leq
C\|\delta\mathscr{F}_{h}(\Pi_{h}\bar{\varphi})\psi_{h}\|_{Y_{h}}\qquad\forall~{}\psi_{h}\in
X_{h}.$ (4.20)

Before proving this lemma, let

$\displaystyle\hat{X}_{h}$ $\displaystyle:=\Big{\\{}\hat{\psi}_{h}\in
X:~{}\hat{\psi}_{h}|_{[\alpha_{k-1},\alpha_{k}]}\text{ is a cubic polynomial
for}~{}1\leq k\leq M$ $\displaystyle\hskip 56.9055pt\quad{\rm
and}~{}\hat{\psi}_{h}~{}{\rm
satisfying}~{}~{}\hat{\psi}^{\prime}_{h}(\alpha_{k})=(\tilde{D}_{h}\Pi_{h}\hat{\psi}_{h})_{k}\quad{\rm
for}~{}0\leq k\leq M\Big{\\}}$ (4.21)

equipped with the norm $\|\cdot\|_{X}$. We see that dim$(\hat{X}_{h})=N(M-1)$
and any $\hat{\psi}_{h}\in\hat{X}_{h}$ can be determined by their values at
the nodes $\alpha_{k}~{}(1\leq k\leq M-1)$. Thus
$\Pi_{h}:\hat{X}_{h}\rightarrow X_{h}$ is an isomorphism. Denote its inverse
operator by $P_{h}:X_{h}\rightarrow\hat{X}_{h}$ satisfying

$P_{h}\psi_{h}(\alpha_{k})=\psi_{h,k}\qquad{\rm for}~{}k=0,\cdots,M.$ (4.22)

Given $\hat{\psi}_{h}\in\hat{X}_{h}$, let

$f_{h}=\Pi_{h}\delta\mathscr{F}(\bar{\varphi})\hat{\psi}_{h}.$ (4.23)

To show Lemma 4.2, we will first bring in a basis set of $\hat{X}_{h}$ to
derive a spectral representation of (4.23) and obtain some linear systems.
Then by analyzing the coefficient matrics we will show that the operator
$\Pi_{h}\delta\mathscr{F}(\bar{\varphi}):\hat{X}_{h}\rightarrow Y_{h}$ is
invertible and
$\left\|\big{(}\Pi_{h}\delta\mathscr{F}(\bar{\varphi})\big{)}^{-1}\right\|_{\mathscr{L}(Y_{h},\hat{X}_{h})}$
is uniformly bounded with respect to $h$. Finally, Lemma 4.2 follows from a
perturbation argument.

For $\alpha\in[0,1]$, we consider the eigenvalue problem

$\Big{(}P^{\perp}_{\bar{\varphi}^{\prime}(\alpha)}\nabla^{2}E\big{(}\bar{\varphi}(\alpha)\big{)}P^{\perp}_{\bar{\varphi}^{\prime}(\alpha)}\Big{)}\xi_{i}(\alpha)=z_{i}(\alpha)\xi_{i}(\alpha)\qquad{\rm
for}~{}i=0,1,\cdots,N-1$ (4.24)

with $\\{z_{i}(\alpha)\\}_{i=0}^{N-1}$ the eigenvalues and
$\\{\xi_{i}(\alpha)\\}_{i=0}^{N-1}$ the corresponding eigenfunctions. From
[11], we can order the functions $z_{i}(\alpha)$ such that $z_{0}(\bar{s})\leq
z_{1}(\bar{s})\leq\cdots\leq z_{N-1}(\bar{s})$ and $z_{i}\in
C^{2}([0,1];\mathbb{R}),~{}\xi_{i}\in C^{2}\left([0,1];\mathbb{R}^{N}\right)$
for $i=0,1,\cdots,N-1$. Moreover, we have $z_{0}(\alpha)\equiv 0$,
$\xi_{0}(\alpha)\equiv\bar{\varphi}^{\prime}(\alpha)/\lvert\bar{\varphi}^{\prime}(\alpha)\rvert$
and

$\displaystyle z_{j}(\alpha)>0\qquad{\rm
for}~{}j=1,\cdots,N-1,~{}\alpha=0,\bar{s},1.$ (4.25)

Let $g_{h,k,j}\in\hat{X}_{h}~{}(1\leq k\leq M-1,~{}0\leq j\leq N-1)$ be a set
of functions satisfying

$\displaystyle
g_{h,k,j}(\alpha_{t})=\delta_{kt}~{}\xi_{j}(\alpha_{t})\qquad{\rm
for}~{}t=1,\cdots,M-1.$ (4.26)

Note that $\big{\\{}g_{h,k,j}\big{\\}}$ gives a complete basis set for the
finite dimensional space $\hat{X}_{h}$, since any function
$\hat{\psi}_{h}\in\hat{X}_{h}$ can be determined by their values at the nodes
$\alpha_{k}~{}(1\leq k\leq M-1)$.

Now we can represent $\hat{\psi}_{h}$ in (4.23) by the basis $\\{g_{h,k,j}\\}$
($1\leq k\leq M-1,~{}0\leq j\leq N-1$)

$\hat{\psi}_{h}(\alpha)=\sum_{k=1}^{M-1}\sum_{j=0}^{N-1}\gamma_{k,j}g_{h,k,j}(\alpha)$
(4.27)

with $\\{\gamma_{k,j}\\}$ the unknown coefficients. As the operator
$\delta\mathscr{F}(\bar{\varphi})$ (4.1) has completely different behavior in
the direction $\bar{\varphi}^{\prime}$ and in the subspace
$\bar{\varphi}^{\prime\perp}$, we will project (4.23) into these two subspaces
respectively. Define $\gamma_{M,j}=\gamma_{0,j}=0$ for $j=0,\cdots,N-1$.
Denote $\gamma^{\perp}\in(\mathbb{R}^{N-1})^{M+1}$ and
$\gamma_{0}\in\mathbb{R}^{M+1}$ with

$\gamma^{\perp}_{k}:=\\{\gamma_{k,j}\\}_{j=1}^{N-1},\qquad(\gamma_{0})_{k}:=\gamma_{k,0}\qquad{\rm
for}~{}k=0,\cdots,M,$ (4.28)

Denote $F^{\perp}\in(\mathbb{R}^{N-1})^{M-1}$ and $F_{0}\in\mathbb{R}^{M-1}$
with

$F^{\perp}_{k}:=\left\\{\xi_{j}(\alpha_{k})^{T}f_{h,k}\right\\}_{j=1}^{N-1},\qquad(F_{0})_{k}:=\xi_{0}(\alpha_{k})^{T}f_{h,k}\qquad{\rm
for}~{}k=1,\cdots,M-1.$ (4.29)

We first consider the problem (4.23) in the subspace perpendicular to the
tangent of MEP
$\xi_{0}=\bar{\varphi}^{\prime}/\lvert\bar{\varphi}^{\prime}\rvert$. By
substituting (4.27) into (4.23), taking $k$-th component, multiplying both
sides with $\xi_{j}(\alpha_{k})^{T},~{}1\leq j\leq N-1$, using (4.16), (4.24)
and (4.26), we obtain by a direct calculation that

$\displaystyle
A(\alpha_{k})\gamma^{\perp}_{k}-\bar{\lambda}(\alpha_{k})(\tilde{D}_{h}\gamma^{\perp})_{k}$
$\displaystyle=F^{\perp}_{k}\qquad{\rm for}~{}k=1,\cdots,M-1,$ (4.30)

where

$A(\alpha):={\rm
Diag}(z_{1}(\alpha),z_{2}(\alpha),\cdots,z_{N-1}(\alpha))=\begin{pmatrix}z_{1}(\alpha)&&&\\\\[4.30554pt]
&z_{2}(\alpha)&&\\\\[4.30554pt] &&\ddots&\\\\[4.30554pt]
&&&z_{N-1}(\alpha)\end{pmatrix}$ (4.31)

for $\alpha\in[0,1]$. Note that the linear system (4.30) is independent of
$\gamma_{0}$.

To better estimate $\|\hat{\psi}_{h}\|_{X}$ by $\|f_{h}\|_{Y_{h}}$, we
transform the linear system (4.30) into another form. We know from (4.6) that
$(\tilde{D}_{h}\gamma^{\perp})_{k_{\bar{s}}}=\frac{1}{2}\big{(}(\tilde{D}_{h}\gamma^{\perp})_{k_{\bar{s}}+1}+(\tilde{D}_{h}\gamma^{\perp})_{k_{\bar{s}}-1}\big{)}$.
Taking $k=k_{\bar{s}}-1,~{}k_{\bar{s}},~{}k_{\bar{s}}+1$ in (4.30), we obtain
the following linear system:

$P_{h,\bar{s}}\begin{pmatrix}(\tilde{D}_{h}\gamma^{\perp})_{k_{\bar{s}}-1}\\\\[4.30554pt]
(\tilde{D}_{h}\gamma^{\perp})_{k_{\bar{s}}+1}\\\\[4.30554pt]
\gamma^{\perp}_{k_{\bar{s}}}\end{pmatrix}=\begin{pmatrix}\frac{F^{\perp}_{k_{\bar{s}}-1}-F^{\perp}_{k_{\bar{s}}}}{-h}\\\\[4.30554pt]
\frac{F^{\perp}_{k_{\bar{s}}+1}-F^{\perp}_{k_{\bar{s}}}}{h}\\\\[4.30554pt]
F^{\perp}_{k_{\bar{s}}}\end{pmatrix},$ (4.32)

where

$\displaystyle
P_{h,\bar{s}}:=\begin{pmatrix}B_{k_{\bar{s}}-1}-\frac{\bar{\lambda}(\alpha_{k_{\bar{s}}})}{2h}I&-\frac{\bar{\lambda}(\alpha_{k_{\bar{s}}})}{2h}I&\frac{A(\alpha_{k_{\bar{s}}})-A(\alpha_{k_{\bar{s}}-1})}{h}\\\\[4.30554pt]
-\frac{\bar{\lambda}(\alpha_{k_{\bar{s}}})}{2h}I&B_{k_{\bar{s}}+1}-\frac{\bar{\lambda}(\alpha_{k_{\bar{s}}})}{2h}I&\frac{A(\alpha_{k_{\bar{s}}+1})-A(\alpha_{k_{\bar{s}}})}{h}\\\\[4.30554pt]
-\frac{\bar{\lambda}(\alpha_{k_{\bar{s}}})}{2}I&-\frac{\bar{\lambda}(\alpha_{k_{\bar{s}}})}{2}I&A(\alpha_{k_{\bar{s}}})\end{pmatrix}\qquad{\rm
with}$ (4.33) $\displaystyle
B_{k_{\bar{s}}-1}:=I+h\frac{A(\alpha_{k_{\bar{s}}-1})}{\bar{\lambda}(\alpha_{k_{\bar{s}}-1})}\qquad{\rm
and}\qquad
B_{k_{\bar{s}}+1}:=I-h\frac{A(\alpha_{k_{\bar{s}}+1})}{\bar{\lambda}(\alpha_{k_{\bar{s}}+1})}.$

For any $\delta\in(0,\frac{\bar{s}}{2})$, there exists a
$k_{\delta}\in\mathbb{Z}_{+}$ such that

$k_{\delta}h\leq\delta<(k_{\delta}+1)h.$ (4.34)

Taking $k=1,\cdots,k_{\bar{s}}-2$ in (4.30) and combining with

$-h\sum_{t=k}^{k_{\bar{s}}-1}(\tilde{D}_{h}\gamma^{\perp})_{t}+\gamma^{\perp}_{k_{\bar{s}}}=\gamma^{\perp}_{k}=h\sum_{t=0}^{k-1}(\tilde{D}_{h}\gamma^{\perp})_{t},\quad
h\sum_{t=0}^{k_{\delta}}(\tilde{D}_{h}\gamma^{\perp})_{t}=\gamma^{\perp}_{k_{\delta}+1},$

we obtain the following linear systems:

$\displaystyle\frac{B_{k}}{\alpha_{k}-\alpha_{k_{\bar{s}}}}(\tilde{D}_{h}\gamma^{\perp})_{k}=\frac{A(\alpha_{k})-A(\alpha_{k_{\bar{s}}})}{\alpha_{k}-\alpha_{k_{\bar{s}}}}\gamma^{\perp}_{k_{\bar{s}}}+\frac{\bar{\lambda}(\alpha_{k_{\bar{s}}})}{\alpha_{k}-\alpha_{k_{\bar{s}}}}(\tilde{D}_{h}\gamma^{\perp})_{k_{\bar{s}}}$
$\displaystyle\hskip
85.35826pt-\frac{F^{\perp}_{k}-F^{\perp}_{k_{\bar{s}}}}{\alpha_{k}-\alpha_{k_{\bar{s}}}}-\frac{hA(\alpha_{k})}{\alpha_{k}-\alpha_{k_{\bar{s}}}}\sum_{t=k+1}^{k_{\bar{s}}-1}(\tilde{D}_{h}\gamma^{\perp})_{t}\qquad{\rm
for}~{}k_{\bar{s}}-k_{\delta}\leq k\leq k_{\bar{s}}-2,$ (4.35)
$\displaystyle\hskip
34.14322ptB_{k}(\tilde{D}_{h}\gamma^{\perp})_{k}=-hA(\alpha_{k})\sum_{t=k+1}^{k_{\bar{s}}-1}(\tilde{D}_{h}\gamma^{\perp})_{t}+A(\alpha_{k})\gamma^{\perp}_{k_{\bar{s}}}-F^{\perp}_{k}$
(4.36)

for $k_{\delta}+1\leq k\leq k_{\bar{s}}-k_{\delta}-1$ and

$P_{h,\delta}\begin{pmatrix}(\tilde{D}_{h}\gamma^{\perp})_{0}\\\\[4.30554pt]
(\tilde{D}_{h}\gamma^{\perp})_{1}\\\\[4.30554pt] \vdots\\\\[4.30554pt]
(\tilde{D}_{h}\gamma^{\perp})_{k_{\delta}}\end{pmatrix}=\begin{pmatrix}\gamma^{\perp}_{k_{\delta}+1}\\\\[4.30554pt]
\frac{F^{\perp}_{1}}{\alpha_{1}}\\\\[4.30554pt] \cdots\\\\[4.30554pt]
\frac{F^{\perp}_{k_{\delta}}}{\alpha_{k_{\delta}}}\end{pmatrix},$ (4.37)

where

$\displaystyle B_{k}$
$\displaystyle:=I+h\frac{A(\alpha_{k})}{\bar{\lambda}(\alpha_{k})}\qquad{\rm
for}~{}1\leq k\leq k_{\bar{s}}-2\hskip 56.9055pt{\rm and}$ (4.38)
$\displaystyle P_{h,\delta}$
$\displaystyle:=\begin{pmatrix}hI&hI&hI&\cdots&hI\\\\[4.30554pt]
\frac{hA(\alpha_{1})}{\alpha_{1}}&-\frac{\bar{\lambda}(\alpha_{1})}{\alpha_{1}}I&\\\\[4.30554pt]
\frac{hA(\alpha_{2})}{\alpha_{2}}&\frac{hA(\alpha_{2})}{\alpha_{2}}&-\frac{\bar{\lambda}(\alpha_{2})}{\alpha_{2}}I&&&\\\\[4.30554pt]
\vdots&\vdots&\vdots&\ddots&&\\\\[4.30554pt]
\frac{hA(\alpha_{k_{\delta}})}{\alpha_{k_{\delta}}}&\frac{hA(\alpha_{k_{\delta}})}{\alpha_{k_{\delta}}}&\frac{hA(\alpha_{k_{\delta}})}{\alpha_{k_{\delta}}}&\cdots&-\frac{\bar{\lambda}(\alpha_{k_{\delta}})}{\alpha_{k_{\delta}}}I\end{pmatrix}.$
(4.39)

Note that we have rewritten (4.30) into (4.32), (4.2), (4.36) and (4.37),
which are used to bound $(\tilde{D}_{h}\gamma^{\perp})_{k}$ at the saddle,
near the saddle, between the saddle and the minimizer $y^{A}_{M}$, and near
the minimizer $y^{A}_{M}$, respectively. Note that for $k>k_{\bar{s}}+1$ we
have the similar linear systems. To estimate
$|\tilde{D}_{h}\gamma^{\perp}|_{\infty}$, we only need to analyze the
coefficient matrix in these linear systems.

We then consider the problem in the tangent direction
$\bar{\varphi}^{\prime}$. By substituting (4.27) into (4.23), taking $k$-th
component, multiplying both sides with $\xi_{0}(\alpha_{k})^{T}$, we obtain

$\displaystyle\alpha_{k}(\alpha_{k}-1)(\alpha_{k}-\bar{s})\left((\tilde{D}_{h}\gamma_{0})_{k}-\sum_{i=0}^{N-1}\sum_{t=1}^{M-1}\gamma_{t,i}c_{h,t,i}\right)=(F_{0})_{k}\qquad{\rm
for}~{}k=1,\cdots,M-1,$ (4.40)

where

$\displaystyle c_{h,t,i}$
$\displaystyle:=\int_{0}^{1}g_{h,t,i}(s)^{T}\xi^{\prime}_{0}(s)~{}{\rm
d}s\quad{\rm for}~{}\alpha\in[0,1].$ (4.41)

We see that $\gamma_{0}$ depends on the component $\gamma^{\perp}$. In our
analysis, we will first focus on the problem (4.30) in the subspace
$\bar{\varphi}^{\prime\perp}$, and then put things back into (4.40) in tangent
direction $\bar{\varphi}^{\prime}$.

So far, we have obtained a spectral representation of (4.23). Now we define a
suitable $\delta$ and turn to the analysis of the coefficient matrics in
(4.32), (4.2), (4.36) and (4.37). Define

$\displaystyle\sigma_{j}:=\frac{z_{j}(0)}{\bar{\lambda}^{\prime}(0)}\qquad{\rm
for}~{}1\leq j\leq N-1\qquad{\rm and}$ (4.42) $\displaystyle
n_{j}(\alpha):=\frac{z_{j}(\alpha)}{\bar{\lambda}(\alpha)}-\frac{\sigma_{j}}{\alpha}\quad\quad{\rm
for}~{}\alpha\in\big{(}0,\frac{\bar{s}}{2}\big{]},~{}1\leq j\leq N-1$ (4.43)

We have from assumption (B) that $\sigma_{j}>1$ for $1\leq j\leq N-1$. Since

$\displaystyle\left\lvert\lim_{\alpha\rightarrow
0^{+}}n_{j}(\alpha)\right\rvert$
$\displaystyle=\left\lvert\lim_{\alpha\rightarrow
1^{-}}\frac{z_{j}(\alpha)\alpha-\sigma_{j}\bar{\lambda}(\alpha)}{\alpha\bar{\lambda}(\alpha)}\right\rvert=\left\lvert\lim_{\alpha\rightarrow
1^{-}}\frac{z^{\prime}_{j}(\alpha)\alpha+z_{j}(\alpha)-\sigma_{j}\bar{\lambda}^{\prime}(\alpha)}{\alpha\bar{\lambda}^{\prime}(\alpha)+\bar{\lambda}(\alpha)}\right\rvert$
$\displaystyle=\left\lvert\lim_{\alpha\rightarrow
0^{+}}\frac{2z^{\prime}_{j}(\alpha)+z^{\prime\prime}_{j}(\alpha)\alpha-\sigma_{j}\bar{\lambda}^{\prime\prime}(\alpha)}{\alpha\bar{\lambda}^{\prime\prime}(\alpha)+2\bar{\lambda}^{\prime}(\alpha)}\right\rvert=\left\lvert\frac{2z^{\prime}_{j}(0)-\sigma_{j}\bar{\lambda}^{\prime\prime}(0)}{2\bar{\lambda}^{\prime}(0)}\right\rvert<C$
(4.44)

with $C$ a constant depending only on $\bar{\lambda}$ and $z_{j}$, we have
that $n_{j}$ is bounded on $(0,\frac{\bar{s}}{2}]$. Define

$\bar{\sigma}:=\frac{1}{2}\left(1+\max_{1\leq j\leq N-1}\sigma_{j}\right).$
(4.45)

From the fact that $\sigma_{j}>1$ for $1\leq j\leq N-1$ , we have
$\bar{\sigma}>1$. Therefore, there exists $\delta\in(0,\frac{\bar{s}}{2})$
such that $z_{j}(\alpha)>0$ for
$\alpha\in[\bar{s}-\delta-2h,\bar{s}+h],~{}1\leq j\leq N-1$,
$\frac{1}{2}\bar{\lambda}^{\prime}(0)\leq\lVert\bar{\lambda}^{\prime}\rVert_{C([0,\delta];\mathbb{R})}\leq
2\bar{\lambda}^{\prime}(0)$ and

$\frac{\sigma_{j}-\bar{\sigma}}{\delta}\geq\|n_{j}\|_{C((0,\frac{\bar{s}}{2}];\mathbb{R})}\qquad{\rm
for}~{}1\leq j\leq N-1.$ (4.46)

This together with (4.43) implies

$\frac{z_{j}(\alpha)}{\bar{\lambda}(\alpha)}=\frac{\sigma_{j}}{\alpha}+n_{j}(\alpha)\geq\frac{\bar{\sigma}}{\alpha}\qquad\forall~{}\alpha\in(0,\delta].$
(4.47)

######  Lemma 4.3.

Under the assumption of Lemma 4.2, let $\bar{\sigma}$, $\delta$, $k_{\delta}$
and $B_{k}$’ be given by (4.45), (4.46), (4.34) and (4.38), respectively. Then
for sufficient small $h$, $B_{k}$’s are invertible and there exist a constant
$C>0$ depending only on $\bar{\lambda}$ and $A$ such that

$\displaystyle\lVert
B_{k}^{-1}\rVert_{\infty}\leq\frac{\bar{\sigma}}{k+\bar{\sigma}}\qquad$
$\displaystyle{\rm for}~{}1\leq k\leq k_{\delta},$ (4.48) $\displaystyle\lVert
B_{k}^{-1}\rVert_{\infty}\leq 1+Ch\qquad$ $\displaystyle{\rm
for}~{}k_{\delta}+1\leq k\leq k_{\bar{s}}-k_{\delta}-1,$ (4.49)
$\displaystyle\lVert
B_{k}^{-1}\rVert_{\infty}<\frac{1}{\alpha_{k_{\bar{s}}}-\alpha_{k}}\qquad$
$\displaystyle{\rm for}~{}k_{\bar{s}}-k_{\delta}\leq k\leq k_{\bar{s}}-2,$
(4.50)

where the constants $C$ depend only on $\bar{\lambda}$ and $A$. Moreover, we
have

$\lVert B_{k}^{-1}-I\rVert_{\infty}\leq\frac{C}{k}\qquad{\rm for}~{}1\leq
k\leq k_{\delta}.$ (4.51)

###### Proof.

Since $B_{k}$ is a diagonal matrix, we have

$\lVert B_{k}^{-1}\rVert_{\infty}=\max_{1\leq j\leq
N-1}\left|\frac{1}{1+\frac{z_{j}(\alpha_{k})h}{\bar{\lambda}(\alpha_{k})}}\right|\qquad{\rm
and}\qquad\lVert B_{k}^{-1}-I\rVert_{\infty}=\max_{1\leq j\leq
N-1}\left|\frac{\frac{z_{j}(\alpha_{k})h}{\bar{\lambda}(\alpha_{k})}}{1+\frac{z_{j}(\alpha_{k})h}{\bar{\lambda}(\alpha_{k})}}\right|$
(4.52)

for $1\leq k\leq k_{\bar{s}}-1$. For $1\leq k\leq k_{\delta}$, using (4.45)
and (4.47) yields

$\displaystyle 1+\frac{z_{j}(\alpha_{k})}{\bar{\lambda}(\alpha_{k})}h\geq
1+\frac{\bar{\sigma}h}{\alpha_{k}}\geq 1+\frac{\bar{\sigma}}{k}\qquad{\rm
for}~{}1\leq j\leq N-1.$

This together with (4.52) leads to

$\displaystyle\lVert
B_{k}^{-1}\rVert_{\infty}\leq\frac{k}{k+\bar{\sigma}}\qquad{\rm
and}\qquad\lVert
B_{k}^{-1}-I\rVert_{\infty}\leq\frac{Ch}{\bar{\lambda}(\alpha_{k})}\leq\frac{C}{k}\frac{\alpha_{k}}{\int_{0}^{\alpha_{k}}\lvert\bar{\lambda}^{\prime}(s)\rvert~{}{\rm
d}s}\leq\frac{C}{k},$

where $C$ depends on $\bar{\lambda}$ and $A$. This completes the proof for
(4.48) and (4.51).

For $k_{\delta}+1\leq k\leq k_{\bar{s}}-k_{\delta}-1$, since
$\lvert\bar{\lambda}(\alpha_{k})\rvert$ have both upper bound and lower bound,
we have from (4.52) that

$\displaystyle\lVert B_{k}^{-1}\rVert_{\infty}\leq\frac{1}{1-Ch}\leq 1+Ch,$

where $C$ depends on $\bar{\lambda}$ and $A$. This completes the proof for
(4.49).

For $k_{\bar{s}}-k_{\delta}\leq k\leq k_{\bar{s}}-2$, using the fact that
$z_{j}(\alpha_{k})>0$ and $\bar{\lambda}(\alpha_{k})>0$, we obtain

$\left\lVert B_{k}^{-1}\right\rVert_{\infty}=\max_{1\leq j\leq
N-1}\frac{1}{1+\frac{z_{j}(\alpha_{k})h}{\bar{\lambda}(\alpha_{k})}}<1<\frac{1}{\alpha_{k_{\bar{s}}}-\alpha_{k}}.\qed$

######  Lemma 4.4.

Under the assumption of Lemma 4.2, let $\delta$, $k_{\delta}$, $P_{h,\bar{s}}$
and $P_{h,\delta}$ be given by (4.46), (4.34), (4.33) and (4.39),
respectively. Then for sufficiently small $h$, $P_{h,\bar{s}}$ and
$P_{h,\delta}$ are invertible and there exist constants $C>0$ depending on
$\bar{\lambda}$ and $A$ such that

$\lVert P_{h,\bar{s}}^{-1}\rVert_{\infty}\leq C\qquad{\rm and}\qquad\lVert
P_{h,\delta}^{-1}\rVert_{\infty}\leq C.$ (4.53)

###### Proof.

Since $\displaystyle\lim_{h\rightarrow 0^{+}}\alpha_{k_{\bar{s}}}=\bar{s}$ and
$\bar{\lambda}(\bar{s})=0$, we have

$\displaystyle P_{\bar{s}}:=\lim_{h\rightarrow
0^{+}}P_{h,\bar{s}}=\begin{pmatrix}B_{\bar{s}}&&A^{\prime}(\bar{s})\\\\[4.30554pt]
&B_{\bar{s}}&A^{\prime}(\bar{s})\\\\[4.30554pt]
&&A(\bar{s})\end{pmatrix},\qquad
P_{\bar{s}}^{-1}=\begin{pmatrix}B_{\bar{s}}^{-1}&&-B_{\bar{s}}^{-1}A^{\prime}(\bar{s})A(\bar{s})^{-1}\\\\[4.30554pt]
&B_{\bar{s}}^{-1}&-B_{\bar{s}}^{-1}A^{\prime}(\bar{s})A(\bar{s})^{-1}\\\\[4.30554pt]
&&A(\bar{s})^{-1}\end{pmatrix},$

where
$B_{\bar{s}}:=I-\frac{A(\bar{s})}{\bar{\lambda}^{\prime}(\bar{s})}\in\mathbb{R}^{(N-1)\times(N-1)}$.
From Lemma 4.1 (ii) and the fact that $z_{j}(\bar{s})>0$ for $1\leq j\leq
N-1$, we obtain

$\lVert B_{\bar{s}}^{-1}\rVert_{\infty}\leq\max_{1\leq j\leq
N-1}\frac{1}{1-\frac{z_{j}(\bar{s})}{\bar{\lambda}^{\prime}(\bar{s})}}\leq
1,\qquad\lVert A(\bar{s})^{-1}\rVert_{\infty}\leq\max_{1\leq j\leq
N-1}\frac{1}{z_{j}(\bar{s})}.$

This implies

$\lVert P_{\bar{s}}^{-1}\rVert_{\infty}\leq\lVert
B_{\bar{s}}^{-1}\rVert_{\infty}+\lVert B_{\bar{s}}^{-1}\rVert_{\infty}\lVert
A^{\prime}(\bar{s})\rVert_{\infty}\lVert A(\bar{s})^{-1}\rVert_{\infty}+\lVert
A(\bar{s})^{-1}\rVert_{\infty}\leq C,$

where $C$ depends only on $\bar{\lambda}$ and $A$. Thus $P_{h,\bar{s}}$ is
inveritible and $\lVert P_{h,\bar{s}}\rVert_{\infty}\leq C$ for sufficiently
small $h$. This completes the proof for first part of (4.53).

We then investigate the inverse of $P_{h,\delta}$:

$P_{h,\delta}^{-1}=Q_{h,\delta}T_{h,\delta},$ (4.54)

where

$\displaystyle Q_{h,\delta}$
$\displaystyle:=\begin{pmatrix}h^{-1}B_{1}^{-1}\cdots
B_{k_{\delta}}^{-1}&-B_{1}^{-1}&-B_{1}^{-1}B_{2}^{-1}&\cdots&-B_{1}^{-1}\cdots
B_{k_{\delta}}^{-1}\\\\[4.30554pt] h^{-1}(I-B_{1}^{-1})B_{2}^{-1}\cdots
B_{k_{\delta}}^{-1}&B_{1}^{-1}&(B_{1}^{-1}-I)B_{2}^{-1}&\cdots&(B_{1}^{-1}-I)B_{2}^{-1}\cdots
B_{k_{\delta}}^{-1}\\\\[4.30554pt] h^{-1}(I-B_{2}^{-1})B_{3}^{-1}\cdots
B_{k_{\delta}}^{-1}&&B_{2}^{-1}&\cdots&(B_{2}^{-1}-I)B_{3}^{-1}\cdots
B_{k_{\delta}}^{-1}\\\\[4.30554pt] \vdots&&&\ddots&\\\\[4.30554pt]
h^{-1}(I-B_{k_{\delta}}^{-1})&&&&B_{k_{\delta}}^{-1}\end{pmatrix},$
$\displaystyle T_{h,\delta}$
$\displaystyle:=\begin{pmatrix}I&&&\\\\[4.30554pt]
&-\frac{\alpha_{1}}{\bar{\lambda}(\alpha_{1})}I&&\\\\[4.30554pt]
&&\ddots&\\\\[4.30554pt]
&&&-\frac{\alpha_{k_{\delta}}}{\bar{\lambda}(\alpha_{k_{\delta}})}I\end{pmatrix}.$

We see from Lemma 4.1 (iii) that $\lVert T_{h,\delta}\rVert_{\infty}\leq C$
where $C$ only depends on $\bar{\lambda}$ and $\delta$. To estimate $\lVert
Q_{h,\delta}\rVert_{\infty}$, we first present two inequalities. We claim that
there exist constants $C_{1},C_{2},C>0$ depending only on $\bar{\sigma}$ such
that

$\displaystyle
C_{1}\frac{1}{n^{\bar{\sigma}}}\leq\prod_{s=1}^{n}\frac{s}{s+\bar{\sigma}}\leq
C_{2}\frac{1}{n^{\bar{\sigma}}}\qquad$ $\displaystyle{\rm
for}~{}n\in\mathbb{Z}_{+}\qquad{\rm and}$ (4.55)
$\displaystyle\frac{1}{n^{\bar{\sigma}}}\leq
C\left(\frac{1}{(n-1)^{\bar{\sigma}-1}}-\frac{1}{n^{\bar{\sigma}-1}}\right)\qquad$
$\displaystyle{\rm for}~{}n\geq 2,$ (4.56)

The inequality (4.55) results from induction and the fact that

$\lim_{n\rightarrow+\infty}\frac{n+1}{n+1+\bar{\sigma}}=1=\lim_{n\rightarrow+\infty}\frac{n^{\bar{\sigma}}}{(n+1)^{\bar{\sigma}}}.$

Using the monotonicity, there exists a constant $C>0$ depending on
$\bar{\sigma}$ such that

$\frac{1}{n}\leq\frac{1}{2}\leq C(2^{\bar{\sigma}-1}-1)\leq
C\left(\frac{n^{\bar{\sigma}-1}}{(n-1)^{\bar{\sigma}-1}}-1\right)\qquad\forall~{}n\geq
2,$

which leads to (4.56) by multiplying $n^{1-\bar{\sigma}}$ on both sides.

Applying (4.55) and Lemma 4.3, we first estimate

$\left\lVert\prod_{t=k+1}^{k^{\prime}}B_{t}^{-1}\right\rVert_{\infty}\leq\prod_{t=k+1}^{k^{\prime}}\frac{t}{t+\bar{\sigma}}\leq
C\frac{(k+1)^{\bar{\sigma}}}{(k^{\prime})^{\bar{\sigma}}}\leq
C(k+1)^{\bar{\sigma}}\left(\frac{1}{(k^{\prime}-1)^{\bar{\sigma}-1}}-\frac{1}{(k^{\prime})^{\bar{\sigma}-1}}\right)$
(4.57)

for $1\leq k<k^{\prime}\leq k_{\delta}$. Since $k_{\delta}$ scales as
$O(h^{-1})$, using (4.57) and Lemma 4.3, we obtain

$\displaystyle\lVert h^{-1}B_{1}^{-1}\cdots
B_{k_{\delta}}^{-1}\rVert_{\infty}\leq\frac{C}{h(k_{\delta})^{\bar{\sigma}}}\leq
C,$ (4.58)
$\displaystyle\sum_{k=2}^{k_{\delta}}\left\lVert\prod_{t=1}^{k}B_{t}^{-1}\right\rVert_{\infty}\leq
C\sum_{k=2}^{k_{\delta}}\left(\frac{1}{(k-1)^{\bar{\sigma}-1}}-\frac{1}{k^{\bar{\sigma}-1}}\right)\leq
C\left(1-\frac{1}{k_{\delta}^{\bar{\sigma}-1}}\right)\leq C,$ (4.59)

and

$\displaystyle\sum_{t=k+1}^{k_{\delta}}\lVert(B_{k}^{-1}-I)\prod_{t^{\prime}=k+1}^{t}B_{t^{\prime}}^{-1}\rVert_{\infty}$
$\displaystyle\leq\frac{C(k+1)^{\bar{\sigma}}}{k}\sum_{t=k+1}^{k_{\delta}}\left(\frac{1}{(t-1)^{\bar{\sigma}-1}}-\frac{1}{t^{\bar{\sigma}-1}}\right)$
$\displaystyle\leq\frac{C(k+1)}{k}\left(\frac{(k+1)^{\bar{\sigma}-1}}{k^{\bar{\sigma}-1}}-\frac{(k+1)^{\bar{\sigma}-1}}{k_{\delta}^{\bar{\sigma}-1}}\right)\leq
C,$ (4.60) $\displaystyle\lVert
h^{-1}(I-B_{k}^{-1})\prod_{t=k+1}^{k_{\delta}}B_{t}^{-1}\rVert_{\infty}$
$\displaystyle\leq\frac{C}{kh}\frac{(k+1)^{\bar{\sigma}}}{(k_{\delta})^{\bar{\sigma}}}\leq
C$ (4.61)

for $1\leq k\leq k_{\delta}-1$. Combining (4.58)-(4.61), we have from (4.54)
that

$\lVert P_{h,\delta}^{-1}\rVert_{\infty}\leq\lVert
Q_{h,\delta}\rVert_{\infty}\lVert T_{h,\delta}\rVert_{\infty}\leq C,$

where $C$ only depends on $\bar{\lambda}$, $\delta$ and $\bar{\sigma}$. This
completes the proof of (4.53). ∎

Before proving Lemma 4.2, we give the discretization error of the operators
$\tilde{D}_{h},~{}\hat{D}_{h}$ and $L_{h}$. Given $\varphi\in
C^{2}\big{(}[0,1];\mathbb{R}^{N}\big{)}\cap\mathscr{A}$, we see from (4.6) and
(4.9) that there exists a constant $C>0$ independent of $h$, such that

$\displaystyle|\tilde{D}_{h}\Pi_{h}\varphi-\Pi_{h}\varphi^{\prime}|_{\infty}+|\hat{D}_{h}\Pi_{h}\varphi-\Pi_{h}\varphi^{\prime}|_{\infty}\leq
Ch.$ (4.62)

This together with the standard error estimate of the numerical integration
leads to

$|L_{h}(\Pi_{h}\varphi)-L(\varphi)|\leq\Big{|}h\sum_{t=1}^{M}\big{(}|(\hat{D}_{h}\Pi_{h}\varphi)_{t}|-|\varphi^{\prime}(\alpha_{t})|\big{)}\Big{|}+\Big{|}h\sum_{t=1}^{M}\varphi^{\prime}(\alpha_{t})|-\int_{0}^{1}|\varphi^{\prime}(s)|~{}{\rm
d}s\Big{|}\leq Ch.$ (4.63)

We obtain by a similar argument that

$\displaystyle|\delta L_{h}(\Pi_{h}\varphi)\psi_{h}-\delta
L(\varphi)(P_{h}\psi_{h})|\leq Ch\qquad\forall~{}\psi_{h}\in X_{h},$ (4.64)

where $P_{h}$ is defined in (4.22).

Now we are ready to show the stability of the discrete MEP.

###### Proof of Lemma 4.2.

For any $\psi_{h}\in X_{h}$,we have

$\delta\mathscr{F}_{h}(\Pi_{h}\bar{\varphi})\psi_{h}=\big{(}\delta{\mathscr{F}}_{h}(\Pi_{h}\bar{\varphi})\psi_{h}-\Pi_{h}\delta\mathscr{F}(\bar{\varphi})(P_{h}\psi_{h})\big{)}+\Pi_{h}\delta\mathscr{F}(\bar{\varphi})(P_{h}\psi_{h}).$
(4.65)

By comparing $\delta\mathscr{F}_{h}(\Pi_{h}\bar{\varphi})\psi_{h}$ in (4.1)
and $\delta\mathscr{F}(\bar{\varphi})(P_{h}\psi_{h})$ in (4.1), using (4.62),
(4.63) and (4.64), we obtain

$\lVert\delta{\mathscr{F}}_{h}(\Pi_{h}\bar{\varphi})\psi_{h}-\Pi_{h}\delta\mathscr{F}(\bar{\varphi})(P_{h}\psi_{h})\rVert_{Y_{h}}\leq
Ch\lVert P_{h}\psi_{h}\rVert_{X}\leq Ch\lVert\psi_{h}\rVert_{X_{h}},$ (4.66)

where $C$ depends only on $E$ and $\bar{\varphi}$.

Now we turn to the second term in (4.65). For any $\psi_{h}\in X_{h}$, let
$\hat{\psi}_{h}=P_{h}\psi_{h}$ and
$f_{h}=\Pi_{h}\delta\mathscr{F}(\bar{\varphi})\hat{\psi}_{h}$. With the
representation (4.27), using Lemma 4.3 and Lemma 4.4, we have from (4.32),
(4.2), (4.36) and (4.37) that

$|(\tilde{D}_{h}\gamma^{\perp})_{k}|_{\infty}\leq C\lVert
F^{\perp}\rVert_{Y_{h}}\qquad{\rm for}~{}0\leq k\leq k_{\bar{s}}+1,$

where $\gamma^{\perp}$ and $F^{\perp}$ are given by (4.28) and (4.29),
respectively. A similar estimate holds for $k_{\bar{s}}+2\leq k\leq M$ and
then we have

$|\tilde{D}_{h}\gamma^{\perp}|_{\infty}\leq C\lVert
F^{\perp}\rVert_{Y_{h}}\leq C\lVert f_{h}\rVert_{Y_{h}},$ (4.67)

where $C$ depends only on $\bar{\lambda}$ and $A$. We obtain by taking
$k=1,\cdots,k_{\bar{s}}-1,k_{\bar{s}}+1,\cdots,M-1$ in (4.40) that

$|\tilde{D}_{h}\gamma_{0}|_{\infty}\leq\|F_{0}\|_{Y_{h}}+|\gamma^{\perp}|_{\infty}\sum_{t=1}^{M-1}\sum_{j=1}^{N-1}|c_{h,t,j}|+|\gamma_{0}|_{\infty}\sum_{t=1}^{M-1}|c_{h,t,0}|,$
(4.68)

where $\gamma_{0}$, $F_{0}$ and $c_{h,t,i}$ are given by (4.28), (4.29) and
(4.41), respectively. From the definition of $g_{h,t,i}$ in (4.26) and
$\hat{X}_{h}$ in (4.21), we know that supp
$(g_{h,t,i})\subset[\alpha_{t-2},\alpha_{t+2}]$ for $1\leq t\leq M-1,~{}0\leq
i\leq N-1$. Since $\xi_{0}\in C^{2}\big{(}[0,1];\mathbb{R}^{N}\big{)}$ and
$|\xi_{0}(\alpha)|^{2}\equiv 1$, we have

$\displaystyle|c_{h,t,0}|$ $\displaystyle\leq
4h\|\xi_{0}(\alpha_{t})^{T}\xi^{\prime}_{0}(\alpha)\|_{C([\alpha_{t-2},\alpha_{t+2}];\mathbb{R})}\leq
Ch^{2}\qquad{\rm and}$ (4.69) $\displaystyle|c_{h,t,j}|$ $\displaystyle\leq
4h\left\|g_{h,t,j}^{T}\xi^{\prime}_{0}\right\|_{C([\alpha_{t-2},\alpha_{t+2}];\mathbb{R})}\leq
Ch,$ (4.70)

where $C$ depends only on $\\{\xi_{i}\\}_{i=0}^{N-1}$. Combining
(4.67)-(4.70), we have for sufficiently small $h$ that

$\lvert\tilde{D}_{h}\gamma_{0}\rvert_{\infty}\leq
C\lvert\gamma^{\perp}\rvert_{\infty}+C\lVert F_{0}\rVert_{Y_{h}}\leq C\lVert
f_{h}\rVert_{Y_{h}}.$ (4.71)

From (4.27) we have

$\displaystyle(\hat{D}_{h}\psi_{h})_{k}$
$\displaystyle=h^{-1}\sum_{i=0}^{N-1}\big{(}\gamma_{k,i}\xi_{i}(\alpha_{k})-\gamma_{k-1,i}\xi_{i}(\alpha_{k-1})\big{)}$
$\displaystyle=\sum_{i=0}^{N-1}\left(\frac{\gamma_{k,i}-\gamma_{k-1,i}}{h}\xi_{i}(\alpha_{k})+\gamma_{k-1,i}\frac{\xi_{i}(\alpha_{k})-\xi_{i}(\alpha_{k-1})}{h}\right)\qquad{\rm
for}~{}1\leq k\leq M.$

This together with (4.67) and (4.71) implies that for sufficient small $h$,
there exists a $C>0$ depending only on $\bar{\lambda},~{}A$ and
$\\{\xi_{i}\\}_{i=0}^{N-1}$ such that

$\lVert\psi_{h}\rVert_{X_{h}}\leq
C\left|\begin{pmatrix}\tilde{D}_{h}\gamma^{\perp}\\\\[4.30554pt]
\tilde{D}_{h}\gamma_{0}\end{pmatrix}\right|_{\infty}\leq C\lVert
f_{h}\rVert_{Y_{h}}.$

This indicates

$\|\psi_{h}\|_{X_{h}}\leq
C\|\Pi_{h}\delta\mathscr{F}(\bar{\varphi})(P_{h}\psi_{h})\|_{Y_{h}}\qquad\forall~{}\psi_{h}\in
X_{h}.$ (4.72)

Combining (4.65), (4.66) and (4.72) yields our stated result. ∎

###  4.3 Convergence analysis

Now we are ready to prove the convergence and convergence rate of the
approximation $\bar{\varphi}_{h}$ with respect to the mesh size $h$.

######  Theorem 4.1.

Let $\bar{\varphi}\in C^{3}\big{(}[0,1];\mathbb{R}^{N}\big{)}$ be a solution
of (2.1). Assume that (A) and (B) are satisfied. Then for sufficiently small
$h$, there exist a solution $\bar{\varphi}_{h}\in\mathscr{A}_{h}$ of (4.13)
and a constant $C$ depending only on $\bar{\varphi}$ and $E$, such that

$\displaystyle\lVert\bar{\varphi}_{h}-\Pi_{h}\bar{\varphi}\rVert_{X_{h}}\leq
Ch.$

###### Proof.

The proof is divided into three parts. We will first show the consistency and
stability of the discrete MEP equation (4.13), and then derive the convergence
by inverse function theorem.

1\. Consistency. By comparing (4.11) and (4.1), we have from (4.62) and (4.63)
that

$\displaystyle\lVert\mathscr{F}_{h}(\Pi_{h}\bar{\varphi})\rVert_{Y_{h}}=\lVert{\mathscr{F}}_{h}(\Pi_{h}\bar{\varphi})-\Pi_{h}\mathscr{F}(\bar{\varphi})\rVert_{Y_{h}}\leq
Ch,$ (4.73)

where $C$ only depends on $E$ and $\bar{\varphi}$.

2\. Stability. It follows directly from Lemma 4.2 that there exists a constant
$C$ depending on $E$ and $\bar{\varphi}$ such that

$\lVert\delta\mathscr{F}_{h}(\Pi_{h}\bar{\varphi})\psi_{h}\rVert_{Y_{h}}\geq
C\lVert\psi_{h}\rVert_{X_{h}}\qquad\forall~{}\psi_{h}\in X_{h}.$ (4.74)

3\. Application of the inverse function theory. With the consistency (4.73)
and stability (4.74), we can apply the inverse function theorem [13, Lemma
2.2] to obtain that, for $h$ sufficiently small, the solution
$\bar{\varphi}_{h}$ of (4.13) exists and

$\lVert\bar{\varphi}_{h}-\Pi_{h}\bar{\varphi}\rVert_{X_{h}}\leq Ch.\qed$

Finally, we can prove Theorem 2.1 by building a connection between the
equations (2.8) and (4.13).

###### Proof of Theorem 2.1.

Let $\bar{\varphi}_{h}$ be a solution of (4.13). Note that the saddle $y_{S}$
reaches the energy maximum along the MEP. Then we have from the convergence
result in Theorem 4.1 that for sufficiently small $h$,

$E(\bar{\varphi}_{h,0})<E(\bar{\varphi}_{h,1})<\cdots<E(\bar{\varphi}_{h,k_{\bar{s}}})>E(\bar{\varphi}_{h,k_{\bar{s}}+1})>\cdots>E(\bar{\varphi}_{h,M}).$

By comparing (2.2) and (4.6), we see that these two first order difference
schemes of $\bar{\varphi}_{h}$ are exactly the same, i.e.,
$\tilde{D}_{h}\bar{\varphi}_{h}=D_{h}\bar{\varphi}_{h}$. This together with
(4.14) yields that $\bar{\varphi}_{h}$ is a solution of (2.8) and it remains
to show the estimates (2.9) and (2.10).

From (4.62) and Theorem 4.1, we estimate

$\displaystyle\max_{0\leq k\leq
M}\Big{(}|(D_{h}\bar{\varphi}_{h})_{k}-\bar{\varphi}^{\prime}(\alpha_{k})|+|\bar{\varphi}_{h,k}-\bar{\varphi}(\alpha)|\Big{)}\leq
C\big{(}|\tilde{D}_{h}\bar{\varphi}_{h}-\Pi_{h}\bar{\varphi}^{\prime}|_{\infty}+|\bar{\varphi}_{h}-\Pi_{h}\bar{\varphi}|_{\infty}\big{)}$
$\displaystyle\hskip 85.35826pt\leq
C|\tilde{D}_{h}(\bar{\varphi}_{h}-\Pi_{h}\bar{\varphi})|_{\infty}+C|\tilde{D}_{h}\Pi_{h}\bar{\varphi}-\Pi_{h}\bar{\varphi}^{\prime}|_{\infty}+C|\bar{\varphi}_{h}-\Pi_{h}\bar{\varphi}|_{\infty}$
$\displaystyle\hskip 85.35826pt\leq
C|\hat{D}_{h}(\bar{\varphi}_{h}-\Pi_{h}\bar{\varphi})|_{\infty}+Ch+C|\bar{\varphi}_{h}-\Pi_{h}\bar{\varphi}|_{\infty}\leq
Ch,$

which completes the proof of (2.9). Then we turn to the error estimate for the
energy barrier. Using the fact that $\nabla
E\big{(}\bar{\varphi}(\bar{s})\big{)}=0$ and $E\in C^{4}(\mathbb{R}^{N})$, we
have

$\big{|}\delta E(\bar{\varphi})-\delta
E(\bar{\varphi}_{h})\big{|}=\big{|}E\big{(}\bar{\varphi}(\bar{s})\big{)}-E(\bar{\varphi}_{h,k_{\bar{s}}})\big{|}\leq
C|\bar{\varphi}(\bar{s})-\bar{\varphi}_{h,k_{\bar{s}}})|_{\infty}^{2},$ (4.75)

where $C$ is a constant depending only on $E$ and $\bar{\varphi}$. Theorem 4.1
and the fact that $\bar{\varphi}\in C^{3}\big{(}[0,1];\mathbb{R}^{N}\big{)}$
implies

$|\bar{\varphi}(\bar{s})-\bar{\varphi}_{h,k_{\bar{s}}})|_{\infty}\leq|\bar{\varphi}(\bar{s})-\bar{\varphi}(\alpha_{k_{\bar{s}}})|_{\infty}+|\bar{\varphi}(\alpha_{k_{\bar{s}}})-\bar{\varphi}_{h,k_{\bar{s}}})|_{\infty}\leq
Ch,$

which together with (4.75) yields (2.10) and completes the proof. ∎

## References

  * [1] N. Berglund.  Kramers’ law: Validity, derivations and generalisations.  Markov Process. Related Fields, 19:459–490, 2013.
  * [2] J. Braun and C. Ortner.  Sharp uniform convergence rate of the supercell approximation of a crystalline defect.  SIAM J. Numer. Anal., 58:279–297, 2020.
  * [3] M. Cameron, R. Kohn, and E. Vanden-Eijnden.  The string method as a dynamical system.  J. Nonlinear Sci., 21:193–230, 2011.
  * [4] W. E, W. Ren, and E. Vanden-Eijnden.  String method for the study of rare events.  Phys. Rev. B, 66:052301, 2002.
  * [5] W. E, W. Ren, and E. Vanden-Eijnden.  Simplified and improved string method for computing the minimum energy paths in barrier-crossing events.  J. Chem. Phys., 126:164103, 2007.
  * [6] P. Hänggi, P. Talkner, and M. Borkovec.  Reaction-rate theory: fifty years after kramers.  Rev. Modern Phys., 62:251–341, 1990.
  * [7] G. Henkelman, B. P. Uberuaga, and H. Jónsson.  A climbing image nudged elastic band method for finding saddle points and minimum energy paths.  J. Chem. Phys., 113:9901–9904, 2000.
  * [8] H. Jónsson, G. Mills, and K. W. Jacobsen.  Nudged elastic band method for finding minimum energy paths of transitions.  In Classical and Quantum Dynamics in Condensed Phase Simulations, pages 385–404. Citeseer, 1998.
  * [9] B. Koten and M. Luskin.  Stability and convergence of the string method for computing minimum energy paths.  Multiscale Model. Simul., 17:873–898, 2018.
  * [10] G. D. Leines and J. Rogal.  Comparison of minimum-action and steepest-descent paths in gradient systems.  Phys. Rev. E, 93:022307, 2016.
  * [11] X. Liu, H. Chen, and C. Ortner.  Stability of the minimum energy path.  ArXiv e-prints, 2204.00984, 2022.
  * [12] K. Müller and L. D. Brown.  Location of saddle points and minimum energy paths by a constrained simplex optimization procedure.  Theor. Chim. Acta, 53:75–93, 1979.
  * [13] C. Ortner.  A priori and a posteriori analysis of the quasinonlocal quasicontinuum method in 1d.  Math. Comp., 80:1265–1285, 2011.
  * [14] C. Ortner, J. Kermode, et al.  JuLIP: Julia Library for Interatomic Potentials.  https://github.com/JuliaMolSim/JuLIP.jl.
  * [15] W. Ren and E. Vanden-Eijnden.  A climbing string method for saddle point search.  J. Chem. Phys., 138:134105, 2013.
  * [16] G. H. Vineyard.  Frequency factors and isotope effects in solid state rate processes.  J. Phys. Chem. Solids, 3:121–127, 1957.